Estimation of dynamic geologic CO2 storage resources in the Illinois Basin, including effects of brine extraction, anisotropy, and hydrogeologic heterogeneity

Since the vast majority of carbon dioxide (CO₂) storage resources in the United States are in deep saline aquifers, optimizing the use of these saline storage resources could be crucial for efficient development of geologic CO₂ storage (GCS) resources and basin- or larger-scale deployment of GCS in the country. Maximum CO₂ injection rates can be enhanced by extracting brine from the CO₂ storage unit. However, disposal of the extracted brine is both a technological and economic challenge. The lowest-cost option would likely be reinjection of the extracted brine into another formation above or below the CO₂ storage unit. Therefore, it is important to estimate brine injectivity as it will constrain the potential to increase CO₂ injectivity at an injection site that has access to multiple geologic storage units where either CO₂ or brine can be injected. Using a simulation-optimization framework, coupled with a non-isothermal, multiphase CO₂-water-salt equation-of-state module, we developed a computationally efficient method for evaluating optimization of simultaneous CO₂ injection, brine extraction, and brine (re)injection at hypothetical injection sites deployed across a geologic basin. The Illinois basin is ideal for testing our methodology because it contains multiple geologic storage units with seals in between them to isolate injection of CO₂ in one unit from interfering with the injection of either brine or CO₂ in another unit above or below it. In addition, we investigated the relative effects of variation in key geologic parameters as well as two reservoir structures (hydrogeologic heterogeneity/anisotropy and homogeneity/isotropy) on CO₂ injectivities and enhancement of CO₂ injectivity through extracting brine. Results suggest that permeability, depth, and especially thickness of the storage unit could be the most influential parameters determining CO₂ injectivity. They also suggest that only injecting CO₂ into the storage unit with the greatest injectivity, enhancing that unit’s injectivity by extracting brine, and disposing of the produced brine in other suitable units could maximize total CO₂ injectivity in limited regions of the basin. At the majority of simulated injection sites, however, we found that injecting CO₂ into all of the accessible and suitable storage units was more likely to maximize the CO₂ storage resource.

Introduction

A number of studies have suggested that carbon capture, utilization, and storage (CCUS) is expected to play a large role in efforts to reduce greenhouse gas emissions (e.g., National Petroleum Council, 2019; Intergovernmental Panel on Climate Change, 2023). Including projects under development in 2020, the International Energy Agency (2021) projected that 40 million metric tons per year (Mt/yr) of carbon dioxide (CO₂) could be captured, and that this figure would need to rise to 1,670 Mt/yr by 2030 and to 7,600 Mt/yr by 2050 in order for the world to be on a cost-effective pathway to net-zero emissions. If carbon capture rises to these levels, it would far exceed the current capacity for economic utilization of it (including for CO₂-enhanced oil recovery), and the vast majority would need to be stored in underground reservoirs through a method known as geologic carbon sequestration (GCS). The U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team (2013) and the U.S. Department of Energy (2015) have assessed the potential to store captured CO₂ in the Earth’s subsurface. These studies have estimated large and widespread volumetric storage resources for GCS in the United States, mostly within deep saline aquifers. However, there was no consideration of the time dimension in these assessments, and modeling of how fast CO₂ can be injected into the assessed volumes of rock pore space at basin-scale is still needed to enable optimal utilization of these CO₂ storage resources.

Since the costs and expected revenues from a CO₂ storage operation depend on the rate at which CO₂ can be safely injected into a given unit of subsurface geologic media (i.e., the dynamic CO₂ storage capacity, or injectivity), the estimated volume of rock pore space available for storage is not sufficient for assessing the economic feasibility of deployment of GCS ‘at scale’ (Anderson and Jahediesfanjani, 2019; 2020; National Petroleum Council, 2019). With respect to estimating “safe” rates of injection, a critical issue to consider is that injection into geologic media causes fluid pressure to increase within the storage unit,¹ which could increase the risks of induced seismicity (Zoback and Gorelick, 2012) or leakage of CO₂ toward the ground surface (Plampin and Merrill, 2024). Increasing the rate of injection increases reservoir pressure buildup (Plampin et al., 2023), which could increase these risks. To be able to effectively consider risk, methods developed to estimate injectivity need to control for the in-situ pressurization that is likely to occur when CO₂ is injected. In addition to managing risk, CO₂ storage operators may also need to limit the rate of injection and pressure buildup to avoid (pressure) interference with neighboring utilization of the subsurface (Anderson, 2017; Duggan Jr et al., 2024; Gasda et al., 2024; Gresham et al., 2010; Gresham and Anderson, 2010).

Pressure buildup can be mitigated by extraction of native fluids (in the case of saline aquifers, brine) from the CO₂ storage unit, thus increasing effective injectivity. However, the produced fluid must then be managed, usually by reinjection into a subsurface formation that is isolated from the CO₂ storage unit. Since brine transportation on the ground surface could be cost-prohibitive, reinjection of extracted brine would likely need to occur into a different unit that is accessible at the same geographic location (i.e., at the same well pad). Therefore, the maximum possible rate at which CO₂ could be sequestered at any given location could depend not only on the injectivity of the selected CO₂ storage unit, but also on that of other subsurface formations above or below it that are suitable for brine injection.

Some previous studies have developed numerical and analytical methods to estimate CO₂ injectivity in saline aquifers (e.g., Birkholzer and Zhou, 2009; Bandilla et al., 2012; Anderson and Jahediesfanjani, 2019; 2020), but to our knowledge, only Plampin et al. (2023) have rigorously quantified the impacts of potential constraints on enhancing basin-scale CO₂ injectivity that could be imposed by limitations on brine injectivity in sedimentary layers above or below the CO₂ storage unit. Once this important constraint is quantified, it is possible to optimize CO₂ storage resources across multiple storage units that are accessible from individual injection sites. However, this optimization requires assessment of whether more CO₂ can be injected over a given period of time if there is simultaneous injection of only CO₂ into all accessible and suitable storage units, or whether total CO₂ storage resources would be maximized by focusing injection on the unit(s) with the highest injectivity, further enhancing the injectivity of the selected storage unit(s) by extracting in-situ brine, and disposing of the produced brine onsite by reinjecting it into the other accessible and suitable units that are not used for CO₂ storage.

In this study, we have developed fully numerical simulation methods that can be used to solve this CO₂ storage resource optimization problem, and we apply these methods to analyze hypothetical deployment of 175 injection sites where either CO₂ or brine (produced to manage pressure in selected CO₂ storage units) can be injected into multiple subsurface formations. We expand upon methods developed by Plampin et al. (2023) by quantifying the constraints and impacts on CO₂ storage without relying on empirical equations that represent fewer site-specific conditions (Jahediesfanjani et al., 2019). We then compare the results of applying the more scientifically rigorous approach developed in this study to results based on empirical equations. In addition, this study applies correlation analysis to determine which input parameters are most influential on injectivity, using a combination of rank correlation and spatial cross-correlation analyses (Chen, 2015).

Materials and methods

The Illinois basin has a history of GCS that dates back more than a decade. The Illinois Basin - Decatur Project captured CO₂ from a fuel ethanol production plant and injected it into a deep saline aquifer at a rate of 1,000 tons/day from 2011 to 2014 (Finley, 2014). A second injection project (known as the Illinois Industrial Carbon Capture and Storage Project) began in 2017, and site closure is not expected until later in 2025 (Greenberg, 2021). Both projects utilize(d) the Mount Simon Sandstone as the CO₂ injection reservoir. However, several shallower subsurface formations exist within the Illinois basin that could also potentially be used for CO₂ sequestration, including the Potosi Dolomite (Leetaru et al., 2014) and the St. Peter Sandstone (Will et al., 2014).

Hydrogeological setting delineation

Figure 1 displays the relative locations of the Mount Simon Sandstone, Potosi Dolomite, St. Peter Sandstone, and Silurian-Devonian storage units (shaded in blue or yellow) as well as the sealing formations (shaded in brown or green) within an idealized stratigraphic column. Sealing formations exist above all the potential CO₂ storage units considered in this study, isolating them from each other and even-shallower freshwater aquifers. Note that the Mount Simon Sandstone, Potosi Dolomite, and St. Peter Sandstone storage units are single formations, but the Silurian-Devonian storage unit consists of multiple formations that are not isolated from each other by any intervening sedimentary layers with low enough permeability to be considered a sealing formation.

Figure 1

Figure 1. Section of a representative stratigraphic column of the Illinois basin, based on Swezey (2009). Green and brown shading indicates sealing formations, and the four permeable units that were investigated as possible CO₂ storage or produced-brine disposal units are colored in yellow for the sandstone units (e.g., Mount Simon Sandstone) or blue for the calcareous units (e.g., Potosi Dolomite). Abbreviations used: Fm., Formation; Ls., limestone; Mbr., Member; and Ss., Sandstone.

Plampin et al. (2023) mapped volumes within four permeable subsurface units that meet the conditions for potential injection of CO₂ or brine. All of these volumes were limited to being within the lateral geologic boundaries of the Mount Simon Sandstone CO₂ Storage Assessment Unit (SAU), as defined by the U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team (2013). The viable injection volumes were also limited by the existence of an adequate low-permeability sealing formation above the potential injection volume, and the non-existence of freshwater within the storage unit at the injection location. The salinity limit for freshwater was set to 10,000 ppm of total dissolved solids (TDS) in the formation waters, as codified by the U.S. Environmental Protection Agency (2014). When considering CO₂ injection, the storage units were further limited to locations where the top depth of the unit was at least 3,000 ft, because the typical temperatures and pressures at this depth and below are sufficient to maintain the injected CO₂ in a supercritical state and maximize the storage resource per volume injected (Brennan et al., 2010).

The outlines of the lateral boundaries of all four storage units are displayed in Figure 2. In the legend of Figure 2 and hereafter in the text, “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, “Silurian-Devonian” refers to the Silurian-Devonian, and “St. Peter” refers to the St. Peter Sandstone storage units. The outline of the volume within the Mt. Simon that Anderson and Jahediesfanjani (2019) defined as a potential “sweet spot” for CO₂ storage is also displayed in Figure 2. It is useful to note that within the Mt. Simon, the volume that could potentially accommodate enhanced CO₂ injection (i.e., extraction of brine and reinjection into overlying units) was limited to the region outlined by whichever boundary of the Potosi, St. Peter, or Silurian-Devonian storage units was outermost, because reinjection of brine extracted from the Mt. Simon into another unit could not occur where there are no other units that meet the criteria for brine injection and (vertically) overlap with the Mt. Simon. For the Potosi, St. Peter, or Silurian-Devonian, the entire areal extent of the storage unit overlaps with at least one of the other units that meet the criteria for reinjection of extracted brine.

Figure 2

Figure 2. Outlines of the shallower injection units investigated in this study (colored lines) are superimposed on the shaded area, which represents the extent of the Mount Simon Sandstone Storage Assessment Unit from the USGS CO₂ storage resource assessment (U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team, 2013). “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone. Within the Mt. Simon, the outline of a potential CO₂ storage “sweet spot”, as defined by Anderson and Jahediesfanjani (2019), is also displayed. Dots represent modeled injection-well locations with a uniform 25 km × 25 km well-spacing (from Plampin et al., 2023).

At every point (hypothetical injection well) on the 25 km × 25 km well-spacing grid (shown by the black dots in Figure 2), Plampin et al. (2023) estimated the depth and thickness of each of the four storage units. This hydrostratigraphic architecture was reused in this study and is displayed in Figure 3. Top depth and thickness were the main parameters for which spatially differentiated data rasters were available and varied between individual simulations. So, the spatial distributions of these parameters within each storage unit are displayed in Figure 3.

Figure 3

Figure 3. Contour plots of depth to top (left column) and total thickness (right column) of the four units, as approximated on a 25 km × 25 km grid. Each plot has different contour intervals, and the definitions of these intervals are provided by the colored numbers. “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone (adapted from Plampin et al., 2023).

Modeling domain setup

The numerical simulation methodology developed in this study builds upon methods developed by Plampin et al. (2023). The major advance in this study is that we performed direct numerical simulations to estimate CO₂ injection rates and the impacts of brine extraction from the storage units on those rates, rather than relying on empirical equations based on estimated injectivity across multiple formations and basins. To do this, we defined adjacent square “injection sites” surrounding hypothetical injection wells mapped on a 25 km × 25 km grid that roughly covered the areal extent of the Mount Simon Sandstone SAU (U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team, 2013). At each individual site throughout the basin, we modeled a well injecting CO₂ into any of the four units that were accessible and met the criteria for CO₂ sequestration at that location. The simulations estimated maximum injection rates that could be achieved without violating limitations on pressure buildup, including not allowing pressure interference with neighboring injection sites. Figure 4 displays the general modeling conceptualization of six injection sites that penetrate a hypothetical subsurface layer that dips and thickens from left to right. Each injection volume x is represented by a square prism that measures 25 km × 25 km (A_x = 125 km²) in the lateral dimension, while Z_top,x and B_x represent the depth to the top and thickness, respectively, of injection volume x. Injection was simulated at the center of each square prism, and limitations on pressure buildup were imposed at both the injection location (i.e., bottomhole pressure at the well, P_well,x) and the nearest boundary of the prism (P_boundary,x) in the simulations.

Figure 4

Figure 4. Schematic showing the general modeling conceptualization through a set of six injection volumes in a hypothetical subsurface layer that dips and thickens from left to right. The qualitative dip is exaggerated in this schematic. Variables Z_top, B, P_well, P_boundary, and A indicate the top depth, thickness, bottomhole well pressure, outer boundary pressure, and areal extent of the hypothetical injection sites indicated by subscripts 1 through 6 (adapted from Plampin et al., 2023).

Model parameterization

The model domain was parameterized first with a set of homogeneous/isotropic rock properties of porosity and permeability associated with each storage unit considered. The homogeneous parameter values of permeability and porosity for each unit were assumed to be equal to those defined in Plampin et al. (2023), and these are summarized at the top of Table 1. Below that summary of parameter values, the table provides the number of injection volumes within each unit that met the criteria for the purpose listed in the second column. Since all regions of the storage units we modeled met the criteria for CO₂ or brine injection, the number of volumes listed for each purpose were mostly equal to the total number of injection wells penetrating each unit. The only exception was the Mt. Simon, because there were a number of sites on the periphery of this storage unit that were located outside of the lateral boundaries of any of the other units (Figure 2). Under the assumptions of our model, onsite disposal of brine produced from the Mt. Simon would not be possible at these locations, which prohibited brine extraction to enhance CO₂ injection.

Table 1

Table 1. Constant (unit average) values^a of porosity and permeability that were applied in the homogeneous/isotropic simulations. “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone. The number of injection wells that penetrated each storage unit were tallied on the 25 km × 25 km grid shown in Figure 2 Φ is porosity, k_h is horizontal permeability, and B_x is unit thickness.

We first performed simulations to estimate CO₂ injectivity. This first round of simulations is referred to as the “base case”, because no injectivity enhancement via brine extraction was simulated. The alternative case, “enhanced CO₂ injectivity” was analyzed via simulations at geographic locations where CO₂ could be injected into at least one storage unit and brine could be injected into at least one other unit that was accessible at the same location. Multiphase flow parameters 1/P₀, λ, S_lr, S_ls, and S_gr, which represent inverse air entry pressure, a retention function shape parameter, residual liquid saturation, saturated liquid saturation, and residual gas saturation, respectively, were used in the constitutive models of van Genuchten (1980) and Mualem (1976). These parameters were uniformly set to 5.10 × 10⁻⁵, 0.457, 0.0, 1.0, and 0.05, respectively, in all simulations. This was done to limit the number of independent variables that could affect the results, which allowed the analysis of the potential impacts of variation in geologic parameters on CO₂ injectivity to be more tractable.

In addition to developing an improved numerical simulation approach to estimate CO₂ injection rates, a main goal of this study was to investigate the effects of hydrogeologic heterogeneity and anisotropy on injectivity. Unlike for depth and thickness, spatially differentiated data rasters, published structure contours, or detailed data were not available for estimating localized parameter values of permeability and porosity. For the heterogeneous/anisotropic simulations, localized estimates of porosity (Φ) were obtained based on equations from the literature that relate porosity to depth (z), and then these calculated porosities were used as inputs into other equations from the literature that relate permeability to porosity in order to estimate the associated permeabilities in those same locations. These equations and the sources of them are provided in Table 2. The carbonate properties ($\mathit{\Phi }$_carbonate and k_h^carb) were used for both the Potosi and the Silurian-Devonian, because both of these storage units are composed of carbonate rocks (Casey, 1996). Anisotropy was applied by using a constant ratio to convert the heterogeneous horizontal permeabilities (k_h) into the associated vertical permeabilities (k_v) within the same column beneath each 25 km × 25 km injection site. Thus, the parameter values for porosity and permeability that were used as inputs in the heterogeneous/anisotropic simulations varied between the injection volumes according to their different depths.

Table 2

Table 2. Equations used in the parameterization of the heterogeneous/anisotropic simulations. “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone. Variables in the equations are porosity (Φ), depth (z), horizontal permeability (k_h), and vertical permeability (k_v).

Numerical simulation methodology

To increase computational efficiency, we took advantage of symmetry by using a computational mesh that represented one-eighth of the volume of a (25 km × 25 km × B_x) square prism like those depicted in the schematic in Figure 4. This mesh, which takes the shape of an extruded right triangle, is shown in more detail in Figure 5. The mesh contained 10 rows of grid blocks, each of which was set to a thickness equal to one-tenth of the unit thickness (B_x) at the hypothetical injection site being modeled (see Figure 3 for a geographical depiction of formation thicknesses across the basin). The simulations involved fifty years of CO₂ injection into one of the bottom corners of a symmetric cell that formed a 45-degree angle in the view from the top of the cell (Figure 5). In simulations of enhanced injectivity, fifty years of brine extraction was assumed to be from the opposite-bottom corner of the symmetric cell. The simulations monitored pressure-buildup at the CO₂ injection location as well as at the corner of the grid block domain that takes the shape of a right angle in the view from the top of the cell (i.e., the corner of the grid block that takes the shape of a square prism in Figure 5).

Figure 5

Figure 5. Computational mesh used in the numerical simulations. Using symmetry, we modeled every 25 km × 25 km injection domain as an eighth of its actual volume. In all three directions, the domain was discretized into ten evenly distributed grid blocks. CO₂ injection was simulated into the bottom cell at one corner of the symmetric cell. In simulations of enhanced CO₂ injectivity, brine extraction was simulated from the opposite-bottom corner of the domain.

For every modeled injection site, we assigned the depth and thickness of each storage unit accessible at that geographic location according to Figure 3. Depth was incorporated into the simulations by assigning an initial pressure at the top of the domain equal to 10,000 Pa per meter (P_a/m) times the unit top depth (Z_top) at that location, and the same hydrostatic gradient through all the rows of grid blocks below the top row. The maximum pressure condition applied at the injection location was set to 80% of the fracture pressure, or 12,900 P_a/m times the depth (z) of the bottom row of grid blocks. Hydrologic properties were assigned uniformly to all the simulated injection volumes within a given storage unit in the homogeneous/isotropic simulations.

The injection rate was constant over the entire 50-year injection period. The simulation algorithm determined the maximum injection that would not violate any of the limiting conditions on pressure buildup. “No-flow” boundary conditions were imposed in all simulations. In the case of the top and bottom boundaries on the modeled domain, this assumption was motivated by the units targeted for injection being significantly more permeable than the units above or below them. In the case of the outer (lateral) boundaries, the no-flow boundary conditions were motivated by the symmetry of the well configuration, combined with the condition that only minimal pressure interference between neighboring injection units was allowed. This condition was assured by imposing the pressure constraint at that boundary and limiting injection rates at the center of each storage subdomain to not violate that constraint. Specifically, we allowed a maximum pressure increase of only 100 kilopascals (kP_a) at the boundary of the injection volume nearest the injection location. If operators deploy wells at a spacing that is different from the generic, regular configuration assumed here, they could adjust injection rates to be higher than those estimated in this study for wider well-spacings or lower for tighter well-spacings in order to maximize the rate of injection at each well without violating the same pressure constraints we have imposed.

In the base-case simulations for this study, we assumed that only CO₂ would be injected into all four potential storage units being considered. This was done to demonstrate the maximum amount of CO₂ that could be injected without performing any pressure management via brine extraction. For the enhanced cases, brine extraction was simulated from the bottom grid block (as depicted in Figure 5), and along the entire thickness of the model domain. This attempt to demonstrate the effect of different extraction-well-screening strategies was partially motivated by a finding by Plampin and Merrill (2024) that (in the case of CO₂ leakage into shallower aquifers) hydrological plume mitigation by pumping can be affected by the length of the screened interval of the extraction well. In simulations that included brine extraction, the total rate of brine extraction from the CO₂ storage unit was set to the sum of the maximum brine injection rates estimated for all of the other storage units that were accessible at the same geographic location of the injection site (Plampin et al., 2023).

In summary, we defined scenarios that included a base case (without brine extraction) and an enhanced case (with brine extraction); we considered two reservoir structures (homogeneous/isotropic and heterogeneous/anisotropic); and for the enhanced cases, we considered extraction of brine from the bottom vs brine screening throughout the total thickness of the storage unit. For each injection volume and case (base or enhanced), and each reservoir structure (homogeneous/isotropic or heterogeneous/anisotropic), we performed an optimization simulation to estimate the maximum CO₂ injection rate. To do this, we used the multiphase, multicomponent flow and transport simulator iTOUGH2-ECO2N (Pruess, 2005; Finsterle, 2007; Finsterle et al., 2017). The parameter estimation capabilities of this tool were re-purposed to automatically adjust the CO₂ injection rate to the maximum (constant) value it could achieve without violating the pressure buildup constraints listed above at any time during the 50-year injection period. The applicability of the two types of reservoir structure, and consideration of the two cases (with or without brine extraction) and the two brine extraction strategies for the simulations of all of the storage sites penetrating all four of the storage formations are provided in Table 3. The simulation numbers in the table correspond to the subfolders in the U.S. Geological Survey (USGS) model archive data release that accompanies this paper (Plampin et al., 2025).

Table 3

Table 3. Scenarios for the multi-optimization simulations. “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone. Each number in the interior of the table identifies the simulation that corresponded to a scenario that considered one of the four storage units, one of the two reservoir structures, one case (with or without brine extraction), and (for the enhanced cases) one of the two brine extraction strategies shown in the column and row headers.

Benchmark empirical injectivity calculations

For all homogeneous/isotropic simulations, the base-case and enhanced CO₂ injection rates estimated using the comprehensive numerical simulation methods developed in this study were compared to benchmark estimates of basin-scale injectivity based on applying empirical relationships developed by Jahediesfanjani et al. (2019). To do this, we first estimated base-case CO₂ injectivities for every storage unit that vertically overlapped each injection location using Equation 1:

$V_{inj,t}=a{\varphi }^b{H}^c{A}^d{k}^e{D}^f{T}^g(1)$

In Equation 1, V_inj,t is the maximum CO₂ injection rate, Φ is porosity (in percent, %), A is the areal extent of the site around each hypothetical injection well (a uniform 125 km² for all injection sites in this study), D is the depth to the top of a possible CO₂ storage unit (in meters, m), H is the net thickness of the storage unit (in m), k is the reservoir permeability (in millidarcies, mD), T is the injection duration (consistently set at fifty years in this study), and a, b, c, d, e, f, and g are dimensionless curve-fitting constants. Based on data from eleven separate CO₂ storage units (including some in the Illinois basin), Jahediesfanjani et al. (2019) estimated a, b, c, d, e, f, and g to be equal to 1.67 × 10⁻⁹, 0.66, 0.85, 0.64, 0.18, 1.25, and −0.31, respectively.

We then estimated enhanced CO₂ injectivities using Equation 2:

$V_{inj,t}^{ *}=\frac{1.53-V_{inj,t}}{25,000}\left(V_{ext,t}-\frac{25,000\times V_{inj,t}}{V_{inj,t}-1.53}\right)(2)$

In Equation 2, $V_{inj,t}^{*}$ is enhanced CO₂ injectivity and V_ext,t is the associated brine extraction rate, which was set to the sum of the brine injectivities for the (maximum of three) other (non-CO₂ storage) units accessible at that same location, as estimated by Plampin et al. (2023). Finally, we calculated “enhancement” as the percent increase from the base-case CO₂ injection rate to the enhanced rate at each site location in each scenario.

Correlation analysis

Rather than fitting a predictive relationship as in Jahediesfanjani et al. (2019), we estimated correlation coefficients to determine which model input parameters were more or less influential on injectivity. We tested the significance of the individual (bivariate) correlations between each of the input parameters (initial pressure, thickness, porosity, and horizontal permeability) and the corresponding CO₂ injection rate. We tested these individually for each storage unit, accounting for reservoir structure (either homogeneous/isotropic or heterogeneous/anisotropic), as well as aggregated over stacked formations that were accessible at the same well-site and both reservoir structures. All correlation coefficients (e.g., Pearson, Kendall, or spatial/partial correlations) take values between −1 and one and describe the sign and magnitude of the correlation between two random variables. The Pearson correlation coefficient may not be appropriate for our analysis because: 1) it assesses the degree of linearity between the two variables, and 2) it does not account for spatial autocorrelation. Given the complexity of the reservoir fluid flow dynamics approximated by the numerical methods used in this paper, a linear relationship between the input parameters and estimated injection rates could be considered unlikely. Still, we constructed scatterplots and used them to provide visual evidence of whether there might be some degree of a linear relationship between estimated injection rates and either depth or thickness. In addition, we estimated Kendall rank correlations, which is an alternative measure that can capture nonlinear relationships between variables (Kendall, 1938).

However, we mostly focused on determining whether the input parameters and injection rates estimated in this study could be spatially auto- or cross-correlated. Spatial autocorrelation is correlation of a variable at a given location with the same variable at other (nearby) locations, whereas spatial cross-correlation quantifies the correlation between one spatially located random variable and another at a collocated or proximal location. Spatial correlation can be interpreted visually using bubble plots/heatmaps, and these visual interpretations can be supported by the use of geostatistical (spatially autocorrelated) models to predict/simulate both the thickness and depth of geological formations for resource assessments (Olea and Luppens, 2014). Methods to account for spatial relationships when estimating the cross-correlation between two spatially indexed variables have been developed, and these methods often rely on decomposing the covariance matrix into a spatial and non-spatial (partial) component (Chen, 2015). The simple correlation coefficient (SCC), R₀, is the sum of the spatial cross-correlation index (SCI), R_c, and the partial spatial cross-correlation coefficient (PSCC), R_p. Separating out the PSCC component from the SCC is of specific interest in this study because it isolates the correlations between the input variables and injectivity from the spatially autocorrelated component. In order to separate out the PSCC component of the SCC, we used a spatial-weights (adjacency) matrix with first-order neighbors at every point (hypothetical injection well location) on the 25 km × 25 km well-spacing grid (Figure 2).

In summary, we approached the question of which input parameters influence injectivity using an analysis of rank correlations, and we combined this with an analysis of partial spatial correlations. We then used permutation tests to calculate p-values for each correlation coefficient, and we adjusted the final p-values using the Bonferroni correction for the number of correlation tests performed (Bonferroni, 1936; Miller, 1991). All correlation analyses were implemented in R version 4.4.1 using the tidyverse package (R Core Team, 2024; Wickham et al., 2019).

Results

Our first goal was to estimate CO₂ injectivities across all four storage units individually, as well as the enhancement that could potentially be achieved by extracting brine from a single CO₂ storage unit and injecting it into all of the other units accessible above or below the storage unit at that geographic location (i.e., the enhanced-injectivity cases). The results from this analysis are displayed graphically, and ranges of the numerical results are provided in Figure 6. For enhanced-injectivity cases, we only present the results in which brine extraction was simulated from the bottom cell of the model domain (as shown in Figure 5), because there was no significant difference between these results and the results of simulating brine extraction from the entire thickness of the CO₂ storage unit.

Figure 6

Figure 6. CO₂ injectivities estimated by the simulation-based optimization algorithm. Each row of graphs presents results for CO₂ injection into one of the four storage units considered. “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone. The ranges and spatial distribution of base-case injection rates (without enhancement via brine extraction), the percent increase in base-case injection rates caused by extracting brine from the CO₂ storage unit at the maximum rate determined by the feasible rate of reinjecting it into all of the other accessible units, and the final enhanced CO₂ injection rates are displayed, respectively, in plots (a–c) for the Silurian-Devonian; (d–f) for the St. Peter; (g–i) for the Potosi; and (j–l) for the Mt. Simon. The totals in the bottom-right corners of plots a,d,g,j (left column) are the sums of the base-case CO₂ injection rates of all of the wells penetrating each storage unit, the totals in plots b,e,h,k (middle column) are the average percent increases in base-case injection rates across the same sets of wells; and the totals in plots c,f,i,l (right column) are the sums of the enhanced CO₂ injection rates for the same sets of wells.

The base-case CO₂ injectivities (per well) for each storage unit (displayed in the bottom-left corner of each plot in the left column of Figure 6) were positively correlated with the volume of the porous interval penetrated by the injection well. Since every injection-well site had identical areal dimensions of 25 km × 25 km, the injection volume depended entirely on unit-thickness at each site. Therefore, the overall ranking of average unit-injectivity per well was the same as the ranking of average unit-thickness, which was highest for the Mt. Simon, followed, respectively, by the Silurian-Devonian, Potosi, and St. Peter, although the St. Peter and especially the Potosi were still found to have notable base-case injectivity. In addition, the spatial distribution of injectivity within each storage unit (shown by the colored shading) appeared to be more sensitive to the thickness of the unit than the depth at each hypothetical injection location (See Figure 3 for the basis of these comparisons.).

The results presented in the middle column of Figure 6 suggest that the CO₂ injectivity of each storage unit could be significantly enhanced via brine extraction at the maximum rate constrained by injectivity into the other accessible units. The average percent enhancement per unit (the “Total” in the bottom-right corner of each plot in the middle column) varied between 12% (for the Mt. Simon) and 311% (for the St. Peter). The higher percent increases in injectivity (middle column) mostly occurred in the regions within each storage unit that had lower base-case injectivity (left column). The geospatial distributions of relatively lower or higher enhanced-injectivity regions and the final numerical ranges of enhanced injection rates are presented in the right column of Figure 6. Because of enhancement through brine extraction, the potential “CO₂ injectivity sweet spot” of the Mt. Simon shifted from being focused in the northern part of the basin toward the southern part.

The range of enhanced injectivity per well in the new sweet spot of the Mt. Simon was estimated to be about 0.4–0.5 Mt/yr CO₂ (Figure 6l) compared to the base-case injectivity in roughly the same region of the same storage unit of about 0.2–0.3 Mt/yr CO₂ (Figure 6j). Vertically above this possible enhanced-injectivity sweet spot of the Mt. Simon, the base-case injectivities per well of the other three storage units were estimated to be in the ranges of about 0.1–0.3 Mt/yr CO₂ in the Silurian-Devonian, 0 to 0.04 Mt/yr CO₂ in the St. Peter, and 0 to 0.2 Mt/yr CO₂ in the Potosi (Figures 6a,d,g, respectively). If an injection site were located where the base-case injectivities of all four storage units were at the lowest ends of the estimated ranges (displayed in Figures 6a,d,g,j), our results suggest that the total injectivity (if injecting CO₂ into all suitable units) could be about 25% less than injecting CO₂ only into the Mt. Simon at a rate of about 0.4 Mt/yr CO₂ (the lowest end of the range of enhanced injectivity in the new sweet spot in Figure 6l). However, if an injection site were located where the base-case injectivities of all four storage units were at the highest end of each range (displayed in Figures 6a,d,g,j), our results suggest that the total injectivity (if injecting CO₂ into all suitable units) could be about 60% greater than injecting CO₂ only into the Mt. Simon at a rate of about 0.5 Mt/yr CO₂ (the highest end of the range of enhanced injectivity in the new sweet spot in Figure 6l). Given the level of uncertainty in our modeling results, this rough comparison represents the extent of this type of analysis that we could perform. Still, it suggests that there are likely to be locations where total CO₂ injectivity could be maximized by injecting CO₂ only into the Mt. Simon, enhancing its injectivity as much as possible by extracting brine, and only using the other accessible units for disposal of the produced brine. Qualitatively, this is consistent with a similar result in Plampin et al. (2023), who found a limited region in the Illinois basin where total CO₂ injectivity could be optimized by injecting CO₂ only into the Mt. Simon and using the Silurian-Devonian only for disposal of brine extracted to maximize injectivity of the Mt. Simon.

Comparison with benchmark injectivity estimates

To benchmark our fully numerical simulation-based CO₂ injectivity optimization methodology against the empirical injectivity estimates by Plampin et al. (2023), we directly compared the results from the homogeneous/isotropic simulations for this study with the results from that previous paper and present the comparative results in Figure 7. In general, the numerical and empirical methodologies produced similar injectivity estimates for the base cases, as most of the data points were located near the black 1:1 line in the left plot. The discrepancies between the empirically estimated and simulated injectivities increased with increasing injection rates. For the enhancement cases (displayed in the right plot of Figure 7), the empirical injectivity estimates were generally greater than the simulated ones, as most of the data points extend out to the right of the 1:1 line. This was especially true for the injection volumes with the greatest injectivities, which were all within the Mt. Simon. We suggest that the results of the simulation-based method developed for this paper could be more representative of expected injection rates for potential CO₂ storage in the Illinois basin. This is because they were not based on applying empirical equations that were calibrated based on data and simulations for a nationwide set of potential CO₂ storage units located in multiple different (geologic types of) basins. Instead, the data and input values of the parameters used in the simulations for this study were specific to the injection locations and storage units that were considered in this analysis.

Figure 7

Figure 7. Comparisons of CO₂ injectivities estimated via numerical optimization simulations on the y-axis, versus that based on calibrated empirical relationships on the x-axis. Black diagonal lines represent parity between the two methods (i.e., a 1:1 correlation). “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone.

The relationships between the estimated injectivities and the values of storage-unit depth and thickness that were input into the base- and enhanced-case simulations for this study are summarized in Figure 8. The much tighter dispersion of data points (around a line) in the plots in the right column, as compared to those in the left column, suggests that injectivity generally exhibited a stronger linear relationship with thickness than with depth. However, for the Silurian-Devonian (displayed as red triangle markers in the plots), there appeared to be a stronger linear relationship between injectivity and depth than for the other storage units, especially in the enhanced cases. Still, it is important to recognize that depth and thickness were also closely correlated with each other in the Silurian-Devonian, which could imply that the more linear relationship between injectivity and depth of the Silurian-Devonian may be attributable to some correlation with thickness as well as depth.

Figure 8

Figure 8. Scatterplots showing the relationships between estimated CO₂ injectivities and the input values of the depth and thickness parameters. “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, and “St. Peter” refers to the St. Peter Sandstone.

Effects of heterogeneity and anisotropy

The results of the heterogeneous/anisotropic simulations of base-case CO₂ injectivities are displayed in comparison to those of the homogeneous/isotropic simulations in the scatter plot in Figure 9. For the base case in Figure 9, the estimated injectivity was typically higher in the results of the heterogeneous/anisotropic simulations than in those of the homogeneous/isotropic simulations. Since permeability has been found to be an important determinant of injectivity (Raza et al., 2016), this result (in part) is likely owing to the heterogenous values of permeability (estimated using the equations in Table 2) used as inputs into the heterogeneous/anisotropic simulations being significantly greater, on average, than the unit- or reservoir-average values of permeability (based on results of basin-scale CO₂ storage resource assessments and other large-scale studies) used in the homogeneous/isotropic simulations.

Figure 9

Figure 9. Scatter plot of base-case CO₂ injectivities based on results of the heterogeneous/anisotropic vs the homogeneous/isotropic simulations. “MtSimon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, “SilDev” refers to the Silurian-Devonian, and “StPeter” refers to the St. Peter Sandstone.

Comparison of the increases in base-case CO₂ injectivities as a result of brine extraction in the results of the heterogeneous/anisotropic simulations to those in the results of the homogeneous/isotropic simulations is displayed in the scatter plot in Figure 10. The estimated percent increases over the base-case injectivity estimates achieved through enhancing CO₂ injectivities by extracting brine to relieve pressure constraints was typically greater in the results of the homogeneous/isotropic simulations than in the results of the heterogeneous/anisotropic simulations. This can be observed in Figure 10 by noting the far greater number of points corresponding to percent enhancement results to the right of the 1:1 line. It should also be noted here that varying the well-screen length for brine extraction did not have a significant impact on enhanced CO₂ injectivity, and this could be owing to the brine extraction rates being constrained by the maximum rates of brine disposal in the other three storage units to be identical for both screen-length assumptions. Thus, we did not obtain results that were amenable to analysis of the potential effect of different brine-extraction-well-screening lengths on enhancement of CO₂ injectivities in this study.

Figure 10

Figure 10. Scatter plot of the percent enhancement of base-case CO₂ injectivities based on the results of heterogeneous/anisotropic vs homogeneous/isotropic simulations. “MtSimon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, “SilDev” refers to the Silurian-Devonian, and “StPeter” refers to the St. Peter Sandstone.

Since the results presented in Figure 9 indicated typically higher base-case injectivity estimates as a result of the heterogeneous/anisotropic simulations (which used higher input values of permeability, on average), there could have been less capacity for enhancing those already higher injectivity estimates by simulating the extraction of brine compared to the percent enhancement in the homogeneous/isotropic simulations. Without more detailed field data (from core samples) on the potential heterogeneity of permeability in these storage units, we could not determine whether our relatively higher estimates of injectivity and lower estimates of the extent that extracting brine might enhance injectivity were strictly owing to accounting for heterogeneity/anisotropy. Spatially differentiated rasters of heterogeneous permeability data were not available for the saline aquifers considered in this study. This data issue could be less of a problem if our model were applied to estimate the CO₂ injectivity of storage units for which there exist more extensive drilling core samples and permeability data.

Correlation analysis

We first tested rank correlations (displayed in Figures 11a,b) between each input variable (horizontal axis) and injectivity. Correlations were estimated for each of the four storage units individually and combined (vertical axis) and by reservoir structure (left and right columnar panels). Warm colors indicate positive correlations and cool colors negative correlations, with an asterisk indicating statistical significance (permutation test, p < 0.05 with Bonferroni correction).

Figure 11

Figure 11. Panels (a,b) provide Kendall rank correlations between input parameters (horizontal axes) and injectivity by storage unit (vertical axis) and reservoir structure (left and right panels). Panels (c,d) provide simple, partial, and spatial correlations (horizontal axes) of injectivity with thickness. In all panels, an asterisk indicates statistical significance (permutation test, p < 0.05 with Bonferroni correction). “Mt. Simon” refers to the Mount Simon Sandstone, “Potosi” refers to the Potosi Dolomite, “St. Peter” refers to the St. Peter Sandstone, and “Init. Press.” refers to Initial pressure.

In the results of the homogeneous/isotropic simulations (Figure 11a), there was a consistent positive correlation between injectivity and both initial pressure (i.e., depth) and thickness. Higher pressure could lead to a smaller specific volume of CO₂ and, thus, a smaller pressure build-up induced by injecting it into the formation, which could allow for higher injection rates without violating the maximum pressure constraints. Furthermore, a thicker storage unit with the same homogeneous characteristics as a thinner one could accommodate a higher total volume of CO₂. Correlations of injectivity with thickness were typically stronger than with initial pressure (i.e., depth).

Comparing the influence of the two different reservoir structures overall, there were weaker and insignificant correlations (p ≥ 0.05) of injectivity with porosity and initial pressure and a slightly stronger correlation with thickness in the results of the heterogeneous/anisotropic simulations (Figure 11b) than in the results of the homogeneous/isotropic simulations (Figure 11a). In the results of the heterogeneous/anisotropic simulations, the correlations of injectivity with initial pressure and with porosity for the St. Peter were of the opposite sign of that for the other storage units. However, this result should be interpreted with caution, because the St. Peter had the lowest sample size of all the four storage units considered in this study (Table 1).

For the homogeneous/isotropic simulations, correlations of porosity and permeability with injectivity were not estimated, because the values of these input parameters were assumed constant (at reservoir-average values) in this set of simulations. Parameterization of porosity and horizontal permeability (based on depth) resulted in correlations of porosity with injectivity (third columns of panels a and b in Figure 11) that were of equal magnitude but opposite sign of the correlations between initial pressure and injectivity (first columns of panels a and b in Figure 11) under both reservoir structures and for all storage units, except for the homogeneous/isotropic simulations of the injectivity of the Mt. Simon.

For the heterogeneous/anisotropic simulations, heterogenous porosity and permeability values were estimated based on depth using the equations in Table 2. Since these equations resulted in estimates of higher porosity and permeability at greater depth, a significant correlation of depth with thickness could impact the reliability of inferences concerning the relative dependence of injectivity on depth vs thickness. The mappings in Figure 3 suggest that there could be vertical overlap between some of the thicker and deeper portions of the Silurian-Devonian and the St. Peter, as well as with the thicker and shallower portions of the Mt. Simon. Upon further analysis, however, only the Silurian-Devonian had a strong and significant rank correlation between depth and thickness (r = 0.68*), and all other storage units had weaker and insignificant correlations (r ≤ 0.16) between these two input variables. The strong and significant correlation between the depth and thickness input parameters for the Silurian-Devonian suggests that the (unaccounted for) independent correlation between depth and thickness could significantly impact the results for the correlations between each of these input parameters and estimated injection rates. Therefore, the results for the correlations of injectivity with thickness and initial pressure for the Silurian-Devonian in Figure 11 should be interpreted with the significant independent correlation between depth and thickness of the storage unit in mind.

Visual inspection of the plots on the left in Figure 8 suggests that our injectivity estimates might exhibit a linear relationship with the values of the thickness input parameter, and the results in panels a and b in Figure 11 indicate that our injectivity estimates were most strongly correlated with thickness. Therefore, results just for the correlations of injectivity with thickness are presented in panels c and d of Figure 11. Panels c and d of Figure 11 display the results of the decomposition of the simple correlation coefficient between injectivity and thickness into the partial and spatial correlations between thickness and injectivity for each of the four storage units individually and combined. The significance and strength of the simple correlations (left column of panels c and d in Figure 11) are consistent with the rank correlations with thickness in panels a and b.

In panels c and d of Figure 11, the simple correlations indicate the total strength of the association between thickness and injectivity, accounting for the direct (partial) correlation without distance effect and the indirect (spatial) correlation based on the distance decay effect. In the second and third columns of the lower panels (c and d), the simple correlation coefficients are decomposed into partial and spatial correlations, respectively. In general, the partial component of the simple correlations for the Potosi and St. Peter was larger than the spatial component, and the spatial component was larger than the partial component for the Mt. Simon and Silurian-Devonian. Thus, we interpret the strong, simple correlation between thickness and injectivity as a mostly spatial (indirect) correlation in the Mt. Simon and Silurian-Devonian, and a mostly non-spatial (direct) correlation in the St. Peter and Potosi. These results indicate to what degree injection rates may covary spatially and directly with thickness, and therefore spatial/partial decompositions such as shown in Figures 11c,d provide more detail when analyzing cross-correlations between spatially indexed data compared to traditional (Pearson) correlations.

Discussion

Estimating the maximum feasible injection rate at any selected geographic location can enable optimal utilization of subsurface CO₂ storage resources. Pressure buildup caused by CO₂ or brine injection could impact feasibility, as excessive pressurization could increase risks such as induced seismicity or leakage of CO₂ or brine from the storage unit. In the case of basin-scale deployment of geologic carbon sequestration (GCS), closely spaced injection wells could also interfere with one another via pressurization of a shared storage unit. Pressure buildup can be managed by extraction of native brine from the storage unit, but the brine would need to be disposed of, most likely by injection into a different unit at the same geographic location but isolated by at least one low-permeability layer from the CO₂ storage unit. Therefore, brine injectivity into other formations above or below the CO₂ storage unit could constrain maximum feasible CO₂ injection rates. In order to find optimal injection strategies that deal with these complicated circumstances, we developed an approach that provides insight into whether total CO₂ injectivity could be maximized by injecting CO₂ into a single storage unit, extracting brine to maximize injectivity in the selected unit without violating limitations on pressure buildup, and reinject the produced brine into other suitable units (also without violating pressure limits), or whether total CO₂ injectivity could be maximized by simply injecting CO₂ into all accessible and suitable storage units.

Using a small number of input parameters, we applied our new methodology to numerically estimate CO₂ injectivity in the Illinois basin. This was a proof-of-concept study, and a greater number of input parameters could be incorporated into this framework in future investigations. A uniform well-spacing on a 25 km × 25 km grid across the areal extent of possible storage units was assumed. For each hypothetical injection site, a computational mesh was used to model a representative injection pressurized volume created by each injection well. A simultaneous fifty years of CO₂ injection was simulated at each hypothetical injection site, and pressure was monitored at both the bottom of the injection well and the nearest boundary of the domain (representing the boundary between two neighboring injection sites). Inverse modeling capabilities were used to automatically adjust the CO₂ injection rate in the simulations to achieve the maximum value that would not cause exceedance of pressure limitations at the monitoring locations for each site.

Based on the same uniform well-spacing for four potential storage units in the Illinois basin, we first performed “base-case” CO₂ injection simulations without brine extraction. We then estimated “enhanced” CO₂ injectivities for each of the four storage units by performing optimization simulations that included extraction of in-situ brine from a selected CO₂ storage unit at a rate that equaled the maximum brine injectivity into the remaining storage units that were accessible (onsite) at that geographic location. We also benchmarked our new numerical optimization simulation method against a previously developed empirical method for calculating injectivity and found comparable results. The new method developed here is more rigorous than the empirical method, in that it performs numerical simulations of multiphase flow in the subsurface for the purposes of injectivity estimation. More complex hydrogeological architecture can be analyzed using this numerical simulation domain than by using the empirical method, allowing for more localized injectivity estimates.

Our approach allowed us to analyze and compare the correlations of geologic characteristics and evaluate their relative importance for estimating injectivity. Our correlation analysis suggests that CO₂ injectivity depends strongly on the thickness of the storage unit, and to a lesser degree on the depth, porosity, and permeability. Significantly higher base-case CO₂ injection rates (without brine extraction) were estimated as a result of the heterogeneous/anisotropic simulations compared with those as a result of the homogeneous/isotropic simulations. A related finding was that the incremental enhancement (via extraction of brine) of the already higher base-case CO₂ injection rates estimated as a result of the heterogeneous/anisotropic simulations was typically less than the percent enhancement of the relatively lower base-case injection rates estimated as a result of the homogeneous/isotropic simulations. Further analysis revealed that these results may have been partially owing to the permeabilities estimated using the equations in Table 2 and input into the heterogeneous/anisotropic simulations being greater, on average, than the reported reservoir- or unit-average values of permeability input into the homogeneous/isotropic simulations. If field data (from drill-core samples) is available to estimate heterogeneous values of permeability in saline aquifers at basin-scale, the methods developed in this study could be used to better estimate whether there could be relatively lower base-case CO₂ injectivity and greater potential for enhancement through brine extraction in homogeneous/isotropic storage units vs in more heterogeneous/anisotropic units, or whether this difference in reservoir structure is not likely to significantly impact CO₂ injectivity.

The results in Figure 6 indicated that brine extraction from any of the CO₂ storage units significantly increased its total injectivity, and the mappings also suggest that brine extraction could expand and/or shift the spatial distribution of the highest injectivity zones within each unit. Depending on the degree of overlap of the highest enhanced injectivity zones in one storage unit with lower base-case injectivity zones in the other three units, the results suggest that there may be some limited regions in the Illinois basin where maximum CO₂ injectivity could be obtained by injecting CO₂ only into the storage unit with the highest base-case injectivity (the Mt. Simon), extracting brine from the Mt. Simon to maximize its injectivity, and disposing of the produced brine in the other accessible units. Follow-on economic analysis of the results of this study could be used to determine whether the value of being able to inject CO₂ at the enhanced rates predicted by our model would be greater than the additional costs to extract brine from the targeted CO₂ storage unit and reinject it for disposal in another unit that is accessible at a given injection site location. If present in sufficient concentrations, part of the economic incentives to produce brine to manage pressure during CO₂ storage operations could include the opportunity to extract valuable mineral commodities from the brine before reinjecting it for disposal. At most of the hypothetical injection sites simulated in this study, however, estimated CO₂ injectivity was maximized by injecting CO₂ into all accessible storage units without any brine extraction/disposal.

The methods developed in this paper could also be used to analyze situations where brine needs to be extracted purely to manage risk or it may not be feasible to inject CO₂ into all accessible units at a given injection site. If it is more difficult to obtain a permit or more costly to inject CO₂ into some of the potential storage units being considered than in others, if brine produced to enhance CO₂ storage meets the criteria for disposal in the units that are not feasible for CO₂ storage, and the units overlap at a selected injection location, then this could contribute to focusing on development of just one storage unit to sequester CO₂, extracting brine to enhance the injectivity of the selected CO₂ storage unit, and only developing the other accessible units to dispose of the produced brine.

Data availability statement

The model input and output files generated for this study are available as a U.S. Geological Survey (USGS) Model Archive Data Release (Plampin et al., 2025). Some data used to construct the numerical models for this study were not gathered by the USGS and may have limited external availability because the private entities that compiled the data have a proprietary interest in them. For more information on the sources of these data, please contact the USGS Geology, Energy & Minerals Science Center in Reston, VA.

Author contributions

MP: Conceptualization, Investigation, Methodology, Data curation, Formal Analysis, Visualization, Writing – original draft. SA: Conceptualization, Formal Analysis, Supervision, Visualization, Writing – original draft, Writing – review and editing. SF: Software, Methodology, Writing – review and editing. AW: Formal Analysis, Visualization, Writing – original draft, Writing – review and editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Acknowledgments

Acknowledgements

We acknowledge the helpful comments provided on this product by Emil Attanasi, Matthew Merrill, Jeremy Ray, Karl Haase, Jenna Shelton, Ole Kaven, and Clinton Scott of the U.S. Geological Survey and two reviewers selected by the journal. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflict of interest

Author SF was employed by Finsterle GeoConsulting, LLC.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹In this study, a storage unit is defined as a mappable volume of rock consisting of a storage formation (reservoir for CO₂ storage) or a group of storage formations that has a sealing formation (a confining rock unit) above it.

References

Anderson, S. T. (2017). Risk, liability, and economic issues with long-term CO₂ storage—a review. Nat. Resour. Res. 26, 89–112. doi:10.1007/s11053-016-9303-6

CrossRef Full Text | Google Scholar

Anderson, S. T., and Jahediesfanjani, H. (2019). Estimating the pressure-limited dynamic capacity and costs of basin-scale CO₂ storage in a saline formation. Int. J. Greenh. Gas. Control 88, 156–167. doi:10.1016/j.ijggc.2019.05.031

CrossRef Full Text | Google Scholar

Anderson, S. T., and Jahediesfanjani, H. (2020). Estimating the net costs of brine production and disposal to expand pressure-limited dynamic capacity for basin-scale CO₂ storage in a saline formation. Int. J. Greenh. Gas. Control 102, 103161. doi:10.1016/j.ijggc.2020.103161

CrossRef Full Text | Google Scholar

Bandilla, K. W., Kraemer, S. R., and Birkholzer, J. T. (2012). Using semi-analytic solutions to approximate the area of potential impact for carbon dioxide injection. Int. J. Greenh. Gas. Control 8, 196–204. doi:10.1016/j.ijggc.2012.02.009

CrossRef Full Text | Google Scholar

Birkholzer, J. T., and Zhou, Q. (2009). Basin-scale hydrologic impacts of CO₂ storage: regulatory and capacity implications. Lawrence Berkeley National Laboratory (LBNL-1716E). Berkeley, CA. doi:10.2172/951197

CrossRef Full Text | Google Scholar

Bonferroni, C. E. (1936). “Teoria statistica delle classi e calcolo delle probabilità,” in Issue 8 of Pubblicazioni del R. Istituto superiore di scienze economiche e commerciali di Firenze, Florence R. Istituto superiore di scienze economiche e commerciali(Seeber). doi:10.4135/9781412961288.n455

CrossRef Full Text | Google Scholar

Brennan, S. T., Burruss, R. C., Merrill, M. D., Freeman, P. A., and Ruppert, L. F. (2010). A probabilistic assessment methodology for the evaluation of geologic carbon dioxide storage. U.S. Geological Survey Open-File Report 1127. Reston, Virginia. doi:10.3133/ofr20101127

CrossRef Full Text | Google Scholar

Casey, G. D. (1996). Hydrogeologic framework of the midwestern basins and arches region in parts of Indiana, Ohio, Michigan, and Illinois. U.S. Geological Survey Professional Paper 1423-B. Denver, Colorado. doi:10.3133/pp1423B

CrossRef Full Text | Google Scholar

Chen, Y. (2015). A new methodology of spatial cross-correlation analysis. PLoS ONE 10, e0126158. doi:10.1371/journal.pone.0126158

PubMed Abstract | CrossRef Full Text | Google Scholar

Duggan, Jr., J. E., Ogland-Hand, J. D., Anderson, S. T., and Middleton, R. S. (2024). “Managing basin-scale carbon sequestration: a tragedy of the commons approach,” in Proceedings of the 17th Greenhouse Gas Control Technologies Conference (GHGT-17), 20–24 October 2024, Calgary, Canada. doi:10.2139/ssrn.5030762

CrossRef Full Text | Google Scholar

Finley, R. J. (2014). An overview of the Illinois Basin – Decatur Project. Greenh. Gas. Sci. Technol. 4 (5), 571–579. doi:10.1002/ghg.1433

CrossRef Full Text | Google Scholar

Finsterle, S. (2007). iTOUGH2 user’s guide. Lawrence Berkeley National Laboratory (LBNL-40040). Berkeley, CA. Available online at: https://itough2.lbl.gov/itough2-users-guide/.

Google Scholar

Finsterle, S., Commer, M., Edmiston, J. K., Jung, Y., Kowalsky, M. B., Pau, G. S. H., et al. (2017). iTOUGH2: a multiphysics simulation-optimization framework for analyzing subsurface systems. Comput. Geosciences 108, 8–20. doi:10.1016/j.cageo.2016.09.005

CrossRef Full Text | Google Scholar

Gasda, S. E., Sandve, T. H., Tveit, S., Landa-Marbana, D., Pettersson, P., Krumscheid, S., et al. (2024). “Quantifying the impact of regional-scale pressure interference on commercial CO₂ storage targets for multiple licenses,” in Proceedings of the 17th greenhouse gas control technologies conference (GHGT-17), 20–24 October 2024, Calgary, Canada. doi:10.2139/ssrn.5053633

CrossRef Full Text | Google Scholar

Greenberg, S. E. (2021). An assessment of geologic carbon sequestration options in the Illinois basin: phase III. U.S. Department of Energy Technical Report (DOE-UIUC-42588). Champaign, Illinois: Illinois State Geological Survey. doi:10.2172/1806192

CrossRef Full Text | Google Scholar

Gresham, R. L., and Anderson, O. L. (2010). Legal and commercial models for pore-space access and use for geologic CO₂ sequestration. U. Pitt. L. Rev. 72, 701. doi:10.5195/lawreview.2011.170

CrossRef Full Text | Google Scholar

Gresham, R. L., McCoy, S. T., Apt, J., and Morgan, M. G. (2010). Implications of compensating property owners for geologic sequestration of CO₂. Environ. Sci. Technol. 44, 2897–2903. doi:10.1021/es902948u

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoholick, J. D., Metarko, T., and Potter, P. E. (1984). Regional variations of porosity and cement: St. Peter and Mount Simon sandstones in Illinois basin. Am. Assoc. Pet. Geol. Bull. 68, 753–764. doi:10.1306/ad461381-16f7-11d7-8645000102c1865d

CrossRef Full Text | Google Scholar

Intergovernmental Panel On Climate Change (2023). Climate Change 2022—mitigation of climate Change: working group III contribution to the sixth assessment report of the intergovernmental panel on climate Change. 1st ed. Cambridge, United Kingdom: Cambridge University Press. doi:10.1017/9781009157926

CrossRef Full Text | Google Scholar

International Energy Agency (2021). Net zero by 2050. Paris, France: International Energy Agency. Available online at: https://www.iea.org/reports/net-zero-by-2050.

Google Scholar

Jahediesfanjani, H., Anderson, S. T., and Warwick, P. D. (2019). Improving pressure-limited CO₂ storage capacity in saline formations by means of brine extraction. Int. J. Greenh. Gas. Control 88, 299–310. doi:10.1016/j.ijggc.2019.06.009

CrossRef Full Text | Google Scholar

Kendall, M. G. (1938). A new measure of rank correlation. Biometrika 30, 81–93. doi:10.1093/biomet/30.1-2.81

CrossRef Full Text | Google Scholar

Leetaru, H. E., Brown, A. L., Lee, D. W., Senel, O., and Coueslan, M. L. (2012). CO₂ injectivity, storage capacity, plume size, and reservoir and seal integrity of the Ordovician St. Peter Sandstone and the Cambrian Potosi formation in the Illinois basin. U.S. Department of Energy Topical Report (DOE/FE0002068-1). Champaign, IL: Illinois State Geological Survey. doi:10.2172/1064414

CrossRef Full Text | Google Scholar

Leetaru, H. E., Smith, V., Adushita, Y., and Freiburg, J. T. (2014). An integrated approach to evaluating the suitability of the Potosi Dolomite as a carbon sequestration target. Interpretation 2, SF125–SF133. doi:10.1190/INT-2013-0185.1

CrossRef Full Text | Google Scholar

Medina, C. R., Rupp, J. A., and Barnes, D. A. (2011). Effects of reduction in porosity and permeability with depth on storage capacity and injectivity in deep saline aquifers: a case study from the Mount Simon Sandstone aquifer. Int. J. Greenh. Gas. Control 5, 146–156. doi:10.1016/j.ijggc.2010.03.001

CrossRef Full Text | Google Scholar

Medina, C. R., Mastalerz, M., and Rupp, J. A. (2018). Pore system characterization of Cambrian-Ordovician carbonates using a new mercury porosimetry-based petrofacies classification system: An overview of the Illinois reservoirs. Greenh. Gas. Sci. Technol. 8 (5), 932–953. doi:10.1002/ghg.1806

CrossRef Full Text | Google Scholar

Miller, R. G. J. (1991). Simultaneous statistical inference. 2nd ed. New York, NY: Springer.

Google Scholar

Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522. doi:10.1029/WR012i003p00513

CrossRef Full Text | Google Scholar

National Petroleum Council (2019). Meeting the Dual Challenge—A Roadmap to 'At-Scale' Deployment of Carbon Capture, Use, and Storage. Washington, D.C., USA. Available online at: https://dualchallenge.npc.org/.

Google Scholar

Olea, R., and Luppens, J. (2014). Modeling uncertainty in coal resource assessments, with an application to a central area of the Gillette coal field. U.S. Geological Survey Scientific Investigations Report 2014–5196. Reston, Virginia. doi:10.3133/sir20145196

CrossRef Full Text | Google Scholar

Plampin, M. R., and Merrill, M. D. (2024). Hypothetical CO₂ leakage into and hydrological plume management within an underground source of drinking water at a proposed CO₂ storage facility, Kemper County, Mississippi, USA. Environ. Earth Sci. 84 (18), 18. doi:10.1007/s12665-024-11973-9

CrossRef Full Text | Google Scholar

Plampin, M. R., Anderson, S. T., Finsterle, S., and Cahan, S. M. (2023). Dynamic estimates of geologic CO₂ storage resources in the Illinois Basin constrained by reinjectivity of brine extracted for pressure management. Greenh. Gas. Sci. Technol. 13, 31–47. doi:10.1002/ghg.2189

CrossRef Full Text | Google Scholar

Plampin, M. R., Morrissey, E. A., Finsterle, S., Anderson, S. T., and Wiens, A. M. (2025). Multiphase simulations of CO₂ injectivity with and without brine extraction constrained by brine reinjectivity to optimize CO₂ storage in the Illinois Basin. U.S. Geological Survey Data Release. doi:10.5066/P14PBU3J

CrossRef Full Text | Google Scholar

Pruess, K. (2005). ECO2N: a TOUGH2 fluid property module for mixtures of water, NaCl and CO₂. Lawrence Berkeley National Laboratory (LBNL-57952). Berkeley, CA. Available online at: https://tough.lbl.gov/assets/docs/TOUGH2_ECO2N_Users_Guide.pdf.

Google Scholar

R Core Team (2024). R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/.

CrossRef Full Text | Google Scholar

Raza, A., Rezaee, R., Gholami, R., Bing, C. H., Nagarajan, R., and Hamid, M. A. (2016). A screening criterion for selection of suitable CO₂ storage sites. J. Nat. Gas Sci. Eng. 28, 317–327. doi:10.1016/j.jngse.2015.11.053

CrossRef Full Text | Google Scholar

Schmoker, J. W., and Halley, R. B. (1982). Carbonate porosity versus depth: a predictable relation for South Florida. Am. Assoc. Pet. Geol. Bull. 66, 2561–6570. doi:10.1306/03b5ac73-16d1-11d7-8645000102c1865d

CrossRef Full Text | Google Scholar

Swezey, C. S. (2009). Regional stratigraphy and Petroleum systems of the Illinois basin. Geological Survey Scientific Investigations Map 3068. Denver, CO. doi:10.3133/sim3068

CrossRef Full Text | Google Scholar

U.S. Department of Energy (2015). Carbon Storage Atlas, 5th edition. National Energy Technology Laboratory (DOE/NETL-2015/1709). Pittsburgh, PA. Available online at: https://www.netl.doe.gov/sites/default/files/2018-10/ATLAS-V-2015.pdf.

Google Scholar

U.S. Environmental Protection Agency (2014). “Part 144 - underground injection control program,” in Title 40 - protection of the environment (Chapter I - U.S. Environmental Protection Agency). Washington, DC, United States. Available online at: https://www.ecfr.gov/current/title-40/chapter-I/subchapter-D/part-144.

Google Scholar

U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team (2013). National assessment of geologic carbon dioxide storage resources—Data (ver. 1.1, September 2013, supersedes ver. 1.0 released June 26, 2013). U.S. Geological Suvey Data Series 774. Reston, VA. doi:10.3133/ds774

CrossRef Full Text | Google Scholar

van Genuchten, M.Th (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. doi:10.2136/sssaj1980.03615995004400050002x

CrossRef Full Text | Google Scholar

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L., François, R., et al. (2019). Welcome to the tidyverse. JOSS 4, 1686. doi:10.21105/joss.01686

CrossRef Full Text | Google Scholar

Will, R., Smith, V. L., and Leetaru, H. E. (2014). Utilization of the St. Peter Sandstone in the Illinois Basin for CO₂ sequestration. U.S. Department of Energy Topical Report (DOE/FE0002068-7). Champaign, IL: Illinois State Geological Survey. Available online at: https://www.osti.gov/biblio/1202234.

Google Scholar

Zoback, M. D., and Gorelick, S. M. (2012). Earthquake triggering and large-scale geologic storage of carbon dioxide. Proc. Natl. Acad. Sci. U.S.A. 109, 10164–10168. doi:10.1073/pnas.1202473109

PubMed Abstract | CrossRef Full Text | Google Scholar