Uncertainty Analysis of the Storage Efficiency Factor for CO₂ Saline Resource Estimation

Abstract

Carbon capture and sequestration (CCS) is a promising technology for reducing CO₂ emissions to the atmosphere. It is critical to estimate the CO₂ storage resource before deploying the CCS projects. The CO₂ storage resource is limited by both the formation pore volume available to store CO₂ and the maximum allowable pressure buildup for safe injection. In this study, we present a workflow for estimating the volume- and pressure-limited storage efficiency factor and quantifying the uncertainty in the estimates. Thirteen independent uncertain physical parameters characterizing the storage formation are considered in the Monte Carlo uncertainty analysis. The uncertain inputs contributing most to the overall uncertainty in the storage efficiency factor are identified. The estimation and uncertainty quantification workflow is demonstrated using a publicly available dataset developed for a prospective CO₂ storage site. The statistical distributions of the storage efficiency factor for the primary storage formation and the secondary storage formation located in deeper depth are derived using the proposed workflow. The effective-to-total porosity contributes most to the overall uncertainty in the estimated storage efficiency factor at the study site, followed by the maximum allowable pressure buildup, the net-to-gross thickness ratio, the irreducible water saturation, and the permeability. While the significant uncertain input variables identified are tailored to the characteristics of the study site, the statistical methodology proposed can be generalized and applied to other storage sites. The influential uncertain inputs identified from the workflow can provide guidance on future data collection needs for uncertainty reduction, improving the confidence in the CO₂ saline storage resource estimates.

1. Introduction

To mitigate the adverse impacts of climate change, technologies are being proposed to reduce the emissions of greenhouse gases such as carbon dioxide. One promising technology for reducing emissions is carbon capture and sequestration (CCS), in which CO₂ is captured from large emission sources, such as electricity power plants, and then injected into deep geologic formations. The types of geologic formations offering the potential for CO₂ storage include deep saline formations, oil reservoirs, natural gas reservoirs, and unmineable coal formations. Deep saline formations are widely spread on Earth and are estimated to have the largest storage potential globally among the geological CO₂ storage options [1]. Estimating the CO₂ storage resource is an important step in CO₂ storage site selection, and the assessment is conducted at different scales, including site-specific, regional, and basin scales.

The CO₂ storage resource ($V_{c{o}_2}$) can be expressed as the effective pore volume ($V_{po}$) times the storage efficiency factor ($E$):

$$V_{c{o}_2}=V_{po}\ast E$$

(1)

The effective pore volume ($V_{po}$) may be estimated as the product of the area, the thickness and porosity of the assessment region, or calculated from the 3D geologic model developed for the assessment region. Depending on the boundary conditions of the CO₂ storage system (i.e., open, closed, or semi-closed), different methodologies were developed for estimating the storage efficiency factor ($E$).

For closed systems in which the formation fluids cannot move out of the assessment domain, such as highly faulted and compartmentalized units (e.g., a salt cavern or a small sandstone compartment), the CO₂ storage resource is determined by the maximum allowable pressure buildup and the volume created by the compressibility of the pore spaces and the formation fluids in response to the pressure buildup. In this situation, the CO₂ storage resource can be estimated using the compressibility approach [2,3]:

$$V_{c{o}_2}=V_{po}\ast C_t\ast \nabla p$$

(2)

$$E=C_t\ast \nabla p={(C}_p+C_w)\ast \nabla p$$

(3)

where $C_t$ is the total compressibility of the formation, which is the sum of the pore compressibility ($C_p$) and the in situ water compressibility ($C_w$). $\nabla p$ is the maximum allowable pressure buildup. The storage efficiency factor ($E$) for closed systems is the product of the total compressibility and the maximum allowable pressure buildup.

Very few aquifers are truly closed, especially for saline formations, which have a relatively large extent. In most cases, saline sequestration sites are considered open systems, which are confined vertically by the sealing caprock but open laterally. For open systems, the formation water is in communication across storage and seal formations, and pressure will propagate from the injection point both outward into the storage formation and into nearly all seal formations. The two most often applied methods in the CCS field for estimating the CO₂ storage resource in open systems are the volumetric methods developed by the U.S. Department of Energy’s (DOE) National Energy Technology Laboratory (NETL) and the Carbon Sequestration Leadership Forum (CSLF) [2,4]. This study focuses on the DOE-NETL method. For a comparison of the DOE-NETL and the CSLF estimation methods, the readers are referred to IEA GHG [5]. The storage efficiency factor in the DOE-NETL method is calculated using the following equation:

$$E={{(A}_n/A_t)\ast (h}_n/h_g)\ast {(\phi }_{eff}/{\phi }_{tot})\ast E_A\ast E_l\ast E_g\ast E_d$$

(4)

where ${(A}_n/A_t)$ is the net-to-total area; ${(h}_n/h_g)$ is the net-to-gross thickness; ${(\phi }_{eff}/{\phi }_{tot})$ is the effective-to-total porosity ratio; $E_A$ is the area displacement efficiency; $E_l$ is the vertical displacement efficiency; $E_g$ is the fraction of the thickness contacted by CO₂ due to density difference between CO₂ and resident brine; and $E_d$ is the microscopic displacement efficiency. A detailed description of each multiplicative term can be found in DOE-NETL [2]. The first three geologic terms are used to define the pore volume of the entire region being assessed, which will be relatively high for site-specific assessment and could be lower for basin-scale assessment. These three parameters are referred to as “geologic terms” because they pertain to the geological characteristics of the formation and are not influenced by CO₂ injection, unlike the other parameters in Equation (4), which are specifically related to the CO₂ injection scenario. As mentioned in Goodman et al. [6], if the pore volume available to CO₂ storage has been identified in estimating the effective pore volume $V_{po}$, the first three geologic terms can be eliminated by assigning the value of one to each of them. The last four terms are used to define the pore volume near a single injector well. This volumetric method assumes that injection wells can be placed throughout the entire region under assessment, so the calculated storage efficiency factor will be the maximum value for the entire region rather than the storage efficiency factor specific to a particular injection scenario.

DOE-NETL [2] conducted six groups of Monte Carlo simulations assuming different distributions for the seven multiplicative terms and found that the resulting storage efficiency factor estimates could range from 0.22% to 10% (at the 15th and 85th percentiles) for regional-scale assessments. IEA GHG [5] estimated the storage efficiency factors at both the formation scale and the site-specific scale via numerical simulations. At the formation scale, the estimated storage efficiency factors range from 1.66% to 5.13% (at the 10th and 90th percentiles) for different lithologies, while at the site-specific scale, the estimated values range from 4% to 17% (at the 10th and 90th percentiles) for different lithologies and depositional environments. Bachu [7] reviewed several studies on the storage efficiency factor and concluded that the values of the storage efficiency factor vary widely from <1% to >10% by a factor of greater than 20, and no value can be universally applied.

As many assumptions and averaging methods are used when estimating the storage efficiency factor and large uncertainty (i.e., variability) exists in the estimates, it is critical to quantify the uncertainty associated with the storage efficiency factor. Effects of several parameters, such as the relative permeability, the irreducible water saturation, and the injection rate, on the storage efficiency factor were investigated by Gorecki and others in the IEA GHG [5] publication, but the relative importance of these input parameters have not been compared. Sarkarfarshi et al. [8] performed a sensitivity analysis to compare the contribution of five physical parameters and four constitutive parameters on predicted uncertainty in the CO₂ plume evolution using a finite-element model. Afanasyev et al. [9] evaluated the CO₂ storage potential for reservoir sectors of West Siberia by conducting numerical reservoir simulations in scenarios with single and ten injection wells. The estimated maximum storage efficiency factor ranges between 0.03 and 0.16 for different depositional environments. Calvo et al. [10] investigated the uncertainties in CO₂ density caused by the uncertainties in subsurface temperature and pressure.

Detailed numerical multi-phase flow simulations have the ability to include a detailed geologic description and fluid flow in the subsurface, but they are resource-consuming and may not be suitable for Monte Carlo uncertainty quantification with multiple uncertain inputs, which requires tens of thousands of simulation runs. Nordbotten et al. [11] developed an analytical solution for the evolution of the CO₂ plume during CO₂ injection. The analytical solution allows for quick estimation of the maximum plume extent and the shape of the CO₂-brine interface, which could be used in Monte Carlo simulations for uncertainty quantification. Building upon the work of Nordbotten et al., Okwen et al. formulated a straightforward analytical equation for estimating the CO₂ storage efficiency [12].

The CO₂ storage resource is limited by both the pore volume available for storing CO₂ and the maximum allowable pressure buildup for safe injection. The commonly applied volumetric estimation methods for open systems do not take the pressure limit into account. This study aims to develop a framework for estimating the storage efficiency factor in the DOE-NETL volumetric method, with the pressure limit taken into consideration, as well as quantifying the uncertainty in the storage efficiency factor estimates considering a substantial number of uncertain parameters. We demonstrate the estimation and quantification framework using a publicly available dataset developed for the Smeaheia site located on the Hordaland Platform in the northern North Sea, offshore Norway [13]. The analytical solution to CO₂ plume evolution during injection developed by Nordbotten et al. [11] is employed for the quick estimation of the storage efficiency factor, which allows for the associated distribution being estimated by 10,000 Monte Carlo simulations with 13 independent uncertain inputs. Important input parameters contributing most to the overall uncertainty in the estimated storage efficiency factor are identified. While the significant uncertain input variables identified are tailored to the characteristics of the study site, the statistical methodology proposed can be generalized and applied to other storage sites.

2. Methods

2.1. Estimation of the Storage Efficiency Factor

The storage formations contributing to the storage volume of the assessment region need to be identified, and then the storage efficiency factor is estimated for each storage formation considered. As shown in Equation (4), the storage efficiency factor in the DOE-NETL volumetric method is a product of seven individual terms. The first three terms are geologic terms for refining the fraction of pore space amenable to CO₂ storage when the effective pore volume ($V_{po}$) is estimated by multiplying the area, thickness, and total porosity of the assessment region. The net-to-total area term $A_n/A_t$ is the fraction of the total assessment area where a suitable storage reservoir is present. If there are faults extending vertically upward and intersecting the sealing caprock formations, the permeability and stability of these faults need to be assessed carefully to ensure safe CO₂ storage in subsurface reservoirs. The areas close to the faults that are unfavorably oriented for reactivation and are likely to create leakage pathways should be excluded from the net area ($A_n$) of the storage reservoir. The net-to-gross thickness term $h_n/h_g$ is the fraction of thickness that has properties (e.g., porosity, permeability, and thickness) suitable for CO₂ storage, which can be estimated based on the net-to-gross ratio of the reservoir interval from wells at the storage site; the effective-to-total porosity ${\phi }_{eff}/{\phi }_{tot}$ is the ratio of the effective porosity (i.e., interconnected pore space) to the total porosity, which can be estimated based on petrophysical analysis of well logs and core analysis.

The area displacement efficiency term $E_A$ and the vertical displacement efficiency $E_l$ represent the fraction of the area and the thickness that can be contacted by the CO₂ plume from a single well due to permeability heterogeneity. To ensure safe injection at the storage site, the pressure buildup should not exceed the maximum allowable pressure buildup above which new fractures could be created or existing fractures could be activated in the caprock. The pressure buildup caused by CO₂ injection is determined by the injection rate, the injection duration, the reservoir permeability, the reservoir thickness, and the mobility (i.e., the ratio of relative permeability to fluid viscosity) of the CO₂ phase and the brine phase [11]. Relative permeability is defined for multiphase flow as the ratio of the effective permeability of a particular fluid at a given saturation to the intrinsic permeability of the rock [14]. The relative permeability to the brine phase is uniformly 1 in the analytical model due to the assumption that the brine zone is fully saturated with brine, while the maximum relative permeability to the CO₂ phase (i.e., the end point relative permeability of CO₂ at the maximum CO₂ saturation) is less than 1 as the result of residual (irreducible) brine saturation behind the CO₂–brine interface. To ensure safe injection and to meet the target injection rate at the same time, the permeability of the CO₂ storage reservoir should be larger than a critical value (${perm}_{crit}$), which can be estimated using Equations (5)–(7) based on Nordbotten et al. [11].

$${perm}_{crit}=\frac{-Q}{2\pi B\nabla p}\ast \left\{\frac{1}{{\lambda }_c}ln\left(\frac{r_{well}}{a_2}\right)+\frac{B}{b_1{\lambda }_c+\left(B-b_1\right){\lambda }_w}ln\left(\frac{a_2}{a_1}\right)+\frac{1}{{\lambda }_w}ln\left(\frac{a_1}{R}\right)\right\}$$

(5)

$$a_1=\sqrt{\frac{{\lambda }_{c}Qt}{{\phi \pi (b}_1{\lambda }_c+\left(B-b_1\right){\lambda }_w)}}$$

(6)

$$a_2=\sqrt{\frac{{\lambda }_{w}Qt}{{\phi \pi (b}_1{\lambda }_c+\left(B-b_1\right){\lambda }_w)}}$$

(7)

where $B$ is the thickness of the storage reservoir; $b_1$ is the thickness of the top layer of the CO₂ plume; $\nabla p$ is the maximum allowable pressure buildup; $r_{well}$ is the radius of the injection well, assumed to be 0.09144 in this study; $Q$ is the volumetric injection rate, which is assumed to be constant in time; $t$ is the injection duration; ${\lambda }_c$ and ${\lambda }_w$ are the mobility of the CO₂ phase and the brine phase, respectively; $\phi$ is the porosity; and $R$ is the expanding outer boundary of the domain, calculated as a function of the injection duration ($t$), the intrinsic permeability ($perm$) and porosity ($\phi$) of the storage reservoir, viscosity ($\mu$) of the fluid, rock compressibility ($C_r$), and the in situ water compressibility ($C_w$) (Equation (8)) based on Nordbotten et al. [15].

$$R=\sqrt{\frac{4U_c\times perm\times t}{\mu \times (C_r+\phi C_w)}}$$

(8)

where $U_c$ is a cutoff value, set at 0.56 in Nordbotten et al. [15]. The water compressibility can be estimated as a function of the in situ temperature and pressure condition using the CoolProp library [16,17].

Once the minimum permeability value (${perm}_{crit}$) required for safe injection at the target injection rate is estimated, $E_A$ and $E_l$ can be estimated by calculating the fraction of the total area and the total thickness where the permeability is larger than ${perm}_{crit}$, provided that the horizontal and the vertical heterogeneities within the reservoir are characterized. In this study, by comparing the permeability values in the top layer of the subsurface 3D heterogeneous geologic model developed for the assessment region with the estimated minimum permeability value (${perm}_{crit}$) required for safe injection, $E_A$ is computed as the fraction of the area with permeability larger than ${perm}_{crit}$ to the total area of the formation top. $E_l$ is estimated as the average ratio of the thickness with permeability larger than ${perm}_{crit}$ to the net thickness considering all the surface locations in the heterogenous geologic model.

The efficiency term $E_g$ further constrains the CO₂ plume volume in the vertical cross section, considering CO₂ rises within the geologic unit due to gravity (buoyance), viscous, and pressure forces. In this work, the analytical solution for the CO₂ injection via a single well into a horizontal, homogenous, and isotropic aquifer developed by Nordbotten et al. [11] is utilized to estimate $E_g$. Based on the analytical solution, the ratio of the CO₂ plume thickness to the total thickness of the formation ($b^{\prime }=\frac{b}{B}$), ranging from 0 to 1 at a particular location, is a function of the density difference between CO₂ and in situ brine ($\mathsf{\Delta }\rho$), the radial distance ($r$) from the injection well, injection rate ($Q$), injection duration ($t$), mobility of the CO₂ phase (${\lambda }_c$) and brine phase (${\lambda }_w$), porosity ($\phi$), permeability ($k$), and reservoir thickness ($B$) (Equations (9)–(13)):

$$-\frac{\lambda -1}{{r^{\prime }\left(\left(\lambda -1\right)b^{\prime }+1\right)}^2}+2\mathsf{\Gamma }{r}^{\prime }{b}^{\prime }+2\mathsf{\Lambda }{r}^{\prime }=0$$

(9)

$$\mathsf{\Lambda }({\lambda -1)}^2-\mathsf{\Gamma }\left(\lambda -1\right)+\mathsf{\Gamma }\lambda \mathrm{l}\mathrm{n}(\frac{\mathsf{\Gamma }+\mathsf{\Lambda }}{\mathsf{\Lambda }\lambda })=\frac{2\lambda {\left[\mathsf{\Lambda }\left(\lambda -1\right)-\mathsf{\Gamma }\right]}^2}{\lambda -1}$$

(10)

where $\mathsf{\Gamma }$ is the dimensionless gravity factor defined as

$$\mathsf{\Gamma }=\frac{2\pi \mathsf{\Delta }\rho g{\lambda }_{w}k{B}^2}{Q}$$

(11)

$\lambda$ is the ratio of mobilities of the two fluid phases:

$$\lambda =\frac{{\lambda }_c}{{\lambda }_w}$$

(12)

and $r^{\prime }$ is a dimensionless variable defined as:

$$r^{\prime }=r\sqrt{\frac{\pi B\phi }{Qt}}$$

(13)

The $\mathsf{\Lambda }$ in Equation (10) is the Lagrangian multiplier, which can be solved by numerical methods for given fluid and rock properties, as well as injection parameters. The ratio of the CO₂ plume thickness to the total thickness of the formation ($b^{\prime }=\frac{b}{B}$) is solved using Equation (9) once the value of $\mathsf{\Lambda }$ is obtained.

The maximum lateral extent of the CO₂ plume at the end of the CO₂ injection period is estimated by Equation (14):

$$r_{max}=\sqrt{\frac{{\lambda }_{c}Qt}{{\lambda }_w\pi B\phi }}$$

(14)

where $r_{max}$ is the maximum lateral extent of the CO₂ plume at time $t$; $Q$ is the volumetric injection rate; $\phi$ is the porosity; ${\lambda }_c$ and ${\lambda }_w$ are the mobility of the CO₂ phase and the brine phase, respectively; and $B$ is the reservoir thickness. We then estimate the efficiency term $E_g$ by taking the average of the calculated $b^{\prime }$ within the maximum lateral extent of the CO₂ plume where $b^{\prime }>0$.

The density and the viscosity of CO₂ are a function of the temperature and the pressure of the storage reservoir, which can be estimated using thermophysical property libraries, such as the NIST website [18] or the CoolProp package [14,15]. In this study, the CoolProp package is used in the estimation as it is open-source and can be included in the code for Monte Carlo uncertainty analysis. A comparison of the estimated CO₂ properties using the CoolProp package and those using the NIST website shows that the estimations obtained from the two libraries are consistent with each other (Figure S1 in the Supplementary Materials). The density and the viscosity of brine, which is a function of the in situ temperature, pressure, and salinity, can be estimated using the empirical relationship developed by Batzle and Wang [19] and Mao and Duan [20,21].

The last efficiency term $E_d$ represents the portion of the pore volume that can be replaced by CO₂. This efficiency term is directly related to the irreducible water saturation ($S_{wir}$) and is computed as $1-S_{wir}$ in this study.

2.2. Monte Carlo Uncertainty Analysis

The uncertainty (i.e., variability) of a variable can be represented by a statistical distribution. The Monte Carlo uncertainty quantification method estimates the distribution of an output variable (i.e., the storage efficiency factor in this study) by drawing thousands of random samples from the distribution of each input parameter to calculate all possible outcomes of the output variable.

To derive the distribution of the storage efficiency factor for uncertainty quantification, 10,000 estimates of the storage efficiency factor for each storage formation at the assessment region, considering the variability of the input parameters, are generated using the analytical approach described above. The parameters that affect the storage efficiency factor estimates can be explicitly identified from the equations in the analytical model. Some parameters in the analytical model are correlated with each other. For example, the density and the viscosity of CO₂ are correlated, and both depend on the temperature and the pressure of the storage reservoir. The in situ temperature and the pressure condition at the storage reservoir, both of which depend on the reservoir depth, can be estimated as a function of the depth of the storage reservoir, the geothermal gradient, and the pressure gradient of the storage site. As independence in uncertain inputs is an assumption in the Monte Carlo uncertainty analysis method, we select 13 independent or weakly correlated variables as uncertain inputs in the Monte Carlo uncertainty analysis, which are all intrinsic physical parameters of the storage formation. The parameters related to the injection operation are treated as fixed parameters. The 13 intrinsic physical parameters of the storage formation include the depth, thickness, porosity, permeability, geothermal gradient, pressure gradient, salinity of the resident brine, irreducible water saturation, maximum relative permeability to CO₂, compressibility of the formation rock, maximum allowable pressure buildup, net-to-gross ratio, and effective-to-total porosity. In this paper, the maximum relative permeability to CO₂ refers to the end point relative permeability of CO₂ when the irreducible water saturation is reached (i.e., at the maximum CO₂ saturation) during drainage.

For each uncertain parameter, the probability density function is estimated by fitting statistical distributions to measurements conducted for the assessment region or reported values for other geologic formations in the literature if the site-specific information is unavailable. In this study, truncated Normal distributions are assumed for all the uncertain input parameters except for the permeability and the compressibility of the formation rock, which are assumed to follow truncated lognormal distributions. For some uncertain parameters considered in this study, such as the salinity of the formation water and the depth of the formation, only range values or few measurements are available. The reported most likely values for these parameters or the average of the reported upper and the lower limits are considered as the mean values. The reported ranges for these parameters are used to derive the standard deviations of these random variables, assuming the estimated ranges are approximately four times the standard deviations. The upper limit and lower limit of the measurements for uncertain input parameters are used in the truncation when generating random variables following the fitted Normal or Log-normal distributions. A total of 10,000 random samples are drawn from the specified distribution for each uncertain input variable, resulting in 10,000 sample vectors $x=\left\{x_1, x_2, \dots , x_n\right\}$ (n = 13 in this study) assuming the n variables are independent. For each sample vector (i.e., each scenario of the CO₂ storage formation), $x_j=\left\{x_{1j}, x_{2j}, \dots , x_{nj}\right\}$, j = 1, 2, 3, …, 10,000, the storage efficiency factor $E_j$ is estimated using the method described in Section 2.1. Based on the 10,000 estimations, statistical measures such as the mean, the standard deviation, and percentiles of the storage efficiency factor are computed. The distribution of the storage efficiency factor at the storage site is estimated and can be presented graphically in histograms and plots of the empirical cumulative density functions. Figure 1 is a flowchart illustrating the workflow for rapid estimation and quantification of the uncertainty in the storage efficiency factor. Besides estimating the overall uncertainty in the storage efficiency factor, the input variables that contribute most to the overall output uncertainty are also identified from the 10,000 Monte Carlo simulations by calculating the correlation coefficient between the storage efficiency factor and each uncertain input. The absolute value of the correlation coefficient provides a measure of the contribution of each input to the overall uncertainty in the storage efficiency factor. The input parameters are ranked by the absolute value of the corresponding correlation coefficient, with the higher absolute value indicating a larger contribution to the overall uncertainty in the estimated storage efficiency factor from that input.

3. Data

We demonstrate the proposed estimation and uncertainty quantification workflow for the storage efficiency factor using a real-world example based on the prospective Smeaheia offshore CO₂ storage site located in the Norwegian North Sea. The Smeaheia dataset, developed by Equinor ASA and Gassnova SF, includes seismic data, well reports, pressure and temperature data, geomechanical and stress data, 3D geologic models, and reservoir simulations [13]. This dataset was made publicly available by SINTEF. The potential storage formations considered at the Smeaheia storage site consist predominantly of sandstones and include the Upper Jurassic Sognefjord, Fensfjord, and Krossfjord formations. The Sognefjord formation is identified as the primary storage unit. The Draupne formation, consisting primarily of impermeable claystones is considered the primary seal for CO₂ storage. The main faults at the Smeaheia site are the Vette fault to the north and west and the Øygarden Fault Complex (ØFC) to the east. As the eastern Øygarden fault is identified as a potential leakage pathway for CO₂ storage, the structure closed in by the Øygarden fault could be excluded from the net suitable storage volumes. Petrophysical interpretations based on well log analysis of seven wells at the study site suggest the effective porosity of the Sognefjord formation ranges from 0.21 to 0.31 with a mean of 0.26. The permeability of the Sognefjord formation ranges from 34.27 mD to 2989 mD with a mean of 1671 mD. Figure 2 shows the spatial distribution of permeability and porosity in the 3D heterogeneous geologic model included in the Smeaheia dataset. Example 2D slices of the 3D permeability field realization are provided in Figure 3. The Sognefjord formation varies in depth from approximately 880 m to 1300 m below the mean sea level, with a thickness range from 68 m to 170 m [13,22]. The depth and thickness of the Fensfjord formation range from 990 m to 1400 m and from 103 m to 230 m, respectively [13,22]. The proposed injection rate is 3.2 million tons per year, and the injection duration is 50 years. The pressure and the temperature gradient are approximately 10,500 Pa/m and 0.0346 °C/m, respectively, based on the temperature and pressure measurements conducted at the storage site [13]. The estimated maximum allowable pressure buildup at the Smeaheia site varies widely from 17.7 bar (based on the minimum horizontal stress assuming pre-existing fractures in the caprock formation) to 75 bar (based on the tensile strength assuming intact caprock formation), and the estimated most likely value is 33 bar [13]. The rock compressibility is uncertain, with the estimated values at the study site ranging from 1.6 × 10⁻⁶ bar⁻¹ to 1.6 × 10⁻⁴ bar⁻¹. Based on the reported salinity values in the northern North Sea [23], the salinity of the formation water is about 60,000 ppm.

4. Uncertain Inputs

We calculate the storage efficiency factor for each storage formation considered at the prospective CO₂ storage site. Due to limited information available about the Krossfjord formation, the storage efficiency factor for the Krossfjord formation is not calculated in this study. Assuming the net area ($A_n$) of the storage formations has been characterized at the study site and the total pore volume $V_{po}$ is computed from the 3D geologic model that represents the net suitable storage volumes with faults taken into consideration, the net-to-total area term $A_n/A_t$ is assigned the fixed value of 1 in the evaluation.

The relative permeability experimental measurements for sandstone compiled in Burnside and Naylor [14] and those conducted by Crandall et al. [24] are used to characterize the uncertainties in the irreducible water saturation ($S_{wir}$) and the maximum relative permeability to CO₂ ($k_{r,co2}$). A scatter plot of the physical experimental measurements of the maximum relative permeability of CO₂ and the irreducible water saturation is provided in Figure S2 of the Supplementary Materials. The measurements of $S_{wir}$ and $k_{r,co2}$ vary widely for sandstone. $S_{wir}$ ranges between 0.2 and 0.7, with a mean of 0.49 and a standard deviation of 0.12. The maximum $k_{r,co2}$ values are within a range of 0.02–0.96, with a mean of 0.3 and a standard deviation of 0.24.

Well logs included in the Smeaheia dataset for wells 32/2-1 and 32/4-1 are used to estimate the distributions of the effective-to-total porosity. Based on Fawad et al. [22], the depth interval of the Sognefjord formation in wells 32/2-1 and 32/4-1 range from 880 m to 990 m and from 1210 m to 1280 m, respectively. The depth interval of the Fensfjord formation in wells 32/2-1 and 32/4-1 are [990 m, 1090 m] and [1340 m, 1570 m], respectively. The measurements of the effective porosity and the total porosity from the two wells at corresponding depth intervals for the Sognefjord formation and Fensfjord formation are extracted from the well logs.

The net-to-gross thickness ratio, thickness, depth, permeability, and porosity data for the Sognefjord formation and Fensfjord formation at seven well locations are obtained from the subsurface evaluation report included in the Smeaheia dataset and Fawad et al. [22]. The net-to-gross thickness ratio was obtained by applying cutoff values for the volume of shale, the effective porosity, and the permeability of 0.3, 0.1, and 20 mD, respectively. The Fensfjord formation is deeper and, on average, thicker than the Sognefjord formation, while the porosity and the permeability are lower than that of the Sognefjord formation. The means of the measured effective-to-total porosity for the Sognefjord formation and Fensfjord formation based on two well logs at the study site are 0.746 and 0.712, respectively. As the formation-specific data are not available for other uncertain parameters considered, the distributions of other uncertain parameters are assumed to be the same for the Sognefjord formation and Fensfjord formation. Typical geothermal gradients and pressure gradients vary in the range of 0.02 to 0.06 °C/m and 10 to 12 kPa/m, respectively [7,25]. Bjørlykke and Gran [23] collected the salinity measurements from well logs in the North Sea and found that the salinity values of formation waters in the northern North Sea vary from 20,000 to 100,000 ppm.

Table 1. Uncertain and fixed parameters assumed in Monte Carlo uncertainty analysis for storage efficiency factors of the Sognefjord formation and the Fensfjord formation. For the statistics of the variables, if the values for the Sognefjord formation and the Fensfjord formation are the same, only one value is listed in the table. The mean and standard deviation values represent parameters of the derived distributions, not statistics of the actual measurements collected at the study site. The upper and lower limits are obtained from actual measurements at the study site or reported values in the literature.

Uncertain Parameter	Distribution Type	Mean a (Sognefjord/Fensfjord)	Standard Deviation a (Sognefjord/ Fensfjord)	Range (Sognefjord/Fensfjord)	Unit	Data/ References
Net-to-gross ratio	Truncated normal	0.748/0.742	0.222/0.178	0.36 to 0.97/0.52 to 0.97	-	[13,22]
Effective-to-total Porosity	Truncated normal	0.746/0.712	0.230/0.277	0.125 to 1.0/0.019 to 1.0	-	[13]
Depth	Truncated normal	1090/1195	105/102.5	880 to 1300/990 to 1400	m	[13,22]
Thickness	Truncated normal	119.0/166.5	25.5/31.75	68 to 170/103 to 230	m	[13,22]
Geothermal gradient	Truncated normal	0.0346	0.01	0.02 to 0.06	°C/m	[13,25]
Pressure gradient	Truncated normal	10,500	500	10,000 to 12,000	Pa/m	[7,13]
Salinity of brine	Truncated normal	0.06	0.02	0.02 to 0.1	weight fraction	[23]
Irreducible water saturation	Truncated normal	0.496	0.120	0.2 to 0.7	-	[14,24]
Maximum relative permeability to CO2	Truncated normal	0.270	0.238	0.018 to 0.96	-	[14,24]
Porosity	Truncated normal	0.259/0.229	0.036/0.030	0.21 to 0.31/0.19 to 0.28	-	[13,22]
Permeability	Truncated Lognormal	1671.3/755.5	6768.6/4520.5	34.27 to 2989/10.58 to 1729	mD	[13,22]
Maximum allowable pressure buildup	Truncated normal	33	14.325	17.7 to 75	bar	[13]
Rock compressibility	Truncated Log10 normal	−4.796 b	0.5 b	−5.796 to −3.796	/bar	[13]
Fixed Parameter		Value		Unit	Data/Reference
Injection rate		3.2		million tonnes/year	[13]
Injection duration		50		year	[13]
Estimated storage efficiency factor for the Sognefjord formation		Mean: 0.0643 Standard deviation: 0.0359 Median (50th percentile): 0.0574 15th percentile: 0.030 85th percentile: 0.100		Estimated storage efficiency factor for the Fensfjord formation	Mean: 0.0817 Standard deviation: 0.0431 Median (50th percentile): 0.0754 15th percentile: 0.039 85th percentile: 0.124

^a The mean and standard deviation values are parameters of the distributions before truncation. ^b The mean and standard deviation values are parameters of the base 10 logarithm of the variable before truncation.

summarizes the key independent variables affecting the estimated storage efficiency factor, including 13 physical parameters of the storage formation and two parameters related to the injection operation (i.e., the volumetric injection rate and the injection duration). The derived distributions for uncertain input parameters and the assigned values for fixed parameters in the Monte Carlo uncertainty analysis for estimating the storage efficiency factor are provided. For all the uncertain input parameters, except for the permeability and the rock compressibility, the truncated normal distributions are fitted to the experimental data and reported typical values. For the compressibility of formation rock, the truncated lognormal distribution in base 10 is assumed. The base 10 logarithm of the reported upper and lower limits of the rock compressibility measurements are used to estimate the mean and the standard deviation of the normal distribution for the base 10 logarithm of the rock compressibility. For the permeability of storage formation, the truncated lognormal distribution is fitted directly to the permeability data obtained from the subsurface evaluation report [13] and Fawad et al. [22].

5. Results

Once the distributions for uncertain input parameters are specified, random samples of the uncertain parameters can be drawn from the specified statistical distributions. Each sample is drawn independently, and the samples for different uncertain input parameters are drawn from different distributions without introducing any correlation or dependence between them explicitly; therefore, the resulting uncertain input variables can be considered independent. We further calculate Spearman’s rank-order correlation coefficient between the generated uncertain input variables, and the correlation coefficient is around $\pm$0.01. Ten thousand runs of Monte Carlo simulations of the storage efficiency factor are performed using the proposed estimation method and the specified statistical distributions for uncertain input parameters. The resulting distributions of the estimated storage efficiency factors for the Sognefjord formation and Fensfjord formation are illustrated in Figure 4 using the histograms and the empirical cumulative density function (CDF) plots. The mean, standard deviation, and median value of the estimated storage efficiency factor $E$ for the Sognefjord formation are 0.064, 0.036, and 0.057, respectively. For the Fensfjord formation, the mean, standard deviation, and median value of the $E$ estimation are 0.082, 0.043, and 0.075, respectively. The estimated storage efficiency factor varies from 0.030 to 0.100 and from 0.039 to 0.124 (at the 15th and 85th percentiles) for the Sognefjord formation and Fensfjord formation, respectively. The mean of the storage efficiency factor for the Fensfjord formation is larger than that for the Sognefjord formation, and the associated uncertainty in the $E$ estimation indicated by the standard deviation value is also higher compared with that for the Sognefjord formation. Spearman’s rank-order correlation coefficient [26] is calculated between the storage efficiency factor and each uncertain input in the Monte Carlo simulations, measuring the contribution of each input to the overall uncertainty in the output. Scatter plots of the calculated $E$ in the 10,000 Monte Carlo simulation runs against the generated net-to-gross ratio, effective-to-total porosity, irreducible water saturation, and depth values for the Sognefjord formation in the 10,000 independent scenarios are included in Figure S3 of the Supplementary Materials. The uncertain inputs are ranked by the absolute value of the rank-order correlation coefficients. The rank of uncertain input parameters and the computed correlation coefficient between the storage efficiency factor and each input are illustrated in Figure 5. The net-to-gross thickness ratio (net_gross), effective-to-total porosity ratio (phe_phit), maximum allowable pressure buildup (diff_P_max), formation thickness, and rock compressibility (Cr) exhibit a positive correlation with the storage efficiency factor ($E$). Conversely, the irreducible water saturation (Swir), geothermal gradient (G_ther), pressure gradient (G_pres), and salinity of brine show a negative correlation with $E$, which are evident from the analytical solution utilized in this study (Equations (4)–(14)). Permeability demonstrates a negative correlation with the storage efficiency factor, primarily due to its influence on CO₂ migration dynamics within the reservoir rock. Higher permeability facilitates more rapid CO₂ migration, leading to a larger spread toward the top of the storage formation and reduced thickness occupied by the CO₂ plume in the vertical cross-section (reflected in the lower $E_g$ efficiency term in Equation (4)). For the Sognefjord formation, the relative permeability of CO₂ (Krc) exhibits a positive correlation with $E$, while for the Fensfjord formation, it shows a negative correlation with $E$. The impact of Krc on $E$ is multifaceted. On one hand, higher relative permeability enhances the mobility of the CO₂ phase, potentially reducing the critical permeability threshold required for safe injection and the target injection rate (as seen in Equation (5)) and allowing for more effective CO₂ retention within pore spaces. However, increased CO₂ mobility may also result in a thinner CO₂ plume in the vertical cross-section due to buoyancy effects, thereby reducing the $E_g$ efficiency term in Equation (4). Ultimately, the overall influence of Krc on $E$ depends on the dominant mechanism within the specific storage reservoir. The formation depth is positively correlated with the storage efficiency factor because, at greater depths, where the density and viscosity of the CO₂ phase are higher, CO₂ mobility tends to decrease, promoting better contact with the rock matrix and potentially leading to a higher $E_g$ efficiency.

For both the Sognefjord formation and Fensfjord formation, the effective-to-total porosity contributes most to the uncertainty in the storage efficiency factor. The top five uncertain parameters contributing most to the uncertainty in the $E$ estimations are the effective-to-total porosity, maximum allowable pressure buildup, net-to-gross thickness ratio, irreducible water saturation, and permeability, although the exact ranking of these parameters is different for the Sognefjord formation, compared with those for the Fensfjord formation. The effective-to-total porosity represents the interconnected pore space contributing to fluid flow to the total porosity, which includes all pore space within the rock regardless of connectivity. The net-to-gross thickness ratio compares the net thickness of the reservoir interval, which contains porous and permeable rock suitable for fluid storage, to the gross thickness, which includes all rock layers within the interval regardless of porosity or permeability. Higher ratios indicate greater connectivity and continuity of porous and permeable rock, suggesting a higher potential for fluid flow. Conversely, lower ratios suggest more heterogeneous or compartmentalized reservoirs with reduced fluid flow potential. As both the effective-to-total porosity and net-to-gross thickness ratio measure the degree of heterogeneity in the permeability and porosity, the uncertainty in the permeability, porosity, allowable pressure buildup, and irreducible water saturation contribute most to the uncertainty in the storage efficiency factor. Reducing the uncertainty in the permeability, porosity, allowable pressure buildup, and irreducible water saturation at the storage site should allow the storage efficiency factor to be estimated with higher confidence.

6. Discussion

A comprehensive uncertainty analysis using 13 independent uncertain input parameters was performed to estimate the uncertainty in the storage efficiency factor and to evaluate the contribution of each uncertain input to the overall uncertainty in the storage efficiency factor. The estimated storage efficiency factor values using the proposed method are within the ranges of previous estimations [5,7,9]. The method introduced in this study enables quick derivation of the statistical distribution of the storage efficiency factor at a storage site, which in turn facilitates statistical quantification of the CO₂ storage capacity for saline aquifers. A substantial number of uncertain parameters were included in the Monte Carlo uncertainty analysis to investigate the contribution of different parameters. The ranks of the uncertain inputs in terms of their individual contribution to the overall uncertainty in the storage efficiency factor depend on the degree of variability associated with the uncertain inputs. A high rank for an uncertain input variable signifies significant variability within that variable. The proposed workflow allows for the distribution of the storage efficiency factor to be estimated, with the variability of all the uncertain input parameters taken into consideration simultaneously. Important uncertain parameters identified from the uncertainty analysis can help decision makers prioritize data collection during aquifer and site characterization and provide guidance on further data collection needs for uncertainty reduction in the storage resource estimates.

As the heterogeneity in the permeability and porosity is a key factor affecting the uncertainty in the storage efficiency factor, it is important to characterize the variability of the porosity and permeability at the assessment region and construct heterogeneous geologic models rather than homogenous models for the storage formations. The irreducible water saturation, $S_{wir}$, is another parameter that contributes a lot to the uncertainty in the storage efficiency factor. There are large variations in the experimental measurements on $S_{wir}$. As mentioned in Burnside and Naylor [14], issues with the experimental set-up and the running conditions may produce experimental biases on the $S_{wir}$ measurements and systematic differences between laboratories. The experimental flow rate has a significant effect on the $S_{wir}$ measurements, as evidenced by measurements conducted by Crandall et al. [24]. To reduce the degree of variability in $S_{wir}$ and the overall uncertainty in the storage efficiency factor, the formation-specific distribution of $S_{wir}$ at in situ conditions needs to be derived for the assessment region, representing heterogeneity in the reservoir formation and suppressing experimental biases. It is recommended that core samples from multiple wells within the same formation are collected and analyzed under identical experimental conditions matching the in situ injection conditions when conducting the $S_{wir}$ measurements at the storage site. Multiple core samples from single rock units at the study site and repeated physical experiments on single samples would allow site-specific statistical distributions for the uncertain parameters to be derived, furthering our understanding of the uncertainty in the estimated storage efficiency factor.

The publicly available dataset utilized in this study is developed specifically for the Smeaheia storage site, making the estimated storage efficiency factor values site-specific. However, for parameters lacking values at the study site, we resort to values from other sites as approximations. For some of the uncertain parameters, such as the formation thickness, the maximum allowable pressure buildup, and the salinity of brine, only a few data points are available for estimating the means and the standard deviations of the truncated normal distributions. Given the limited number of actual measurements available at the study site, the standard deviations of the derived distributions can be quite large, resulting in a relatively rough estimation of the uncertainty in these input parameters. When more measurements on these physical properties become available, the distributions of these properties can be updated by fitting to the new experimental data.

We note that the analytical solution to CO₂ plume evolution used in this study assumes the aquifers are horizontal, homogenous, and isotropic. This assumption of highly simplified geology can lead to errors in the estimated maximum extent of the CO₂ plume, resulting in biases in the storage efficiency estimation. At the Smeaheia site, which is a heterogenous and tilted formation (Figure 2), there should be a preferential motion of the CO₂ upward in the tilted strata.

The vertical-to-horizontal permeability ($k_v/k_h$) anisotropy was not accounted for in this study as vertical equilibrium is assumed in the derivation of the analytical solution [11,27]. As mentioned in Nordbotten and Celia [27], the vertical equilibrium assumption becomes more valid as the CO₂ plume expands and the horizontal length scale increases. To account for the effect of $k_v/k_h$ anisotropy on the storage efficiency factor, numerical simulators could be used instead of the analytical solution. However, using the analytical solution allows for key parameters affecting the storage efficiency factor to be identified, as well as the statistical distributions of the storage efficiency factor being produced via tens of thousands of simulations using the Monte Carlo method.

An offshore CO₂ storage site was used in this study to demonstrate the method for estimating the storage efficiency factor and the associated uncertainty. For offshore storage systems, “deep” refers to deep geological strata beneath the seabed [28]. For both onshore and offshore storage systems, the depth of the storage reservoir should be deep enough to keep CO₂ at the supercritical phase, which means that the pressure and the temperature at the depth of the storage reservoir should be higher than 7.39 MPa and 31.1 °C, respectively. The effect of the water column above can be accounted for in calculating the densities of pore fluids (i.e., CO₂ and brine) by the in situ pressure and temperature measurements. Although marine sediments differ from sedimentary basins on land in several aspects, such as the salinity and chemistry of pore fluids, the permeability of formations [28], which result in different ranges and distributions for the input parameters in estimating the storage efficiency factor and the associated distribution, the estimation methodology is the same for onshore and offshore saline storage systems. With proper distributions specified for uncertain input parameters, the proposed workflow could be used for rapid estimation and quantification of the uncertainty in the storage efficiency factor for onshore and offshore saline storage systems.

7. Conclusions

A workflow was developed for rapid estimation and quantification of the uncertainty in the storage efficiency factor in the DOE-NETL volumetric method, with the maximum allowable pressure buildup for safe injection taken into consideration. The physical properties of storage formations and injection parameters affecting the storage efficiency factor are identified. A publicly available dataset for a prospective CO₂ storage site was utilized to demonstrate the workflow. Monte Carlo uncertainty analysis was performed using 10,000 simulations of the storage efficiency factor by randomly sampling from specified distributions for 13 independent uncertain input variables. The statistical distributions of the storage efficiency factor for the primary storage formation and the secondary storage formation located in deeper depth are derived. By computing the correlation coefficient between the output and each input uncertain variable (i.e., a physical property of the storage formation), important physical properties that contribute most to the overall uncertainty in the storage efficiency factor were identified. Our results indicate that for the prospective CO₂ storage site used as an example in this study, the effective-to-total porosity contributes most to the overall uncertainty in the storage efficiency factor, followed by the maximum allowable pressure buildup, the net-to-gross thickness ratio, the irreducible water saturation, and the permeability. Collecting more actual measurements for the identified significant uncertain input variables is strongly recommended. Reducing the degree of variability in the influential uncertain input variables will reduce the overall uncertainty in the storage efficiency factor estimates. Although the identified important uncertain input variables are specific to the study site, the proposed statistical approach can be applied to other storage sites as well. The influential uncertain physical properties of the storage formation identified from the workflow can provide guidance on future data collection needs for uncertainty reduction, improving the confidence in the CO₂ saline storage resource estimates.