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.

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.

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.

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.

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.

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 (Φ_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.

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).

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).

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:Vinj,t=aϕbHcAdkeDfTg(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:Vinj,t *=1.53−Vinj,t25,000Vext,t−25,000×Vinj,tVinj,t−1.53(2)

In Equation 2, Vinj,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.

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).

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.

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.

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.

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.

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.

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).

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.