Probabilistic evaluation of shallow groundwater resources at a hypothetical carbon sequestration site.

Carbon sequestration in geologic reservoirs is an important approach for mitigating greenhouse gases emissions to the atmosphere. This study first develops an integrated Monte Carlo method for simulating CO2 and brine leakage from carbon sequestration and subsequent geochemical interactions in shallow aquifers. Then, we estimate probability distributions of five risk proxies related to the likelihood and volume of changes in pH, total dissolved solids, and trace concentrations of lead, arsenic, and cadmium for two possible consequence thresholds. The results indicate that shallow groundwater resources may degrade locally around leakage points by reduced pH and increased total dissolved solids (TDS). The volumes of pH and TDS plumes are most sensitive to aquifer porosity, permeability, and CO2 and brine leakage rates. The estimated plume size of pH change is the largest, while that of cadmium is the smallest among the risk proxies. Plume volume distributions of arsenic and lead are similar to those of TDS. The scientific results from this study provide substantial insight for understanding risks of deep fluids leaking into shallow aquifers, determining the area of review, and designing monitoring networks at carbon sequestration sites.

Geologic carbon sequestration (GCS) in deep reservoirs is considered to be a viable long-term solution for mitigating greenhouse gases in the atmosphere1. But there are substantial regulatory and technical challenges that need to be overcome before GCS can be widely deployed. A major regulatory concern that has garnered much research attention is the potential leakage of CO2 and brine from deep storage reservoirs to shallow groundwater resources through preferential pathways, such as faults, wells, and local high-permeability zones in caprocks. The leaked CO2 and brine may cause a decrease in groundwater pH, an increase in total dissolved solids (TDS), and a potential mobilization of trace metals from aquifer materials into groundwater234567891011121314151617181920. While the potential risks to shallow aquifers at GCS operations are known, no defensible methodology for quantifying the risks has been reported to date34.

As part of the National Risk Assessment Partnership (NRAP) project21, we are developing approaches that can be used to quantify long-term risks at GCS sites. Current efforts are focused on analysis of risk proxies (i.e., impacts and their likelihoods) that will ultimately be combined with consequence analyses to create rigorous risk assessments. NRAP utilizes an Integrated Assessment Model (IAM) to predict long-term performance of GCS site. The IAM is developed using a system modeling approach that captures the physical and geochemical processes that occur within and among various sub-components of the sequestration sites, e.g., from sequestration reservoir to leaky well, then to underground sources of drinking water (USDW) and/or the atmosphere. As an integrated model, the IAM simulates the coupling processes from the point of CO2 injection to the interactions among CO2, brine, groundwater, and aquifer materials in case of a leak. Coupling the complex sub-component models involves significant computational costs, since the inherent heterogeneities and uncertainties in geologic systems require modeling a large number of realizations to bracket and quantify uncertainties associated with likelihoods, impacts, and risks. To meet the requirements of modeling coupled processes and performing large number of realizations, NRAP is developing and using reduced order models which are computationally efficient versions of the full process models that capture the underlying physical and chemical processes at reduced computational costs21.

This study demonstrates the developed methodology by using the highly heterogeneous unconfined Edwards aquifer as an example (Texas, USA). This aquifer was used because there is considerable data available, and because carbonate aquifers represent an important class of drinking water sources in the USA. Furthermore, the Edwards aquifer overlies potential sequestration sites including oil and gas reservoirs23. The Edwards aquifer is one of the most productive Karst aquifers in the world. It serves the diverse agricultural, industrial, recreational, and domestic needs of almost two million users in south central Texas22. The total groundwater withdrawal from the Edwards aquifer in 2000 was 2.8 × 106 m3/d of which withdrawal for public supply was about 1.6 × 106 m3/d. The Edwards aquifer consists of an underground layer of porous, fractured, honeycombed, water-bearing limestone rock with highly heterogeneous structures. In the recharge zone, the aquifer is unconfined. In the transition and discharge zones, it is an artesian aquifer confined by a layer of very low permeability sediment.

We select an unconfined zone of the Edwards aquifer as an example in order to evaluate the impact of the leaked CO2 and brine on the USDW and the potential CO2 discharge rate into the atmosphere (Note that other NRAP groups use the confined High Plain aquifer as an example to cover different types of shallow aquifers in the USA21). The study area is located north of San Antonio, where the unconfined carbonate aquifer has a total thickness between 100 and 200 m, with a mean thickness of 150 m. The groundwater flows from northwest to southeast with a hydraulic gradient around 0.00087 and the permeability structure of the unconfined carbonate is highly heterogeneous24. By reviewing the previous study of the permeability data2425262728, we summarize the statistical and spatial correlation parameters such as mean, variance, and integral scale in Table 1. The vertical integral scale is assumed to be linealy correlated to the horizontal integral scale with a factor of 0.01. The existing permeability data from Lindgren et al.2425 are used as pilot-point data for geostatistical simulations. A pilot-point-based Gaussian simulation method is modified from GSLIB29 to a geotatistical model GEOST3031 for simulating permeability fields. PSUADE32 and GEOST are coupled to generate nearly 500 realizations of the process model with Latin Hypercube sampling and geostatistical modeling for an integrated Monte Carlo simulation of CO2 reactive transport in the unconfined carbonate aquifer.

Table 1

Uncertain parameters and the ranges

#	Parameters	Min	Max
P1	Permeability variance	0.01	0.8
P2	Integral scale (km)	1.0	4.0
P3	Anisotropy (kx/kz)	1.0	50
P4	Mean permeability (m2)	10−14	10−10
P5	Mean Porosity (positively correlated to permeability)	0.05	0.35
P6	Max. CO2 pressure (MPa)	10	40
P7	Max. CO2 saturation	0.1	1.0
P8	Wellbore permeability (m2)	10−16	10−12
P9	Cl concentration (mol/l)	0.5	5.4

For each realization the uncertain parameters are sampled first and then the heterogeneous permeability fields for the unconfined aquifer are simulated with the pilot-point-based Gaussian method. Figure 1 shows one realization of the generated permeability fields in three dimensions. The model y axis is set along the groundwater flow direction and the lateral flow rate is varied with the sampled permeability and porosity. The numerical model size is 5000 × 8000 × 150 m with 64315 elements. The possible CO2 and brine leakage point is set at the upstream with a coordinate of (2500, 7000, 0), where the grid is highly refined with minimum grid sizes Δx, Δy, Δz of 2, 4, and 1 m, respectively. Away from the leakage point, the numerical grids become coarse (Figure 1).

Figure 1

The simulated heterogeneous permeability (logm2) field in a 3-D view.

The maximum CO2 pressure and saturation in the sequestration reservoirs and the wellbore permeability in the caprocks listed in Table 1 are three parameters we used to estimate the CO2 and brine leakage rates from the deep reservoirs to the shallow aquifer with a reduced order model developed as part of NRAP28. The reduced order model computes CO2 and brine leakage source terms using CO2 pressure and saturation in sequestration reservoir and the wellbore permeability. One example of the CO2 and brine leakage functions is shown in Figure 2, in which the CO2 injection time and rate in the deep reservoirs were assumed to be 50 years and 5 million tons per year, respectively. When CO2 and brine leaks into the unconfined carbonate aquifer, there could be very complex CO2-induced geochemical reactions in the aquifer, including aqueous equilibrium reactions (or acid-base and complexation reactions), trace metal adsorption or desorption, and mineral dissolution or precipitation111533. After reviewing all those reactions discussed by Bacon33 and Zheng et al.16, for carbonate aquifers such as Edwards aquifer, we incorporated the reactions, which have the most impact on the pH plume development, into our reactive transport simulations:

Figure 2

The simulated CO2 and brine leakage rates in a leaky well from one realization.

CO2 gas dissolution in water to reduce pH,

Aqueous equilibrium reactions and the equilibrium constants at 25°C,

Calcite kinetic dissolution or precipitation,

Calcite is the most plentiful mineral in the carbonate aquifer and can significantly buffer the pH. The forward rate for the kinetically-controlled calcite dissolution or precipitation at 25 °C is 1.5 × 10−6 mol/m2/s11. The specific surface area for calcite is 0.1 cm2/g. Dolomite is another common mineral in carbonate aquifers. In this study area the volume fraction of dolomite is very small relatively to that of calcite. We ignored dolomite in our geochemical calculations. We assumed that the brine is mainly made up of NaCl and treat its concentration as an uncertain variable, ranging from 0.5 to 5.4 mol/L. The geochemical analysis conducted by Bacon33 and Last et al.34 indicates that As, Cd, and Pb concentrations in the brine are proportional to the NaCl concentration. The proportional ratios are 3.16 × 10−7, 3.16 × 10−8 and 1.0 × 10−5 on a molar basis for As, Cd, and Pb, respectively. We simulated the species Cl, Na, and trace metals (As, Cd, and Pb) from the leaked brine as conservative tracers by ignoring the sorption reactions of the trace metals. Since sorption processes generally reduce the concentrations of the trace metals, our approach represents the conservative or worst-case scenarios for the trace metals.

By examining the existing water quality data in the unconfined Edwards aquifer we establish baseline data sets and statistical protocols for determining statistically significant changes between background concentrations and predicted concentrations (after CO2 and brine leakage)3334. The background (or initial) concentrations of the major chemical components, as well as their EPA maximum contaminant levels (MCLs) and no-impact thresholds, are listed in Table 23334. Note that the no-impact thresholds are determined for pH, TDS, As, Cd, and Pb based on their standard deviations of the measured water quality data and are used to identify potential areas (or volumes) of contamination predicted by aquifer numerical models. This study quantitatively evaluates the following five volume-based risk proxies: 1) plume volume of pH less than the MCL or a no-impact threshold in the aquifer; 2) plume volumes of TDS larger than the MCL or a no-impact threshold; and 3) plume volumes of the three trace metals (e.g. As, Cd, and Pb) over their MCLs and no-impact thresholds. Additionally, we calculate the CO2 discharge rate into the atmosphere (It should be noted that currently there is no consequence threshold used for CO2 discharge rate to the atmosphere).

Table 2

Initial values, no-impact thresholds, and MCLs for each component

Component	Initial values	No-impact thresholds	MCLs	Units
Arsenic	0.31	0.55	10	μg/L
Cadmium	0	0.04	5	μg/L
Lead	0.06	0.15	15	μg/L
pH	6.9	6.6	6.5	standard
TDS	330	420	500	mg/L

Results

The integrated MC simulation has been developed in this study for assessing the plume volumes of pH, TDS, and trace metals (e.g. As, Cd, and Pb), and CO2 discharge rates into the atmosphere under a probability framework. In most of the MC realizations the shallow groundwater resources are degraded locally around leakage points by reduced pH, increased total dissolved solids, and trace metals mobilized from the injection reservoir. The simulated pH plumes are the largest of the risk proxies and extend from 100 to 1500 m along the flow direction. Perpendicularly to the flow direction pH plumes extend less than 500 m. Subsequent to the leakage free phase of CO2 gas is formed, which preferentially migrates upwards due to buoyancy and leads to formation of pH plumes to the water table. The amount of CO2 leaving the top of the aquifer ranges between 0–98 of the CO2 leaking into the aquifer.

The global sensitivity results indicate that the volumes of pH less than the MCL and the no-impact threshold are most sensitive to aquifer porosity and CO2 leakage rates, as well as the cumulative CO2 mass and aquifer mean permeability. The volumes of TDS larger than the MCL and the no-impact threshold are most sensitive to the cumulative brine mass, brine leakage rate, and aquifer porosity. The CO2 discharge rates (into the atmosphere) are positively correlated to the leakage rate from deep reservoirs.

For all of the output variables pH plume volumes (measured with MCL or no-impact threshold) have the largest mean and median while the trace metal Cd plume volumes are the smallest. When measured with no-impact thresholds, the plume volume distributions of As and Pb are similar to those of TDS. Since the pH plume sizes are the largest among the five output variables, we will take it as one of the major risk proxies and use it as the reference to determine the area of review and to design the monitoring network in the unconfined shallow aquifer.

Discussion

The numerical results from the multi-phase reactive transport simulator FEHM35 show that when CO2 leaks into the aquifer, there is a small fraction of CO2 gas (less 10%) dissolved into water, which leads immediately to a pH reduction at the leak points. The groundwater lateral flow transports the gaseous and dissolved liquid CO2 downstream to extend the pH plumes. At the same time, the reduced pH enhances calcite dissolution to consume the increasing H+, which buffers the pH. Because the leaked gaseous CO2 has a much smaller density than water in the aquifer, the gaseous CO2 also transports upwards by buoyancy, which leads to pH plumes to extend vertically until the water table (see Figures 3A, 3B and 3C). Finally the gaseous CO2 leaves the top of the unconfined aquifer into the atmosphere. In most of the MC realizations, the simulated pH < 6.5 plumes extend from 100 to 1500 m along the flow direction. Perpendicularly to the flow direction pH plumes extend less than 500 m (Figure 3D).

Figure 3

Simulated pH plumes in the aquifer at different times from a YZ plane view at 40 (A), 100 (B), and 200 years (C), and a XY plane view (D) at four different times.

Global sensitivity analysis

By using the Monte Carlo simulation results as the input for PSUADE32, we conduct a global sensitivity analysis with the Main Effect method. At this step the three parameters that define the CO2 and brine leakage functions (e.g. maximum CO2 pressure and saturation in the squestration reservoirs and the wellbore permeability in the caprock) are converted to time-dependent leakage rates and the accumulated CO2 mass computed from the leakage functions for each realization. The results plotted in Figure 4 show that the plume volumes measured with MCLs and no-impact thresholds have the similar sensitivity to the input parameters, but different risk proxies (e.g. pH and TDS volumes) are sensitive to different parameters. For example, the volumes of pH less than the MCL and no-impact threshold are most sensitive to aquifer porosity and CO2 leakage rates, as well as the cumulative CO2 mass and aquifer mean permeability (Figure 4A). Since the aquifer porosity is linearly correlated to the volumes of pH plumes, it is most sensitive to that parameter. The volumes of TDS larger than the MCL and no-impact threshold are most sensitive to the cumulative brine mass, brine leakage rate and aquifer porosity (Figure 4B). The sensitivities of the volumes of the tracer metals (As, Cd, and Pb) less than their MCLs and thresholds are similar to those of the volumes of TDS. The CO2 discharge rate (into the atmosphere) is positively correlated to the leakage rate from the deep reservoirs.

Figure 4

Global sensitivity of the volumes of pH (A) and TDS plume volumes (B) to the 10 input parameters.

Probability measures for the risk proxies

By using the postprocessing results of the original 460 (200 year) MC simulations, we conduct a statistical analysis of the output variables with their MCLs and no-impact thresholds. The 5th, 50th (median), and 95th percentiles of the outputs are computed at 20 different time steps. The statistical analysis results of the CO2 discharge rate to the atmosphere are shown in Figure 5.

Figure 5

The statistical analysis results for the 5 groundwater risk proxies and atmospheric leakage of CO2 with their 5th, 50th (median), and 95th percentiles.

Since the MCLs and no-impact thresholds for pH and TDS are very close to each other (Table 2), the computed pH and TDS plume volumes with thresholds are slightly higher than their volumes with MCLs (Figures 5A and 5B). The MCLs and no-impact thresholds for trace metals (As, Cd, and Pb) are quite different and their percentile curves are also quite different (Figures 5C, 5D, and 5E). With high standard MCLs the plume volumes of the three trace metals are 0 at the 5th percentiles. The median plume volume of Cd is also 0 when measuring with its MCL (Figure 5D). The trace metal plume volumes measured with no-impact thresholds are significantly larger than those measured with MCLs. Surprisingly, the computed 5th, 50th (median), and 95th percentiles of the trace metal plume volumes with no-impact thresholds are very similar to the TDS volumes. Note that all these plume volumes still keep increasing slowly even though the CO2 and brine leakage decreases after 50 years. The time-dependent trend of the of the CO2 discharge rate to the atmosphere is similar to the CO2 leakage rate to the aquifer, but the discharge rate has a larger tail at the late time, which means that when the leakage rate (into aquifer) goes down after 50 years the gaseous CO2 in the unconfined aquifer still keeps releasing into the atmosphere with a gradually reduced rate. The 5th percentile of the CO2 discharge rate is 0 (a value of 10−8 in the log scale to represent 0 in Figure 5F).

Figure 6 shows histograms of the plume volume distributions of the 5 risk proxies at 200 years. Note that the log volume of 0 represents no plume. Among all these 5 risk proxies pH plume volumes (both with MCLs and thresholds) have the largest mean and median while the Cd plume volumes are the smallest. For As and Pb, 33% of the 460 realizations do not exceed the MCL threshold and 80% of the realizations do not exceed the MCL for Cd since Cd concentration in the leaked brine source is about one order of magnitude lower than those of As and Pb. When measured with no-impact thresholds, the plume volume distributions of As and Pb are similar to those of TDS. The volumes of the risk proxies measured with thresholds represent the potential areas (or volumes) of contamination (but not over the MCLs) in the shallow aquifer. Since the pH plume sizes are the largest among the five risk proxies, they will be used as the reference for determining the area of review and to design the monitoring network for the unconfined shallow aquifer.

Figure 6

The histograms of the computed plume volumes for pH, TDS, As, Cd and Pb measured with their MCLs and no-impact thresholds at 200 years.

Implications

An integrated Monte Carlo (MC) simulation under a probability framework has been developed in this study for assessing risk proxies (plume volumes and probabilities) related to shallow groundwater resources in the unconfined carbonate aquifer. Risk proxies are compared to both MCL and no-impact consequence thresholds. Differences in risk proxies between the two example consequence thresholds (MCL and no-impact) are significant, especially for the trace metals. This implies that careful consideration of consequence threshold definitions must be made.

Even though our study is for a hypothetical site with CO2 and brine leakage, results of our study provide substantial insights into the impact of the potential CO2 and brine leakage on shallow groundwater quality for the range of parameters studied. Results presented here will be used to develop reduced order models for shallow aquifers which will be incorporated in IAMs and subsequently used in risk assessment calculations. When more sitespecific aquifer parameters are obtained from laboratory-scale to field-scale tests, calculations of uncertain parameters must be revisited to refine estimates of all risk proxies. Finally, since the pH plume sizes are the largest among the five risk proxies, they will be used as the reference for determining the area of review and to design effective monitoring network for shallow aquifers.

Methods

Previous studies in the Edwards aquifer provide plentiful information about the ranges and distributions of the aquifer heterogeneous parameters, such as the aquifer thickness, hydraulic gradient, conductivity, and porosity24252627. Based on the statistics of these parameters, we develop a comprehensive methodology for conducting an integrated MC simulation of water and CO2 flow and reactive transport in the unconfined aquifer. The integrated MC simulation sequentially uses the uncertainty quantification code PSUADE32, the geostatistical modeling code GEOST3031 (modified from GSLIB29) and the multi-phase reactive transport simulator FEHM35. PSUADE is used to sample all of the uncertain parameters based on the statistical analysis results of the field data at or near the study area, conduct global sensitivity analysis of the output variables in relation to the input parameters, and evaluate the statistical distributions of the pH, TDS, and trace metal plume volumes. GEOST is used to analyze the heterogeneous structures of the Edwards aquifer and to simulate/generate the heterogeneous fields for permeability and porosity based on the heterogeneity parameters sampled from PSUADE.

The multi-phase reactive transport simulator FEHM is used to model CO2 and water flow and reactive transport in the unconfined aquifer for each generated heterogeneous model for the unconfined carbonate aquifer. The Brooks and Corey3536 relative permeability is used for water/CO2 multiphase flow simulations and the related coefficients are adopted from a literature review35363738. The numerical simulations always start from steady-state flow simulations as the initial conditions and then simulate the CO2 multi-phase reactive transport in the aquifer for 200 years. In the transport model, it is assumed that the molecular diffusion coefficient is 10−9 m2/s and the longitudinal and transverse dispersivities are 800 m and 500 m, respectively. For each simulation a post-process is conducted to evaluate the four risk proxies, including the plume volumes of pH, TDS, and trace metals, and the CO2 discharge rate into the atmosphere at 20 different times, which means that each risk factor will have 20 time-dependent values. During the MC simulations we keep checking the variance and mean of the risk proxies. After 300 MC runs we find that their variance and mean are approaching constants. Then, we stop our MC simulations after completing 460 runs when the MC simulations are considered to be converged. By taking into account the 20 times as an additional dimension, we obtain the total number of outputs of 9200 (realization number (460) times the number of time steps (20)). These data are re-organized in a PSUADE input format for global sensitivity and risk analysis conducted under a probability framework to evaluate quantitatively the chosen risk proxies.

In order to determine the key flow and transport parameters driving CO2 transport behavior in the shallow aquifer, global sensitivity analysis techniques are used for investigating input-output sensitivities that are valid over the entire range of the parameter sampling. The Monte Carlo simulations provide 460 realizations of input transport parameter sets generated with Latin Hypercube sampling and geostatistical modeling. Each realization is propagated through the transport simulator to yield the output response functions represented by the output variables. Global sensitivity analysis entails the comparison of the output distributions to each of the input parameter distributions and identifies the most sensitive parameters for each output variable. The main effect method32 is used to quantify the uncertainty and sensitivity of the output variables to input parameters.

The main effect method is a variance-based analysis and it displays first-order Sobol' indices for the response surface built from the Monte Carlo simulations32. The essence of this analysis is the statistical measure called variance of conditional expectation. The variance-based analysis uses the following equation: where, VCE measures the variability in the conditional expected values of output variable Y when the input parameter X takes on different values. s is the number of distinct values of each input parameter and r is the number of replications. N = sr is the sample size. The sensitivity of the output variable to the input parameters is quantified by Equation (1) and re-ranked with numbers from 0 to 100 to represent the importance of the input parameters32.

Author Contributions

Z.D. designed and carried out the integrated MC simulations, performed data analysis and wrote the draft manuscript. E.K., R.P. and D.B. designed the numerical and geochemical models, discussed the results, reviewed and revised the manuscript. H.V., P.S. and A.J. discussed the results, and revised the manuscript.