The growing importance of subsurface carbon storage for tackling anthropogenic carbon emissions requires new ideas to improve the rate and cost of carbon capture and storage (CCS) project development and implementation. We assessed sandstones from the UK Geoenergy Observatories (UKGEOS) site in Glasgow, UK and the Wilmslow Sandstone Formation (WSF) in Cumbria, UK at the pore scale to indicate suitability for further assessment as CCS reservoirs. We measured porosity, permeability and other pore geometry characteristics using digital rock physics techniques on microcomputed tomographic images of core material from each site. We found the Glasgow material to be unsuitable for CCS due to very low porosity (up to 1.65%), whereas the WSF material showed connected porosity up to 26.3% and permeabilities up to 6040 mD. Our results support the presence of a percolation threshold at 10% total porosity, introducing near full connectivity. We found total porosity varies with permeability with an exponent of 3.19. This provides a reason to assume near full connectivity in sedimentary samples showing porosities above this threshold without the need for expensive and time-consuming analyses.

Supplementary material: Information about the boreholes sampled in this study, additional well logs of boreholes and a summary of the supporting data plotted throughout this article from literature are available at https://doi.org/10.6084/m9.figshare.c.5260074

Thematic collection: This article is part of the Geoscience for CO2 storage collection available at: https://www.lyellcollection.org/cc/geoscience-for-co2-storage

The changing global attitude towards greenhouse gas emissions and their environmental impact has led to investment in innovative methods for reducing anthropogenic CO2 release. One such method is carbon capture and storage (CCS), the process by which CO2 captured from point sources is injected into the subsurface and sequestered for geologically significant periods of time. Following injection, CO2 may be retained on the short term as a fluid either trapped by an impermeable cap rock (Johnson et al. 2004; Song and Zhang 2013), as isolated ganglia through capillary trapping (Blunt et al. 2013; Boot-Handford et al. 2014) or dissolved in the formation water (Ghesmat et al. 2011; Unwin et al. 2016). On the longer timescale, carbonate minerals such as calcite, dawsonite and ankerite form (Ghesmat et al. 2011; Jiang and Tsuji 2014; Zhang and Song 2014), locking carbon away in a thermodynamically stable state through geological carbon storage (GCS).

Global investment in CCS technologies and projects is ever increasing. For example, the EU Commission Innovation Fund plans to spend €10 billion up to 2030 on CCS-related projects in Europe alone (Page et al. 2019) and the Norwegian government have pledged US $1.8 billion for the Longship CCS project (IEA 2020). Such investments will continue to increase following environmental and climate pledges, such as the UK's, to become carbon neutral by 2050 as part of the June 2019 amendment to the Climate Change Act (Government of the United Kingdom 2008). The UK government recognizes CCS as an important component of its decarbonization, and, consequently, announced £800 million in funding for carbon capture utilization and storage (CCUS) projects (IEA 2020) and commissioned research into viable CCUS solutions to be employed at a large scale in the 2030s (CCUS Cost Challenge Taskforce 2018).

One project that provides the opportunity to test some of these concepts on a large scale is the UK Geoenergy Observatories (UKGEOS) project. This project aims to characterize the subsurface and to assess feasibility for CCS at two key sites in Glasgow and Cheshire (Kingdon et al. 2018; Monaghan et al. 2019). The Cheshire site will play a key role in determining the quality of storage reservoirs that could support the proposed Merseyside industrial cluster (CCUS Cost Challenge Taskforce 2018). Research sites, such as the UKGEOS, are important for supporting the development of technologies that the energy industry is increasing investment in. CCS is an important method for industrial offsetting of CO2 emissions, required due to taxes on greenhouse gas emissions to fulfill commitments to global climate agreements (Government of Norway 2017; Mardones and Flores 2018). Similarly, the USA have directly incentivized CCS through implementing an expansion of the 45Q tax credits scheme which puts a value on sequestering CO2 (IEA 2020). Further progress along these lines is necessary to achieve a net zero carbon future (Allwood et al. 2019; Committee on Climate Change 2019).

In order to assess the potential of a geological formation as a storage reservoir it is critical to determine whether: (i) there is enough CO2 storage capacity; (ii) injected CO2 will be able to effectively move through the formation; and (iii) the mineralogy is suitable to induce precipitation reactions. In this study, we address the first two aspects of carbon storage potential of the two sites by analysing the connected porosity, pore network geometry and permeability.

Successful CCS projects typically require highly porous and permeable sedimentary formations, containing up to 30% porosity (Nelson 2004). Consequently, the majority of CCS projects are focused on sandstone and carbonate reservoirs, which commonly provide these desirable characteristics. Deep saline aquifers are of significant interest due to their global ubiquity and large capacity for which there is currently little other use (Bachu 2015; De Silva et al. 2015). CCS is also considered in igneous formations such as in the CarbFix project in Iceland (Matter et al. 2009, 2011, 2016; Gislason et al. 2010), where highly porous and reactive vesicular basalt is used as the host rock (Zahasky et al. 2018).

Hydrocarbon fields are often used either after their production lifetime has expired or during their production phase. Abandoned hydrocarbon fields are well studied and possess pre-existing infrastructure which can be repurposed for injection rather than extraction. This helps to make CO2 injection more economically viable, as is seen in the case of the Sleipner and Goldeneye fields in the North Sea (Chadwick et al. 2002; Spence et al. 2014). 73% of large-scale CCS projects in operation are based on enhanced oil recovery (EOR), such as Boundary Dam CCS, Canada, Petra Nova Carbon Capture, USA and the Jilin Oil Field CO2-EOR project, China (Page et al. 2019). EOR involves CO2 injection to displace hydrocarbons towards an extraction well to improve yield. This has the benefit of reducing the economic cost of CCS, as the initial capital investment required is offset to a degree in the long term. However, the CCS process is environmentally undermined to an extent by the use of the extracted hydrocarbons rather than assisting in net atmospheric carbon reduction.

Assessment of candidate reservoirs using digital rock physics (DRP) to measure porosity and permeability in core samples is becoming increasingly popular. DRP produces comparable results to experimental procedures, such as helium pycnometry and core flooding, without the time and cost of laboratory facilities (Van Geet et al. 2003; Jarzyna et al. 2016; Thomson et al. 2020b). Software packages such as ORS Dragonfly and PerGeos offer the capability for analysis of reconstructed 3D volumes derived from X-ray microcomputed tomographic (μCT) imaging and even the ability to run basic numerical flow simulations. Use of these software packages enables an evaluation to be made of total and connected porosity and permeability, and the production of pore network models from imaged core samples.

Previous studies examining the porosity and permeability in suites of rock samples using DRP have revealed relationships between total and connected porosity that are important to constrain when targeting formations as CCS reservoirs (Thomson et al. 2019, 2020b). Similar relationships have been found in both sandstones (Thomson et al. 2020b) and basalts (Zahasky et al. 2018) with a clear separation between partially and fully connected pore networks at 10 and up to 25%, respectively. This cutoff point has a significant influence on the ability of a geological formation to act successfully as a storage reservoir, as a reduction in connected porosity results in a decline in permeability and space for mineralization. This study aims to further constrain the connectivity cutoff point or percolation threshold, a valuable indicator of formation suitability, with a new suite of sandstone samples.

In this study we used DRP to analyse samples from the Scottish Middle Coal Measures Formation at the Glasgow UKGEOS site and the Wilmslow Sandstone Formation (WSF) in Sellafield, a lateral equivalent of the Permo-Triassic sandstones found at the Cheshire UKGEOS site. The WSF is also comparable in nature to the Utsira Formation sandstone at the Sleipner site in the North Sea (Chadwick et al. 2002). Measurements of porosity and permeability were made from digital μCT reconstructions of small core samples alongside geometrical analysis of the pore and throat structures. We aimed to assess the suitability of the geology at these sites for future CCS implementation based on the characterization of their microscale properties. This is the first detailed porosity–permeability analysis of both suites of samples to be carried out at this scale.

The dataset used in the present study is derived from scientific borehole GGC01 drilled at the UKGEOS Glasgow Geothermal Energy Research Facility Site (GGERFS) (Glasgow, UK) and a site characterization borehole, Sellafield BH13B (Cumbria, UK), located in Figure 1. We collected four core plugs (c. 85 mm in length and 25 mm in diameter) from the core material collected from borehole GGC01 at the UKGEOS site (Kearsey et al. 2019; Monaghan et al. 2019). Here, the borehole cuts through well-consolidated sandstone interbedded with siltstone and mudstone pertaining to the Scottish Middle Coal Measures Formation. We selected clean intervals, avoiding siltstone and mudstone, from a variety of depths described in Table 1.

We subsampled a further seven core plugs (c. 100 mm in length and 20 mm in diameter) from core extracted from Sellafield BH13B. The cored interval pertaining to the Wilmslow Sandstone Formation (Sherwood Sandstone Group) consists of friable fine- to medium-grained cross-bedded sandstones, red in colour (Griffiths et al. 2002; Bloomfield et al. 2006). The deepest sample, SF702, differs from the other material, and is well consolidated and pale grey in colour with a similar grain size. According to previous studies, aeolian sandstone dominates the upper portion of the formation, whereas a silicified fluvial sandstone is found in the lower portion (Evans et al. 1993; Jackson et al. 1997; Ambrose et al. 2014). Thin weakly cemented siltstone and mudstone layers are also present. These are interbedded with the aeolian sand as a result of a sandy sabkha environment, offering more consolidated sediment than the aeolian material (Griffiths et al. 2002; Mountney and Thompson 2002; Bloomfield et al. 2006). The Cheshire Energy Research Facility Site (CERFS), currently under construction as part of the UKGEOS, targets a lateral equivalent of the Wilmslow Sandstone Formation, allowing us to conduct a precursor study. We selected cores from a variety of depths as described in Table 1.

We acquired a mini-plug (10 mm in length and 5 mm in diameter) of each sample for μCT imaging and performed digital image analysis using the commercial software package PerGeos version 1.7.0 (Thermo Fisher Scientific). The methodology employed here is similar to a number of recent studies from the Royal Holloway research group (Thomson et al. 2018, 2019, 2020a, 53b). In addition, we acquired thin sections of each original larger core plug (c. 20 mm in diameter), which allowed for optical validation of different phases that are challenging to distinguish using μCT alone. An example of the material imaged in this work is given in Figure 2a.

X-ray microcomputed tomography

We acquired μCT images at the London Natural History Museum Imaging and Analysis Centre, using a Zeiss Xradia Versa 520 scanner. The scanner uses polychromatic X-rays with an energy of 80 kV to obtain images following a 5 s exposure window for each imaged interval. The image voxel size for each sample is listed in Table 1. Reconstruction of the collected data was carried out using the Zeiss Reconstructor Scout-and-Scan software, where a filtered back-projection algorithm with beam-hardening correction was used to reconstruct the 2400 projections of each sample. The software produced one stack of 16-bit tiff images per sample for segmentation in PerGeos. Figure 2 shows raw images of a single slice of each sample with each phase annotated and matching optical images in both PPL and XPL. Prior to image acquisition, the friable Sellafield samples were impregnated with epoxy. Some of these samples contained trapped air bubbles in the flooded pore space that were easily identifiable as circular features. The densities of the epoxy and trapped air are sufficiently different from the minerals that their greyscale intensities were clearly discernible during phase segmentation.

Image segmentation

The obtained images underwent binary segmentation based on their greyscale intensities in order to make further measurements. We used PerGeos to reconstruct the acquired μCT greyscale image stacks to form a 3D cylindrical volume of material. The largest possible cuboidal subvolume was extracted from each sample, removing exterior voxels not belonging to the sample and slices with significant beam-hardening artefacts. The subvolume dimensions can be found in Table 1.

We processed the raw greyscale images using a non-local means (NLM) filter (Buades et al. 2008, 2010), demonstrated in Figure 3. The filter reduces noise and improves the phase contrast for more successful segmentation. The NLM filter reduces noise by averaging similar pixels across the whole image, taking into consideration the pixels in a window around a given pixel. A threshold value is selected that determines the size of the window. Pixels which are themselves similar but also have a similar surrounding window become averaged accordingly, helping to reduce noise and improve the phase contrast in an image (Buades et al. 2008, 2010).

We used binary segmentation to identify pore space and mineral grains in each stack of filtered images. Two main methods are available to carry out segmentation: manual interactive thresholding; and automated thresholding using an algorithm such as that designed by Otsu (1979). Manual segmentation involves picking the phase threshold based on peaks on a histogram of greyscale intensities representing the X-ray attenuation. The darkest areas are representative of void space, whilst varying brighter shades of grey denote solid material where a speckled texture arises from differences in mineralogy. Since this method is prone to user bias, automatic segmentation has been used through the ‘low Otsu’ option in PerGeos which picks the darkest areas automatically, removing user bias. Despite this option, when filtering is unable to produce a significant enough phase contrast the algorithm fails to reliably separate pore space from darker solid phases.

In the present study, due to a poor phase contrast obtained during filtering, we found the automatic segmentation method unsuitable for the Glasgow samples. Consequently, data from the Glasgow samples were derived from manual segmentation. Conversely, the significantly better-defined phases present in the Sellafield samples after filtering allowed samples to be processed using automatic segmentation.

Porosity measurement

Once each sample is segmented the PerGeos ‘porosity’ tool calculates the volume fraction of each segmented phase, obtaining the total porosity from this measurement. Any phase labelled as a mineral grain was considered to be entirely impervious and non-porous, which means that any microporosity beyond the resolution of the images is not accounted for. As a result, porosity measurements are, in fact, macroporosity and therefore represent a conservative measurement.

We applied the ‘axis connectivity’ tool to the image stack, which removes any disconnected porosity in a specified orientation. The tool analyses the image slices in a stack along a specified axis, identifies unconnected pore-phase voxels – those not belonging to a connected pathway between the top and the bottom slice – and removes them. A neighbourhood value can be given to specify connectivity as six voxels with a common face, 18 voxels with at least one common edge or 26 voxels with at least one common vertex. We used a neighbourhood value of six faces, the most rigorous neighbourhood input in order to provide a connected structure which is not overestimated in size. This tool was applied sequentially in each orientation in order to obtain a fully connected pore network in all directions. Application of the ‘porosity’ tool on the fully connected sample calculates a new volume fraction of each phase, providing the connected porosity of each sample.

Pore network modelling

In order to characterize the pore space of each sample, we produced ball and stick pore network models (PNM) using a series of PerGeos algorithms following the work of Youssef et al. (2007). PNMs were produced for both total porosity and connected porosity. First, a skeletonization algorithm was applied to the image stack to generate a one-voxel-thick skeleton running through the centre of all pore space. Where lines meet, thresholds were assigned to separate individual lines, and were then classified as dead ends, pore lines or channel lines. The classification as pore or channel depends on the ratio of line length and maximum radius. If the largest possible radius of a line was greater than its length it was assigned as a pore, otherwise as a channel. This enables measurement of throat length and calculation of the coordination number of each pore. The pore radii are then measured by determining the maximum radius of a sphere that would fit in each pore section. The equivalent is done in the throat sections, providing measurements of throat radii. For further discussion of the algorithms used in PerGeos for PNM generation we refer the reader to Youssef et al. (2007) and Thomson et al. (2018, 2019). The raw data outputs from generating PNMs of each sample are available in the dataset associated with this study (Payton et al. 2020).

Absolute permeability measurement

We used the connected 3D pore space from the μCT images as the domain on which to test permeability using the ‘absolute permeability experiment simulation’ tool in the PerGeos petrophysics module. Permeability was not measured in the Glasgow samples due to the absence of connected porosity, implying the samples were impermeable. This tool uses a finite-volume solver to simulate single-phase flow of water through the connected pore space according to the Stokes equations:
u=0
(1)
P+μ2u=0
(2)
where u is velocity, P is pressure and μ is fluid viscosity equal to 1  ×  10−3 Pa s for water. The assigned boundary conditions are discussed in depth in Thomson et al. (2018). We used an error tolerance level of 10−6 for the convergence of the L2 norm of the residuals, as this is the optimal tolerance in terms of trade-off between computation time and accuracy of the result (Thomson et al. 2019). The solution obtained is a pressure and velocity field over the pore geometry from which the permeability is calculated using volume averaging. Once the permeability is evaluated from each sample containing connected pore space, we fit our permeability measurements as a function of total porosity using the Kozeny–Carman equation:
K=K0ϕn
(3)
where K is the permeability, ϕ is the error in total porosity and n is the permeability exponent. To estimate the uncertainty in permeability, we used the error in total porosity, reported in parentheses in Table 2, in equation (3) with the fit parameters K0 and n. In addition, we created a global fit to the porosity–permeability relationship, combining our data with previous results from our research group (Thomson et al. 2018, 2019, 2020b). Values of each set of fit parameters are discussed in the Results section.

Porosity

The porosities of the Glasgow samples show much lower values compared to the Sellafield samples. Figure 4 shows a range in total porosity in the Glasgow samples, from 0.04 to 1.65%, and an absence of data for connected porosity (Table 2). In contrast, a range of much larger total porosities are measured in the Sellafield samples, ranging from 9.77 to 26.4%. Each of these samples also displays significant amounts of connected porosity, ranging from 8.89 to 26.3%. The small difference in total and connected porosity for each sample, which is clearly illustrated by the near-equal-sized bars in Figure 4, indicates that the vast majority of pores in these samples form a connected network.

Despite each Sellafield sample showing significant amounts of connected porosity, we found differences in porosity distribution between samples. Figure 5 shows a 3D pore volume of both total and connected porosity extracted from the raw image stack of SF699 and SF702. Figure 5b demonstrates that sample SF702 contains a slightly higher volume of disconnected pore space compared to sample SF699. It is also worth noting that the total pore volume in sample SF699 is larger than the total pore volume in sample SF702 by nearly a factor of 3. We see that in SF699 the disconnected pores manifest as very small individual pores, showing an intragranular nature. In contrast, the disconnected porosity in SF702 is larger and more intergranular in nature, in addition to smaller, seemingly intragranular, porosity.

The relationship between total porosity and connected porosity displays a generally linear relationship, with a departure from this trend at around 10% total porosity. Figure 6 displays this relationship along with a 1:1 linear trend line. The Glasgow data points cluster near the origin due to the absence of connected porosity and small amounts of total porosity. The Sellafield data points plot closely along the 1:1 linear trend down to around 10% total porosity, where the SF702 data point deviates below the trend line on a steeper gradient. This deviation, highlighted by the dashed line, becomes clearer when considering the supporting data from similar samples (Thomson et al. 2018, 2019, 2020b). The change in gradient of the data points at around 10% total porosity is observed most clearly when plotting total porosity against the ratio of connected and total porosity in Figure 7. Similar to Figure 6, the Glasgow data points cluster near the origin due to their overall lack of porosity. In contrast, the Sellafield data create a flat plateau at total porosities >10%, with the SF702 data point showing a slightly lower ratio. The addition of supporting data shows that a significant drop-off occurs below 10% total porosity.

We show a calculated fit to the data collected in this study in Figure 7, which displays a characteristic flattening off beyond 10% total porosity, with a sharp change in gradient below this value. We use the same decaying exponential function as detailed in Thomson et al. (2020b), ϕ=ae(ϕb)/c, where ϕ is the ratio between connected and total porosity, and ϕ is total porosity, to produce the fit line. The values used for a, b and c are included in the plot. The fit using this function appears to underestimate the ratio of connected and total porosity at which the drop-off occurs for the dataset presented in this study. This may be due to the sample suite being heavily influenced by the very-low-porosity Glasgow samples. The shoulder of the fit line sits below the SF702 data point despite a more successful sigmoidal fit function being shown in Thomson et al. (2020b). The dichotomous nature of the complete set of samples analysed in this study appears to distort the fit and reduce the visibility of a sudden drop-off point, demonstrating the need for a good spread of total porosity results to effectively constrain the point at which this drop-off is observed.

Pore network models

The generation of pore network models of each sample for both total and connected porosities allowed us to measure a variety of pore-space geometry characteristics, including pore and throat radii. Histograms of pore and throat radii are shown in Figure 8, with total porosity across the top plots and connected porosity across the bottom plots. We see that in the case of total porosity the pore radii of many samples are heavily skewed towards a lower pore radius, especially in the case of the Glasgow samples. This indicates that the samples contain a very large number of very small pores relative to the number of larger pores.

Variations in average pore radius within the suite of Sellafield samples can be observed in the histograms. We see that samples SF701 and SF702 display a clear second peak at a greater pore radius, highlighted in Figure 8a1. The location of these peaks – fewer, larger pores in sample SF702 compared to sample SF699 – supports the 3D geometry of these pore spaces depicted in Figure 5. We found that 24.8 and 32.9% of pore radii in SF701 and SF702, respectively, measure <10 μm, compared to a range between 88.4 and 92.4% in the other Sellafield samples.

When examining only the pore radii that contribute to the connected pore network (Fig. 8c) the majority of the smaller pores are lost, as noted by the absence of the peaks located at values of log pore radius of c. 0.7 μm. This suggests that almost all of the small pores present in the Sellafield samples do not contribute to the connected network. This is reflected in the suite of Glasgow samples and emphasizes the importance of larger pores in enabling connectivity.

A distinctively different trend is observed in throat radii in Figure 8b and d. When considering the total porosity measurements, a reasonably normal distribution is seen, compared to the total pore radius which is generally bimodal. The Glasgow data obscure this trend to some extent, as they plot towards a lower throat radius, further explaining the low porosity measurements. In contrast, the Sellafield samples maintain a normal distribution, with SF699 tending towards a greater throat radius.

Connected throat radii in Figure 8d show a broadly similar distribution to connected pore radii. A small portion of throats are lost from the rising side of the histogram between 0 and 0.5 μm. The remainder of the data seems to remain similar compared to total porosity. This suggests that the throats with smaller radii are not part of the connected network, whilst the vast majority of throats with a log throat radius >0.5 μm contribute to network connectivity. The clear difference in trends between all pore and throat radii indicates that the two measurements are not completely dependent on one another.

The mean values of pore and throat radii for each sample, along with mean coordination number and mean throat length, were measured from the PNM of each sample in the cases of both total and connected porosity. The results for which can be seen in Table 3, using connected porosity in the case of the Sellafield samples and total porosity for the Glasgow samples. Figure 9 shows these results, allowing a visual comparison to be made between samples and total or connected porosity.

The measurements of mean pore radius show a generally positive correlation with both total and connected porosity, with the exception of SF701 and SF702 total porosity. Figure 9a shows a positive trend for the connected porosities of all samples (distinguished by hollow circles), suggesting that greater pore radii enable a greater degree of connectivity. A flatter plot is observed in total porosity of the Sellafield samples, possessing lower mean pore radii measurements due to the inclusion of smaller, disconnected pores. The mean coordination number of pores, in Figure 9c, also increases with porosity. Amongst the connected porosity results, it is apparent that a greater coordination number facilitates greater connectivity in a similar fashion to a greater pore radius. The seemingly anomalous position of the connected porosity SF701 and SF702 data points occurs again in the coordination number plot. The remaining total porosity Sellafield results show a mean coordination number of <2 that would suggest no connectivity, which is not the case. We found that samples SF701 and SF702 contain just 2 and 8% of pores with a coordination number of <2, respectively. Conversely, the other Sellafield samples have between 82 and 89% of pores with a coordination number of <2. This shows that the Sellafield samples displaying mean coordination numbers <2, and very small mean pore radii contain a large proportion of small intragranular disconnected pores that skew the results accordingly. These pores have very little impact on the measurement of porosity due to their very small volumes, as illustrated in Figure 5.

A positive correlation is observed between mean throat radius and mean throat length with both total and connected porosity, illustrated in Figure 9b and d. The data points for total and connected porosities in Sellafield plot close to one another, following a gentle positive trend. This makes it clear that wider and longer throats facilitate greater porosity measurements. However, as the total and connected porosity results plot closely, it suggests that the influence that throat geometry has on connectivity is minimal compared to pore geometry.

Permeability

Measurements of permeability performed on Sellafield samples reveal a strong positive correlation with total porosity, as shown in Figure 10. Permeability results range from 40 to 6040 mD in the Sellafield samples, with no permeability assumed in the Glasgow samples due to no connected porosity being found (Table 2). We plot supporting data from our research group again here that agrees very well with the measurements made in this study. It is clear that a greater total porosity facilitates a greater permeability.

The errors reported in permeability in Table 2 are out of proportion in the case of SF700 and SF701. This can be explained by the greater error in porosity for these samples as a result of the CT images being more difficult to perform manual segmentation on, which causes the error in calculation. This highlights that small changes in porosity can result in substantial changes in permeability.

We calculated two linear fits: one for the data presented in this work; and a global fit that considers the data published by Thomson et al. (2018, 2019, 2020b), which is the fit displayed in Figure 10. The fit calculated from the results measured in this work gives fit parameters of logk0=6.3 and n=4.7. The incorporation of supplementary data from our research group produces a better constrained global fit where logk0=5.21 and n=3.19. We show additional data from wider literature in Figure 10 which identifies a phenomenon commented on widely where a kink in the overall trend occurs at around 10% total porosity (Mavko and Nur 1997). The literature data (shown as purple diamonds) follow a linear trajectory from high porosity and permeability values down to around the 10% total porosity point. Below this the gradient of the relationship becomes significantly steeper for smaller porosities and permeabilities. The data presented in this study and that from Thomson et al. (2018, 2019, 2020b) agree well with the upper porosity linear trend, with a couple of data points beginning to follow the steeper trend below 10% total porosity. Using the results collected from the sample suites in this study, we are able to offer well-constrained Kozeny–Carman fit parameters that may be applied to future studies examining the relationship of porosity and permeability.

Porosity

Throughout our analyses we found that porosity is controlled by the geometry of both throats and pores, where development of connected porosity is most sensitive to pore geometry. Larger throat radii and greater throat length contribute to a greater porosity in the case of both total and connected porosity, shown clearly in Figure 9 as both hollow and closed symbols plot closely. Similarly, larger pore radii and greater coordination numbers of pores lead to greater total porosity. However, a discrepancy between total and connected porosity data points in the case of pores arises. This shows how the connected porosity has a greater dependence on the pore geometry than throat geometry.

We consistently found SF701 and SF702 to be outliers of sorts when examining their microstructural features and their relationships with porosity in Figures 8 and 9. These two samples display a tendency to have fewer, but larger, pores and a greater proportion of connected pores than the other Sellafield samples. This causes the total and connected porosity representations in Figures 8 and 9 to show greater agreement than that which is seen in the other samples. The difference in SF701 and SF702 is most significantly due to an observed difference in disconnected pore geometries between these samples and the remaining Sellafield material. The presence of many small intragranular pores in SF696–SF700 is illustrated in Figure 5. On average, these samples display around six times as many pores as SF701 and SF702, of which more than 88% are <10 μm in diameter compared to less than 33% in SF701 and SF702. These ubiquitous small pores are significant enough to skew our measurements of mean coordination number and pore radius to be much smaller than expected in the case of total porosity. These small pores, however, are not large enough to have a significant impact on porosity. Previous petrographical analysis of the Wilmslow Sandstone Formation suggests that minor amounts of micropolycrystalline quartz could be present (Kinniburgh et al. 2006), which may give rise to the ubiquitous small pores in SF696–SF700. Alternatively, sulfide mineralization has been studied in the Wilmslow Sandstone Formation, and Naylor et al. (1989) indicated that, during diagenesis, dissolution and removal of minerals is likely to have occurred to form overlying mineral deposits. This may provide another mechanism for the creation of many small pore spaces that are observed in the samples studied here.

Microtomography to well log tie

Our analyses of the Glasgow and Sellafield material characteristics at the microscale are well supported by bulk measurements made by well logging instruments. Figure 11 shows logging data from the GGC01 borehole including gamma ray, bulk density and calculated bulk porosity (Kearsey et al. 2019; Starcher et al. 2019). The variability in the gamma-ray signal suggests that the Formation is sand–clay mixed in composition, with elevated gamma-ray values from around 75 m indicating the possible presence of clay minerals. However, neither the μCT images nor thin sections (Fig. 2) clearly show the presence of significant quantities of clay. This may be as a result of the small sample volume that was selected in a clean and coarse sandy interval. The presence of clay minerals may give rise to microporosity, not analysed at the resolution of this study, which may result in greater porosities than presented in the Results section.

The bulk density is reasonably high, with occasional spikes, showing a background value of around 2.63 g cm−3, indicated by the shaded grey region in Figure 11. We used this value as an estimate of the matrix density due to the near-zero porosity found in all Glasgow samples when calculating bulk porosity. We took a density of 1.1 g cm−3 for the pore fluid and found that the bulk porosity showed suitably low values compared to the results from our microstructural analyses. The subsampled materials came from areas of low porosity, certainly <5% according to bulk porosity. The possibility of microporosity due to the presence of clay minerals would be reflected in a bulk measurement such as this; possibly why some bulk porosity values display around the 10% mark. The suggestion of microporosity requires further work on these samples, making use of techniques capable of resolving microporous material. This could lead to a significant increase in porosity and possibly connectivity, facilitating use in low carbon subsurface technologies.

The observations made at the pore scale are accurate for the given representative elementary volumes (REVs), which may neglect details only apparent at larger or smaller scales. The REV selected for each sample in this study was the largest flat-sided volume possible from the imaged plug. Due to the size of the samples, it is possible that observations made on these volumes could not be fully representative of the sampled material. This study aims to provide detailed observations at the macropore scale but we encourage further work on these materials at variable scales which may provide additional insight into their properties.

The well log measurements made available for Sellafield BH13B (Fig. 12) (Michie and Bowden 1994) confirm the significant amount of porosity we observed in our microstructural analyses, as well as the uniform nature of the formation. The gamma-ray track shows consistent values of around 50 API, suggestive of a clean sandstone with a higher gamma-ray interval between around 180 and 225 m that could be due to pockets of clay, followed by a cleaning-downward trend. This interval also shows increased bulk density and lower porosities. This is reflected in the difference between the locations of SF701 and SF702 on the gamma-ray track. SF701 is located before the gamma-ray increase which denotes a reduction from 24.3 to 9.77% porosity in SF701 and SF702 respectively in our microstructural analyses. SF702 was a more consolidated plug and possessed a greyer colour than all other plugs, which were friable and deep red orange in colour. There was no obvious presence of clay in the μCT images or thin sections (Fig. 2), which may contribute to better consolidation or lower porosity. A characteristic negative separation between bulk density and neutron porosity is observed in the studied interval between 63.8 and 181.39 m, indicative of sandstone. The background bulk density is lower than that measured in the GGC01 borehole (i.e. c. 2.25 g cm−3) and lower than that expected of quartz, due to the substantial porosity rather than a less dense mineralogy. Accordingly, we calculated the bulk density assuming a matrix density of 2.65 g cm−3, equal to that of quartz, and a pore fluid density of 1.1 g cm−3, equal to that of salt water. Both bulk and neutron porosities return values stable between 20 and 30%, closely agreeing with our microstructural analyses of samples SF696–SF701.

The percolation threshold

Our results show that there is a recurring change in fit line gradient at around 10% total porosity. This is true of the relationship between total porosity and connected porosity, the ratio of both porosities and permeability. This phenomenon has previously been attributed to the percolation threshold in sedimentary materials (Mavko and Nur 1997; du Plessis 1999; Gomez et al. 2010; Revil et al. 2014; Thomson et al. 2020b). Total porosity values below this 10% threshold severely inhibit connectivity, leading to partially connected networks or no connection at all, according to percolation theory (Liu and Regenauer-Lieb 2011). Above the percolation threshold a near fully connected network is maintained, allowing for greater permeability to be observed; which is what we see in our measurements.

The idea of a percolation threshold and the relationship with permeability was discussed by Mavko and Nur (1997) using a modification to the Kozeny–Carman equation to only consider porosity in a sample above the percolation threshold. They achieved a successful empirical fit to porosity–permeability data collected by Bourbie and Zinszner (1985) that reconciles the variable trend gradient above and below 10% total porosity illustrated in Figure 10. The majority of work on the percolation threshold in sedimentary material has been based on measurements from very clean sandstones, such as the Fontainebleau Sandstone. This leads to the requirement for measurements to be made on a wider variety of sedimentary materials to better constrain the percolation threshold. All of the samples plotted in Figures 6 and 7 fit with a percolation threshold of around 10% total porosity. These samples are sourced from multiple depositional environments with varying compositions. This suggests that the percolation threshold value and the relationship between total porosity and connected porosity are independent of rock composition and facies. Therefore, we show that this relationship can be used beyond the study of good-quality clean reservoir rocks and can be applied to future studies of less-conventional candidates for subsurface storage. We suggest that further research is needed to investigate this observation from a relatively small suite of samples.

The Sellafield samples, with the exception of SF702, lay beyond the percolation threshold, allowing for a very well-connected network to be observed with high permeability. The SF702 sample is observed to be on or below the percolation threshold throughout our measurements. The significantly lower permeability value measured in this sample strongly suggests that it lies just below the percolation threshold and therefore possesses a partially connected network. The absence of any connectivity and therefore permeability in the Glasgow samples clearly indicates that they are below the percolation threshold with no connected network.

Implications for subsurface carbon storage

The classification of the samples studied in this work, as either above or below the percolation threshold, has significant implications for their viabilities as candidate CCS formations. Based on the studied samples we can strongly rule out the sampled section of the Scottish Middle Coal Measures Formation in Glasgow as a suitable candidate for CCS. Other low carbon technologies requiring porous flow, such as energy storage or geothermal heating (Watson and Westaway 2020), are likely to be hindered due to the lack of porosity at the depths studied. However, further research into the presence of microporosity in these samples could result in some utility for such technologies. Conversely, the significant total porosities measured in the Sellafield samples from the Wilmslow Sandstone Formation indicate that the sampled areas would have enough capacity for CCS. Importantly, all but one of the Sellafield samples measured above the percolation threshold, providing a significant amount of connected porosity and permeability. This makes this formation a good candidate for CCS if the sampled volumes can be deemed representative of the larger area. The goal of this work was to use Sellafield as a proxy measure of the Sherwood Sandstone and Collyhurst Sandstone formations to be explored at the Cheshire UKGEOS site (Kingdon et al. 2018). Based on our observations here, indications are that the Cheshire site has great merit to be explored for CCS possibilities.

We conducted an investigation into the viability of material from the Glasgow UKGEOS site and the Wilmslow Sandstone Formation at Sellafield, as a lateral equivalent to the in-progress Cheshire UKGEOS site, as CCS reservoirs. We found that the Glasgow material possesses little porosity, up to a maximum of 1.65%, of which none is connected; indicating that this material is unsuitable for subsurface low carbon technologies at the studied depths. Conversely, the Sellafield material demonstrates significant porosity up to 26.4%, of which most is connected up to 26.3%; indicating this formation could be a strong candidate for further investigation as a CCS reservoir. We found a percolation threshold of 10% total porosity, above which near full connectivity is observed and below which connectivity rapidly decreases. This value agrees closely with existing reports in the literature. We found that permeability of the Sellafield samples ranged between 40 and 6040 mD, producing a Kozeny–Carman porosity exponent of 3.19 when coupled with supporting data from the literature, which agrees well with existing measurements of other sample suites.

The authors would like to thank Paul-Ross Thomson (RHUL) for technical PerGeos support, Frank Lehane (RHUL) for computational support, Magret Damaschke (BGS) for her comments on the manuscript and Alison Derbenwick Miller (Oracle) for facilitating computational support for the project. Samples and data are derived from the UK Geoenergy Observatories Programme and delivered by the British Geological Survey. This manuscript was published with the permission of the Executive Director of the British Geological Survey (UKRI).

RLP: conceptualization (equal), formal analysis (lead), investigation (lead), methodology (lead), project administration (supporting), validation (lead), visualization (lead), writing – original draft (lead), writing – review & editing (equal); MF: writing – review & editing (equal), resources (equal); BLC: resources (equal); DC: project administration (equal), writing – review & editing (equal), funding acquisition (equal), supervision (equal); AK: conceptualization (equal), visualization (supporting), writing – review & editing (equal), resources (equal), funding acquisition (equal), supervision (equal); SH-M: conceptualization (equal), methodology (supporting), project administration (equal), validation (supporting), writing – original draft (supporting), writing – review & editing (equal), resources (equal), funding acquisition (equal), supervision (equal).

The processed data that support the findings of this study are available from the Royal Holloway, University of London Figshare repository: https://doi.org/10.17637/rh.13401059.v1 (Payton et al. 2020). The raw image data acquired in this study is available from the National Geoscience Data Centre (NGDC) from October 2022. GGC01, Glasgow: https://doi.org/10.5285/cd213efe-a8d2-493c-9e2d-a4e6f7242eb5 and BH13B, Sellafield: https://doi.org/10.5285/dd692e0d-0df5-4eb4-9a91-e5fa0637e945.

This work was funded by the Natural Environment Research Council (GB) (NE/L002485/1), the British Geological Survey and Oracle.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)