Human Health Risk Assessment of CO2 Leakage into Overlying

Article pubs.acs.org/est

Human Health Risk Assessment of CO2 Leakage into Overlying Aquifers Using a Stochastic, Geochemical Reactive Transport Approach Adam L. Atchley, Reed M. Maxwell,* and Alexis K. Navarre-Sitchler Hydrological Science & Engineering Program, Department of Geology & Geological Engineering, Colorado School of Mines, 1500 Illinois St, Golden, Colorado 80401, United States

ABSTRACT: Increased human health risk associated with groundwater contamination from potential carbon dioxide (CO2) leakage into a potable aquifer is predicted by conducting a joint uncertainty and variability (JUV) risk assessment. The approach presented here explicitly incorporates heterogeneous flow and geochemical reactive transport in an efficient manner and is used to evaluate how differences in representation of subsurface physical heterogeneity and geochemical reactions change the calculated risk for the same hypothetical aquifer scenario where a CO2 leak induces increased lead (Pb2+) concentrations through dissolution of galena (PbS). A nested Monte Carlo approach was used to take Pb2+ concentrations at a well from an ensemble of numerical reactive transport simulations (uncertainty) and sample within a population of potentially exposed individuals (variability) to calculate risk as a function of both uncertainty and variability. Pb2+ concentrations at the well were determined with numerical reactive transport simulation ensembles using a streamline technique in a heterogeneous 3D aquifer. Three ensembles with variances of log hydraulic conductivity (σ2lnK) of 1, 3.61, and 16 were simulated. Under the conditions simulated, calculated risk is shown to be a function of the strength of subsurface heterogeneity, σ2lnK and the choice between calculating Pb2+ concentrations in groundwater using equilibrium with galena and kinetic mineral reaction rates. Calculated risk increased with an increase in σ2lnK of 1 to 3.61, but decreased when σ2lnK was increased from 3.61 to 16 for all but the highest percentiles of uncertainty. Using a Pb2+ concentration in equilibrium with galena under CO2 leakage conditions (PCO2 = 30 bar) resulted in lower estimated risk than the simulations where Pb2+ concentrations were calculated using kinetic mass transfer reaction rates for galena dissolution and precipitation. This study highlights the importance of understanding both hydrologic and geochemical conditions when numerical simulations are used to perform quantitative risk calculations.

■

INTRODUCTION AND MOTIVATION

trations, which is in turn related to uncertainty in spreading and migration of contaminants in the subsurface.12,16−20 This uncertainty can be attributed to sparse subsurface data needed to characterize subsurface properties.21−23 To account for these uncertainties, stochastic methods, such as Monte Carlo

Metal release associated with CO2 leakage into groundwater from carbon capture utilization and storage (CCUS) projects can potentially increase human health risk.27,30,10 Quantitative assessment of this potential increased risk is best performed with numerical simulation techniques that account for uncertainty associated with representation of the subsurface physical heterogeneity and geochemical processes1−3 In many cases, uncertainty in human health risk assessment is directly related to uncertainty in groundwater contaminant concen© XXXX American Chemical Society

Received: January 21, 2013 Revised: April 12, 2013 Accepted: April 25, 2013

A

dx.doi.org/10.1021/es400316c | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

transport for three ensembles of 100 equally probable heterogeneous flow realizations with variance of log hydraulic conductivity (σ2lnK) of 1, 3.61, and 16 Risk from Pb2+ concentrations assuming equilibrium with galena is compared to a reactive transport model that simulates reaction rates governing Pb2+ release into groundwater.

approaches, are used to calculate a range of expected contaminant concentrations from ensembles of equally probable subsurface configurations, for example, see refs 17,24−26. This range of concentrations is then sampled, taking into account the probability of a given concentration for calculation of human health risk using variable exposure parameters.10,11,16,19 A joint uncertainty and variability analysis (JUV) is one such approach that quantifies risk through a nested Monte Carlo ensemble where lack of knowledge within the groundwater system is quantified through uncertain parameter distributions in one Monte Carlo loop and results from those groundwater simulations are fed into a second Monte Carlo loop that samples over variability in individual response to contamination exposure within a population.4−12 The power of JUV analysis as a decision making tool is the ability to calculate human health risk as a function of uncertainty in subsurface transport and variability in exposure time, metabolized dose, and low-dose effects within a population.14,15 These established methods have proven robust when considering complex hydrologic uncertainty4−9,13−15,17−19 and relatively simple geochemical processes.10,11,20 However the impact of multicomponent, kinetic geochemical processes such as mineral dissolution on calculated risk remains relatively unexplored, partly due to the computational expense of ensembles of 3D, multicomponent, hydrogeochemical simulations with kinetic mass transfer reactions.44,73 Here we use a streamline technique and employ parallel processing to increase the efficiency of the reactive transport simulations, such that ensembles of multicomponent geochemical simulations with kinetic mass transfer rates are not a limiting factor in risk assessment. Thus, allowing for the evaluation of calculated risk where contaminant concentrations are determined through kinetic mass transfer reactions along a flow path. CO2 leakage from an underground storage reservoir reduces aquifer pH, inducing geochemical reactions in shallow aquifers, and potentially mobilizes EPA regulated contaminants like metals.27−35,74 For example, several studies have indicated that lead concentrations will increase due to galena dissolution in reducing aquifers.27,30 Unlike many other risk assessment scenarios where the contaminant is traced to a single source, in CO2 leakage scenarios the source of metals is the rock matrix itself. Therefore, metal concentration is a function of aquifer mineralogy and the rate at which metals are released. Mineral dissolution rates are a function of mineral saturation state. For metal sulfides like galena, aqueous coupling of sulfur with protons (to form HS−) creates a dependence of mineral saturation on pH. Thus, lowering aquifer pH through dissolution of CO2 into the groundwater and subsequent dissociation of carbonic acid will lead to undersaturation of the groundwater with respect to galena, driving galena dissolution.36,37 Fluid residence times impact the saturation state of the groundwater with respect to reacting minerals39,43,44 resulting in distributions of metals concentrations and mineral saturation states in the plume that reflect distributions of fluid residence times.38−42 Given the influence of uncertainty in concentrations resulting from heterogeneous subsurface flow and reactions in previous risk assessment work,8,10,14−17,19 we hypothesized that the method of geochemical modeling and representation of macrodispersion will impact estimated risk related to increased Pb2+ resulting from a CO2 leak into a reducing aquifer. This study calculates human health risk using previous work,44 which evaluated kinetic mass transfer reactive

■

MATERIALS AND METHODS This paper uses a stochastic, geochemical reactive transport model in a JUV human health risk assessment: a key advancement in risk assessment. Changes in pH and reactive mineral surface areas are accounted for when calculating dissolution and precipitation rates. Streamlines are used to model reactive transport along flow paths between a chemical perturbation and a point of interest such as a pumping well used for drinking water. The streamline approach follows that of Atchley et al,44 where flow paths are mapped through a 3-D flow-field and transformed into 1-D streamlines by a time-offlight formulation(see also refs 48−50). The streamline method is a well-developed approach for simulating contaminant transport. (e.g., refs 23,43,51−59) However, relatively few studies have employed streamlines for geochemical reactive transport studies.43,44,60 Stochastic streamline approaches to characterize aquifer transport have been developed,51,52,56 but to our knowledge this paper presents the only use of streamlines to simulate kinetic geochemistry within a human health risk framework that uses a Monte Carlo approach to explicitly account for subsurface uncertainty.

■

REACTIVE TRANSPORT CrunchFlow2007, a multicomponent reactive transport code61−63 was used to solve kinetic reactive transport along all streamlines. Dissolution and precipitation of mineral species is based on transition state theory36,37 where the rate (rm) of dissolution/precipitation of mineral m is related to the saturation state of the solution as a function of the ion activity product, Qm, of mineral m. This can be expressed by ⎛ Qm ⎞ ⎟ rm = A mkm⎜⎜1 − K m(eq) ⎟⎠ ⎝

(1) 2

where km is the intrinsic rate constant [mol/m /s], Am is the bulk mineral surface area [m2/m3], a dynamic mineral property that changes with dissolution and precipitation reactions.63 Precipitation occurs when rm is positive and dissolution occurs when rm is negative. Saturation Index (SI), a measure of equilibrium between the mineral and aqueous species of the mineral, can then be calculated by

⎛ Q ⎞ m ⎟ SI = log⎜⎜ ⎟ ⎝ K m(eq) ⎠

(2)

where Km(eq) is the corresponding mineral equilibrium constant. A SI of zero where Qm is equal to Km(eq) indicates equilibrium. In oversaturated conditions, SI is greater than zero and mineral precipitation may occur. In undersaturated conditions SI is less than zero and mineral dissolution occurs, releasing ions into solution.

■

HUMAN HEALTH RISK CALCULATIONS Risk calculations are based on the U.S. Risk Assessment Guidance for Superfund sites (RAGS), Volumes I and III2,61 a B



Article

Table 1. Human Health Risk Parametersa parameter

symbol/units

values

distribution

averaging time cancer potency factor, dermal cancer potencey factor, ingestion dermal permeability coefficient in water exposure duration exposure frequency fraction of skin in contact with water ingestion rate per unit body weight shower exposure duration skin surface are per unit body weight unite conversion factor

AT[d] CPFdermal[kg d/mg] CPFingestion [kg d/mg] Kp [m/h] ED [yr] EF[d/yr] fskin[−] IR/BW [L kg/d] Edshower [h/d] SA/BW [m2/kg] CF [L/m3]

25 550 0.0085 0.0085 0.00001−0.01 70 350 0.4−0.9 0.033, 0.013 0.13, 0.009 0.027, 0.0025 0.001

constant constant constant uniform constant constant uniform lognormal lognormal lognormal constant

source EPA 2001 TAC TAC EPA 2004 EPA 2001 EPA 2001 Mckone and Mckone and Mckone and Mckone and

Bogen Bogen Bogen Bogen

1991 1992 1993 1994

a

EPA 2001 US EPA. Risk assessment guidance for superfund. Proscess for conducting probablisitic Risk Assessment (Part A), vol. III. Washington DC; 2001. EPA 2004 US EPA. Risk assessment guidance for superfund. Human health evaluation manual (Part E, supplemental guidance for dermal risk assessment), vol. I. Washington, DC; 2004. Mckone and Bogen 1991 Mckone TE, Bogen KT. Predicting the uncertainties in risk assessmenta California groundwater case-study. Environ. Sci. Technol. 1991, 25(10), 1674−1681. TAC Toxic Air Contaminant document, office of Environmental Health Hazard Assessment (OEHHA).

Figure 1. Joint Uncertainty and Variability Risk Schematic depicting multiple realizations of uncertain subsurface configurations, variable human responses in a nested Monte Carlo framework. Together, uncertain subsurface realizations (Monte Carlo Step 1) and representative human responses (Monte Carlo Step 2) produce a risk curve over a range of uncertainty for each fractile. Collectively, all fractiles of variability produce a three-dimensional risk surface representing both uncertainty and variability. The direction of groundwater flow, source-zone and pumping well location is mapped for each realization depicted. The dimension of the source-zone is 5 m, 10 m, and 13 m in the X, Y, and Z directions. C



Article

Table 2. Geochemical Conditions (Adapted from Atchley et al. (2012).43a mineral content quartz calcite galena H2S(aq) porosity

reaction stoichiometry CaCO3 + H+ ↔ Ca2+ + HCO3− PbS + H+ ↔ Pb2+ + HS− H2S ↔ HS− + H+

percent volume

surface area Am [m2/m3]

rate constant km [mol/m2/s] (log)

64 3 0.1

0.508 0.00456 0.073

−15 −11.5 −11

log Ksp at 25 °C 1.8487 −14.8544 −7.9759

33 condition aquifer aqueous complexes

CaCO3(aq) CaCl+ CaCl2(aq) CaHCO3+ CaOH+ CO2(aq) H2SiO42− H2S(aq) HSiO3−

species +

PbCl PbCl2(aq) PbOH+ Pb(OH)2(aq) Pb(OH)3− Pb(CO3)22− S2− OH− CO3−

condition leak −1

conc. [mol L ] −04

1.00 × 10 4.00 × 10−12 1.00 × 10−10 1.00 × 10−03 charge 8.20 × 10−03 1 × 10−70 [bar] 0.01 [bar] 7

SiO2(aq) Pb2++ HS− Na+ Cl− Ca2++ O2(aq) → O2(g) HCO3− → CO2(g) pH

species

conc. [mol L−1]

SiO2(aq) Pb2+ HS− Na+ Cl− Ca2+ O2(aq) → O2(g) HCO3 → CO2(g) pH

1.00 × 10−04 4.00 × 10−12 1.00 × 10−10 1.00 × 10−03 charge 8.20 × 10−03 1 × 10−70 [bar] 30 [bar] 4.5

a

First layer of Table 2 describes the mineral compositions and associated mineral stoichiometry, percent volume, surface areas, rate constants, and equilibrium constants that determine reaction rates. The second layer of Table 2 lists Aquifer and leak conditions, as well as the aqueous species that describe the first and secondary aqueous complexes.

with water per unit body weight [m2/kg], EF is the standard exposure frequency [d/yr]. Specific to the ADD dermal pathway is Kp,m the dermal permeability coefficient of the contaminant in water [−], fskin is the fraction of skin in contact with water [−], EDshower is the exposure duration of the shower [h/d], and CF is a unit conversion factor (1 × 10−3 L/m3). The inhalation pathway is a result of contaminant volatilization and with a few exceptions is not considered in studies concerning risk due to metals exposure10 and therefore is neglected in this study. Risk from the two considered pathways can be summed to calculate total risk by

thorough description of the Monte Carlo application to risk assessment is provided by Bogen and Spear,4 Maxwell and Kastenberg,8 and Siirila et al.10 JUV can be used for carcinogenic and noncarcogenic risk (e.g., refs 8,10), here we focus on carcinogenic risk. Increased individual lifetime cancer risk is a result of exposure to a carcinogen in up to three pathways: ingestion, dermal absorption, and inhalation. The risk from each exposure pathway is found by risk m = 1 − eCPFm × ADDm × f * m

(3)

where CPFm is the cancer potency factor [kg d/mg] for contaminant m, ADDm is the average daily dose [mg/kg d] of the contaminant, and f *m is the metabolized fraction of the contaminant. CPF and f * are parameters specific to each contaminant. ADDm is a function of the maximum contaminant concentration (Cm) and exposure duration (ED) [yr] recorded over the averaging time (AT) [yr], and a distribution of parameters that represent population variability. In the streamline method, Cm is calculated as a sum of the product of peak contaminant concentration ([Cm]istl) and hydrological flux (Qistl) of each streamline over the total pumping rate (Qtot) to the receptor. Shown here as [Cm]tot =

∑

m=n m=1

(4)

The ingestion and dermal absorptions ADD pathways are calculated from ⎡ IR ⎤ ED × EF ADDm ,ingestion = Cm⎢ ⎣ BW ⎦⎥ AT

m=1

(6)

The human health risk parameters for lead exposure used in this study are listed in Table 1. A log-normal variation for human exposure was used to capture population variance for ingestion rate, shower duration, and skin surface area. Constant health risk parameters such as CPF and Kp are specific to lead. Here, ensembles of 100 equally probable subsurface fields are simulated to capture uncertainty of spatial hydraulic conductivity heterogeneity in the first Monte Carlos step (Figure 1). Variable population parameters are then evaluated in the second Monte Carlo step with each uncertain representation of the subsurface, following the method presented in previous JUV work.4−8,10 From the JUV analysis a surface of risk as a function of both uncertainty and variability is produced. Cuts from this surface are plotted as cumulative distribution functions (CDFs) in either uncertainty or variability. Variability within the human population over a range of sensitivity is distinguished as fractiles of variability. Each fractile can be thought of as an individual within a population of varying susceptibilities controlled by variable physiological parameters. Low fractiles of variability characterize a subset of the population with low susceptibility. Conversely, high fractiles of variability characterize high susceptibility in individuals that are likely to experience adverse health effects, whereas the 50th fractile of variability represents the average risk susceptibility within the human population. In this way the JUV can assess

i [Cm]istl Q stl

Q tot

m=n

∑ risk = ∑ risk i= ingestion + ∑ risk i= dermal

(5a)

and ⎡ SA ⎤ ED × EF K p,mfskin EDshower CF ADDm ,dermal = Cm⎢ ⎣ BW ⎥⎦ AT (5b)

where IR/BW is the ingestion rate of water per unit body weight [L kg/d], and SA/BW is the skin surface area in contact D



Article

uncertainty in subsurface reactive transport). The increased lifetime cancer risk, associated with elevated lead exposure calculated from eq 6 and averaged over each ensemble, is plotted in Figure 2 for the 50th and 95th fractile of variability,

uncertainty for a spectrum of susceptibility within a human population (Figure 1).

■

DOMAIN CONFIGURATION A hypothetical CO2 leak was modeled in an idealized reducing aquifer 900 m up-gradient from a drinking water well pumping 300 m3/day. The CO2 leak with a PCO2 of 30 bar was simulated at the base of the aquifer creating a source-zone with a volume of 135 m3. We assume that all CO2 gas completely dissolved into the groundwater resulting in a moderate leak rate of approximately 4.84e 5 kg CO 2 /yr. The source-zone coordinates are X: 100−105 m, Y: 145−155 m, and Z: 2−15 m. Pumping well coordinates are X: 1000 m, Y: 150 m, and Z: 25−55 m. The bottom casing of the pumping well is 10 m vertically offset from the top of the CO2 source-zone. Therefore, hydrological heterogeneity is needed for a hydrological connection between the CO2 source-zone and the pumping well to occur. For detailed domain setup and boundary conditions see Atchley et al., (2013).44 One hundred realizations with unique spatial distributions of hydraulic conductivity resulting in multiple, equally probable transport scenarios were simulated for each ensemble. Streamlines were initiated from the well-screen area. For each simulation, 3131 streamlines were launched evenly around the pumping well in order to characterize the well capture area. These streamline paths were then traced backward through the domain using particle tracking methods.65,66 Only streamlines that crossed the CO2 source-zone were selected for transport simulation. The streamline definition used here dictates a constant flux in time for each streamline and mixing between streamlines does not occur. Unique streamline fluxes (QSTL) were assigned based on the permeability of the grid cell surrounding the pumping well in which the streamline was initiated. The average streamline flux was 0.034% of the total well flux, but varied between streamlines. The geochemical conditions from the previous work Atchley44 are used to simulate kinetic conditions within an aquifer comprised of 64% quartz by volume, 3% calcite, acting as a pH buffer and 0.1% galena providing a source of lead from mineral dissolution. Kinetic parameters and chemical composition are summarized in Table 2. The geochemical scenario represents simplified mineralogy but is illustrative of more complex systems and is consistent with recent modeling approaches for leakage problems considering potable water aquifers.10,67 The partial pressure of CO2 gas (PCO2) was initially set to 0.01 bar, pH was set to 7 along the entire length of all streamlines at t = 0 and simulated at equilibrium to provide background lead concentrations. A CO2 leak condition (Table 2) with a pH of 4.5 resulting from an increase in PCO2 to 30 bar was applied at the source-zone. The leak condition was applied to all streamlines at the point of exit from the CO2 source-zone. Dissolution/precipitation reactions along streamlines were modeled for the duration of the simulation.

Figure 2. Risk cuts from the 50th (plot a) and 95th fractile (plot b) of variability are plotted for each ensemble. Inset shows the risk cuts for the top 90−100 percentile of uncertainty.

representing an average and a highly exposed individual, respectively. The action level for remediation due to increased cancer risk described by the EPA2 is 1× 10−6. For the hydrogeochemical conditions simulated here, a small increase in cancer risk below the threshold required for remediation is found for all ensembles and as expected risk increased for higher fractiles of variability (Figure 2b). Figure 2 shows that risk increases as σ2lnK increases from 1 to 3.61 and then decreases as σ2lnK increases from 3.61 to 16 for all but the 99th and 100th percentiles of uncertainty. Above the 99th percentile, risk decreases as σ2lnK increases for all cases (Figure 2, inset). The shape of the risk curve for each ensemble is the same for the 50th and 95th fractiles of variability due to the uncertain subsurface characterizations, but differs for each σ2lnK ensemble. The unique risk curve shape and uncertainty for each σ2lnK tested is a reflection of combined hydro-geochemical processes at work, described in the following. Subsurface hydrological and geochemical processes control lead concentrations at the pumping well (Figure 3), which in turn contribute to the estimated human health risk (Figure 2). Here we see that σ2lnK affects the time needed for reactive solutes to reach the pumping well and therefore the geochemical equilibrium at the well, noted as the Damköhler number (Da = tR/tA) (Figure 3d). The Damköhler number is a nondimensional ratio of reaction time (tR) over the time to advect (tA) reactive solute from the point where reactions begin to an established end point68,69 (Figure 3d). Equilibrium is met when the Damköhler number is equal to or greater than one. A Damköhler number less than one represents an advective system that does not reach equilibrium and should be represented by kinetic reactions rates (e.g., refs 68 and 69). The time to reach a new geochemical equilibrium is approximately 500 days (tR), plotted in Figure 4. The state of

■

RESULTS AND DISCUSSION Concentrations of the contaminant at the pumping well and population exposure parameters govern estimated risk (eqs 5a, 5b, and 6). For this study it is assumed that human population exposure parameters are fully characterized. Therefore, estimated risk is a function of both interindividual variability and uncertainty in lead concentrations (which reflect E



Article

Figure 3. Cumulative distribution Function (CDF) of Pb 2+ concentration for all realizations (a). CDF of CO2 chemically impacted groundwater flux to the pumping well for all realizations (b). Cumulative probability numbers on the y-axis for plots a, and b refer to the same realization. Average streamline Pb2+ concentration for each realization (c). Pb2+ concentration plotted in c is not flux weighted, and does not include hydrological controls of Pb2+ concentration. Plot d shows the average Damkölher number for each ensemble, the time to equilibrium (tA) = 500 days. Error bars on plots c and d show the standard deviation.

reaction progress toward chemical equilibrium then determines the average ensemble streamline metals concentration (Figure 3c). In our simulations we see that as σ2lnK increases, reactions are further from equilibrium resulting in an increase in streamline lead concentrations. However, the flux of chemically impacted water advected to the pumping well (Figure 3b) exhibit strong controls over the estimated concentration of lead at the pumping well (Figure 3a) and therefore influence risk estimates. A hydrological connection must be established between the CO2 source-zone and the pumping well before increased human health risk is observed (eq 4 and Figure 3a). 70% of the realizations in the σ2lnK =1 ensemble had no streamlines crossing the CO2 source-zone and therefore resulted in no risk for those realizations (Figure 2 and 3). More than 90% of the realizations in the σ2lnK = 3.61, and 16 ensembles had streamlines connecting the CO2 source-zone to the pumping well and similarly saw an increase in human health risk (Figures 2 and 3) over background levels. The flux from streamlines that connect the source to the well also influences the predicted concentration, which then contributes to human health risk estimates. Of the three ensembles, the σ2lnK of 3.61 ensemble had the greatest flux for 99 realizations, and therefore the largest mean risk, but the σ2lnK = 1 ensemble had the highest flux for one realization so the 100th percentile risk for the σ2lnK = 1 ensemble is the greatest (Figure 3a and b). Furthermore, the σ2lnK = 16 ensemble had the highest average streamline lead concentration for all realizations with a source-to-well connection (Figure 3c), but a lower risk compared to the σ2lnK = 3.61 ensemble, because the source-to-well flux was consistently higher in the σ2lnK = 3.61 ensemble (Figure 3b). Given the importance of the hydrological flux from the CO2 source-zone to the pumping well when estimating pumping well contamination, proper characterization of properties that influence flux will aid risk evaluation. Here the σ2lnK is shown to

Figure 4. Saturation profiles along a streamline for galena and calcite are plotted against residence time. Note that Calcite and Galena are plotted on different axis and scaled to show changes in SI over time. Reaction rate profiles for calcite and galena are also plotted against time. pH and Lead concentrations are plotted against residence time to show changes to solute makeup corresponding to changes in equilibrium conditions. tR is shown in both plots and represents the time to equilibrium at 500 days. Profiles from both plots were taken after 3000 days of simulation time.

influence the source-to-well connectivity, but other factors will also play a role, such as aquifer stratification.10,15,17 The kinetic rate of reactions simulated resulted in dynamic metal concentrations along streamlines until a new chemical equilibrium is met. The reaction path after the intrusion of CO2 and associated drop in pH first triggered galena dissolution, which resulted in a sharp initial increase in lead concentration (Figure 4). However, calcite simultaneously dissolved and buffered the pH, resulting in a gradual increase in pH along the flow path. As the pH increased, groundwater became oversaturated with respect to galena and induced lead precipitation until a new chemical equilibrium was reached at approximately 500 days, when mineral reaction rates approached zero. This new chemical equilibrium was dictated by aquifer mineralogy and the initial conditions of the CO2 leak. The reaction rates of each reactive chemical species and the reactive surface area of the mineralogy determine the time until equilibrium is met37,70 along with the lead concentration at any point along the reaction path. While galena precipitation may not occur in shallow groundwater or surface water systems71,72 lead precipitation could occur in other forms. Therefore, the simple representation of lead precipitation modeled here can illustrate the kinetic response from more F



■

Article

IMPLICATIONS An efficient approach to estimate human health risk is presented that advances the capabilities of probabilistic risk assessments12 by simulating sophisticated hydro-geochemical processes. By employing a streamline hydro-geochemical technique with ensembles of 100 realizations for σ2lnK equal to 1, 3.61, and 16 we show that hydrological and geochemical conditions affect estimated health risk over a range of uncertain subsurface conditions. Human health risk increases as σ2lnK increases from 1 to 3.61 and decreases as σ2lnK increases from 3.61 to 16 for all but the highest percentiles of uncertainty. Above the 96th percentile of uncertainty health risk decreases as σ2lnK increases. We show that the unique risk curve shape for each ensemble is a result of the uncertain hydrological flux represented and geochemical state of reactive solute for each realization, suggesting that proper aquifer characterization is required to assess health risks from CCUS implementation. Lastly, a kinetic representation of geochemical reactions resulted in protective human health risk estimates compared to equilibrium simulations. Here the simple geochemical equilibrium assumption was shown not to be protective, illustrating the importance of considering kinetic mass transfer and transport with uncertainty in human health risk estimates. The estimated risk to human health from CCUS leaks into potable groundwater was found to be low for this hypothetical scenario. However, the purpose of these results are to showcase efficient reactive transport simulations within a human health risk framework, specific conditions to any site in question should be considered for actual health risk assessments. Moreover, recent studies have found increases in more than one carcinogenic metal.31,32 Therefore, a full risk assessment should consider the health risk from a number of host-rock metals rather than just lead. Among the many factors that affect heath risk due to CCUS implementation, coupled hydrogeochemical mechanisms were shown to be important.

complex systems. Additionally, other mineralogal compositions could affect the estimated metals concentration. For example, the lack of a calcite buffer or the presence of other mineral sources of lead such as lead-bearing metals imbedded in calcite can increase lead concentrations in solute and thus increase human health risk.73 Conversely, if no lead-bearing mineral exists, or no dissolution of the lead-baring mineral occurs, then there will be no additional risk to human health. Moreover, sorption and desorption processes, not simulated here can affect lead transport30,11 and will ultimately influence risk estimates. Thus properly representing aquifer mineralogy will determine the accuracy of human health risk estimates. In this study lead concentration at equilibrium with CO2 leak was approximately 0.33 ppb (ug L−1), whereas background lead concentration in equilibrium with the aquifer condition (Table 2) is 0.002 ppb. Dissolution and precipitation reactions in the σ2lnK = 1, and 3.61 ensembles reached equilibrium by the time reactive solute reached the pumping well as indicated by having Damköhler numbers equal to or greater than one. Therefore, the average lead concentration of streamlines in both σ2lnK = 1, and 3.61 ensembles at the pumping well is at or near 3.33 × 10−4 ppm (Figure 3c). Differences between estimated lead concentrations at the pumping well for the σ2lnK = 1, and 3.61 ensembles (Figure 3a) and the estimated human health risk (Figure 2) is then only due to differences in hydrological flux. However, streamlines in the σ2lnK =16 ensemble had a Damköhler number less than one and did not reach the new geochemical equilibrium. Consequently, the average streamline concentration for the σ2lnK =16 ensemble is 5.24 × 10−4 ppm. Therefore, estimated human health risk for the σ2lnK =16 ensemble is different from the σ2lnK = 1, and 3.61 ensemble (Figure 2) because of changes in both hydrological flux and the chemical state of solute arriving at the pumping well. Neglecting kinetic reactions rates in cases where the Damköhler number is less than one will result in a misrepresentation of human health risk. A comparison is made between a kinetic risk curve, produced by a fully coupled hydro-geochemical model and a risk curve that assumes chemical equilibrium is instantly met, where all streamlines have a lead concentration of 3.33 × 10−4 ppm (Figure 5). For the geochemical setup modeled here, assuming that chemical equilibrium is instantly met resulted in an underestimation of human health across all values of risk. For other geochemical scenarios that result in contaminant concentrations increasing as equilibrium is approached, risk might be overestimated.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 303-384-2456; fax: 303-273-3859; e-mail: rmaxwell@ mines.edu. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research has been supported by a grant from the US Environmental Protection Agency’s Science to Achieve Results (STAR) program. Although the research described in the article has been funded wholly or in part by the US Environmental Protection Agency’s STAR program through Grant RD83438701-0, it has not been subjected to any EPA review and therefore does not necessarily reflect the views of the Agency, and no official endorsement should be inferred. The authors would like to further thank the helpful and thorough reviews from John McCray, Magnus Skold, David Benson, and four anonymous reviewers.

Figure 5. Risk cuts for the 50th fractile of variability are plotted for the σ2lnK = 16 ensemble. The Kinetic σ2lnK = 16 cut is the estimated risk with kinetic reactive transport. The Equil σ2lnK = 16 cut is the estimated risk assuming the equilibrium concentration of 0.33 ppb (ug L−1) for all streamlines. G



■

Article

(19) de Barros, F. P. J.; Ezzedine, S.; Rubin, Y. Impact of hydrogeological data on measures of uncertainty, site characterization and environmental performance metrics. Adv. Water Resour. 2012, 36, 51−63, DOI: 10.1016/j.advwatres.2011.05.004. (20) Siirila, E. R.; Maxwell, R. M. Evaluating effective reaction rates of kinetically driven solutes in large-scale, statistically anisotropic media: Human health risk implications. Water Resour. Res. 2012(a, 48, W04527 DOI: 10.1029/2011WR011516. (21) Cirpka, O. A.; Frind, E. O.; Helmig, R. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol. 1999, 40, 159−182. (22) Dentz, M.; Gouze, P.; Carrera, J. Effective non-local reaction kinetics for transport in physically and chemically heterogeneous media. J. Contam. Hydrol. 2011, 120, 222−236, DOI: 10.1016/ j.jconhyd.2010.06.002. (23) Herrera, P. A.; Valocchi, A. J.; Beckie, R. D. A multidimensional streamline-based method to simulate reactive solute transport in heterogeneous porous media. Adv. Water Resour. 2010, 33, 711−727, DOI: 10.1016/j.advwatres.2010.03.001. (24) Yeh, T.-C. J. Stochasitc modeling of groundwater flow and solute transport in aquifers. Hydrological Processes 1992, 6, 369−359. (25) Rubin, Y. Applied Stochastic Hydrogeology, 1st ed.; Oxford University Press: Oxford, U.K., 2003 (26) Maxwell, R. M.; Carle, S. F.; Tompson, A. F. B. Contamination, risk, and heterogeneity: On the effectiveness of aquifer remediation. Environ. Geol. 2008, 54, 1771−1786, DOI: 10.1007/s00254-007-09558. (27) Wang, S.; Jaffe, P. R. Dissolution of a mineral phase in potable aquifers due to CO2 releases from deep formations; effect of dissolution kinetics. Energy Convers. Manage. 2004, 54, 2833−2848. (28) Carroll, S.; Hao, Y.; Aines, R. Geochemical detection of carbon dioxide in dilute aquifers. Geochem Trans. 2009, 10(4), DOI: 10.1186/ 1467-4866-10-4. (29) Kharaka, Y. K.; Thordsen, J. J.; Kakouros, E.; Ambats, G.; Herkelrath, W. N.; Beers, S. R.; Birkholzer, J. T.; Apps, J. A.; Spycher, N. F.; Zheng, L.; Trautz, R. C.; Rauch, H. W.; Gullickson, K. S. Changes in the chemistry of shallow groundwater related to the 2008 injection of CO2 at the ZERT field site, Bozeman, Montana. Environ. Earth Sci. 2010, 60, 273−284, DOI: 10.1007/s12665-009-0401-1. (30) Apps, J. A.; Zheng, L.; Xu, T.; Birkholzer, J. T. Evaluation of potential changes in groundwater quality in response to CO2 leakage from deep geological storage. Transport Porous Media 2010, 82, 215− 246. (31) Keating, E. H.; Fessenden, J.; Kanjorski, N.; Koning, D. J.; Pawar, R. The impact of CO2 on shallow groundwater chemistry: Observations at a natural analog site and implications of carbon sequestration. Environ. Earth Sci. 2010, 60, 521−536, DOI: 10.1007/ s12665-009-0192-4. (32) Little, M. G.; Jackson, R. B. Potential impacts of leakage from deep CO2 geosequestration on overlying freshwater aquifers. Environ. Sci. Technol. 2010, 44 (23), 9225−9232, DOI: 10.1021/es102235w. (33) Wilkin, R. T.; Digiulio, D. C. Geochemical impacts to groundwater from geologic carbon sequestration: Controls on pH and inorganic carbon concentrations from reaction path and kinetic modeling. Environ. Sci. Technol. 2010, 44 (12), 4821−4827. (34) Peter, A.; Lamert, H.; Beyer, M.; Hornbruch, G.; Heinrich, B.; Schulz, A.; Geistlinger, G.; Schreiber, B.; Dietrich, P.; Werban, U.; Vogt, C.; Richnow, H.; Großmann, J.; Dahmke, A. Investigation of the geochemical impact of CO2 on shallow groundwater: Design and implementation of a CO2 injection test in Northeast Germany. Environ. Earth Sci. 2012, DOI: 10.1007/s12665-012-1700-5. (35) Trautz, R. C.; Pugh, J. D.; Varadharajan, C.; Zheng, L.; Bianchi, M.; Nico, P. S.; Spycher, N. F.; Newell, D. L.; Esposito, R. A.; Wu, Y.; Dafflon, B.; Hubbard, S.; Birkholzer, J. T. Effect of Dissolved CO2 on a shallow groundwater system: A controlled release field experiment. Environ. Sci. Technol. 2012, DOI: 10.1021/es3011280t. (36) Lasaga, A. C. Rate laws in chemical reactions. In Kinetics of Geochemical Processes; Lasaga A. C., Kirkpatrick R. L., Eds.; Reviews in

REFERENCES

(1) EPA. Risk Assessment Guidance for Superfund (RAGS): Volume I Human Health Evaluations Manual (HHEM) (Part B, Development of Risk-Based Preliminary Remediation Goals), Technical report EPA 540/ R-92/003; U.S. EPA Office of Emergency and Remedial Response: Washington D.C., 1991. (2) EPA. Risk Assessment Guidance for Superfund (RAGS): Volume IIIPart A, Process for Conducting Probilistic Risk Assessment, Technical report EPA 540/R-02/002,; U.S. EPA Office of Emergency and Remedial Response : Washington DC, 2001. (3) Tartakovsky, D. M.; Nowak, W.; Bolster, D. Introduction to the special issue on uncertainty quantification and risk assessment. Adv. Water Resour. 2012, 36, 1−2, DOI: 10.1016/j.advwatres.2011.12.010. (4) Bogen, K. T.; Spear, R. C. Integrating uncertainty and interindividual variability in environmental risk assessment. Risk Anal. 1987, 7 (4), 427−436. (5) Mckone, T. E.; Bogen, K. T. Uncertainties in health-risk assessment − an integrated case-study based on tetrachloroethylene in California groundwater. Regul. Toxicol. Pharmacol. 1992, 15 (1), 83− 103. (6) Hoffman, F.; Hammonds, J. Propagation of uncertainty in risk assessments: The need to distinguish between uncertainty due to lack of knowledge and uncertainty due to variability. Risk Anal. 1994, 14 (5), 707−712. (7) Frey, H.; Rhodes, D. Characterization and simulation of uncertain frequency distributions: Effects of distribution choice, variability, uncertainty, and parameter dependence. Hum. Ecol. Risk Assess.: Int. J. 1998, 4 (2), 423−468. (8) Maxwell, R. M.; Kastenberg, W. E. Stochastic environmental risk analysis: An integrated methodology for predicting cancer risk from contaminated groundwater. Stochastic Environ. Res. Risk Assess. 1999, 13, 27−47. (9) Daniels, J. I.; Bogen, K. T.; Hall, L. C. Analysis of uncertainty and variability in exposure to characterize risk: Case study involving trichloroethylene groundwater contamination at Beal Air Force Base in California. Water, Air Soil Pollut. 2000, 123, 278−298. (10) Siirila, E. R.; Navarre-Sixtchler, A. K.; Maxwell, R. M.; McCray, J. E. A quantitative methodology to assess the risks to human health from CO2 leakage into groundwater. Adv. Water Resour. 2012, 36 (2), 146−164, DOI: 10.1016/j.advwatres.2010.11.005. (11) Siirila, E. R.; Maxwell, R. M. A new perspective on human health risk assessment: Development of a time dependent methodology and the effect of varying exposure durations. Sci. Total Environ. 2012(b), 431, 221−232, DOI: 10.1016/.scitotenv.2012.05.030. (12) Tartakovsky, D. M. Assessment and management of risk in subsurface hydrology. A review. Adv. Water Resour. 2013, 51, 247− 260, http://dx.doi.org/10.1016/j.advwatres.2012.04.007. (13) Maxwell, R. M.; Pelmulder, S. D.; Tompson, A. F. B.; Kastenberg, W. E. On the development of a new methodology for groundwater-driven health risk assessment. Water Resour. Res. 1998, 34 (4), 833−847. (14) de Barros, F. P. J.; Rubin, Y. A risk-driven approach for subsurface site characterization. Water Resour. Res. 2008, 44, W01414 DOI: 10.1029/2007WR006081. (15) de Barros, F. P. J.; Rubin, Y.; Maxwell, R. M. The concept of comparative information yield curves and its application to risk-based site characterization. Water Resources Research 2009, 45, DOI: 10.1029/2008WR007324. (16) Andričević, R.; Cvetković, V. Evaluation of risk from contaminants migrating by groundwater. Water Resour. Res. 1996, 32 (3), 611−621. (17) Maxwell, R. M.; Kastenberg, W. E.; Rubin, Y. A methodology to integrate site characterization information into groundwater-driven health risk assessment. Water Resour. Res. 1999, 35 (9), 2841−2855. (18) Benekos, I. D.; Shoemaker, C. A.; Stedinger, J. R. Probabilistic risk and uncertainty analysis for bioremediation of four chlorinated ethenes in groundwater. Stochastic Environ. Res. Risk Assess. 2007, 21, 375−390, DOI: 10.1007/s00477-006-0071-4. H



Article

Mineralogy & Geochemistry, Mineralogical Society of America, 1981; pp 135−169. (37) Lasaga, A. C. Chemical kinetics of water-rock interactions. J. Geophys. Res. 1984, 89, 4009−4025. (38) Maher, K. The dependence of chemical weathering rates in fluid residence time. Earth and Planetary Science Letters 2010, 294, 101− 110, DOI: 10.1016/j.epsl.2010.03.010. (39) Bearup, L. A.; Navarre-Sitchler, A. K.; Maxwell, R. M.; McCray, J. E. Kinetic metal release from competing process in aquifers. Environ. Sci. Technol. 2012, 46, 6539−6547, DOI: 10.1021/es203586y. (40) Sudicky, E. A. A natural gradient experiment on solute transport in a sand aquifer: Spatial variability of hydraulic conductivity and its role in dispersion process. Water Resour. Res. 1986, 22 (13), 2069− 2082. (41) Rubin, Y. Transport in heterogeneous porous media: Prediction and uncertainty. Water Resour. Res. 1991, 27 (7), 1723−1738. (42) McLaughlin, D.; Ruan, F. Macrodispersivity and large-scale Hydrologic Variability. Transp. Porous Media 2001, 42, 133−154. (43) Yabusaki, S. B.; Steefel, C. I.; Wood, B. D. Multidimensional, multicomponent, subsurface reactive transport in nonuniform velocity fields: Code verification using an advective reactive streamtube approach. J. Contam. Hydrol. 1998, 30, 299−331. (44) Atchley, A. L.; Maxwell, R. M.; Navarre-Sitchler, A. K. Using streamlines to simulate stochastic reactive transport in heterogeneous aquifers: Kinetic metal release and transport in CO2 impacted drinking water aquifers. Adv. Water Resour. 2013, 52, 93−106, DOI: 10.1016/ j.advwatres.2012.09.005. (45) Harvey, O. R.; Qafoku, N. P.; Cantrell, K. J.; Lee, G.; Amonette, J. E.; Brown, C. F. Geochemical implications of gas leakage associated with geological CO2 storageA qualitative review. Environ. Sci. Technol. 2012, 47, 23−36. (46) Viswanathan, H. S.; Pawar, R. J.; Stauffer, P. H.; Kaszuba, J. P.; Carey, J. W.; Olsen, S. C.; Keating, G. N.; Kavetski, D.; Guthrie, G. D. Development of a hybrid process and system model for the assessement of wellborer leakage at a geologic sequestration site. Environ. Sci. Technol. 2008, 42 (19), 7280−7286. (47) Stauffer, P. H.; Viswanathan, H. S.; Pawar, R. J.; Guthrie, G. D. A system model for geologic sequestration of carbon dioxide. Environ. Sci. Technol. 2009, 43 (3), 565−570. (48) Thiele, M. R.; Batycky, R. P.; Blunt, M. J.; Orr, F. M., Jr. Simulating flow in heterogeneous systems using streamtubes and streamlines. SPE Reservoir Eng. 1996, 11 (1), 5−12. (49) Crane, M. J.; Blunt, M. J. Streamline-based simulations of solute transport. Water Resour. Res. 1999, 35, 3061−3078. (50) Maxwell, R. M.; Welty, C.; Tompson, A. F. B. Streamline-based simulation of virus transport resulting from long term artificial recharge in a heterogeneous aquifer. Adv. Water Resour. 2003, 26, 1075−1096, DOI: 10.1016/S0309-1708(03)00074-5. (51) Simmons, C. S.; Ginn, T. R.; Wood, B. D. Stochasitc-convective transport with nonlinear reaction: Mathematical framework. Water Resour. Res. 1995, 31 (11), 2675−2688. (52) Ginn, T. R.; Simmons, C. S.; Wood, B. D. Stochastic-convective transport with nonlinear reaction: Biodegradation with microbial growth. Water Resour. Res. 1995, 31 (11), 2689−2700. (53) Batycky, R. P.; Blunt, M. J.; Thiele, M. R. A 3D field scale streamline-based simulator. SPE Reservoir Eng. 1997, 12, 246−254. (54) Evaluation of the Hydrologic Source Term from Underground Nuclear Tests in Frenchman Flat at the Nevada Test Site: The CAMBRIC Test, UCRL-ID-132300Tompson, A. F. B., Bruton, C. J., Pawloski, G. A, Eds.; Lawrence Livermore National Laboratory: Livermore, CA., 1999(a). (55) Tompson, A. F. B.; Carle, S. F.; Rosenberg, N. D.; Maxwell, R. M. Analysis of groundwater migration from artificial recharge in a large urban aquifer: A simulation perspective. Water Resour. Res. 1999(b), 35 (10), 2981−2998, DOI: 10.1029/1999WR900175. (56) Ginn, T. R. Stochastic-convective transport with nonlinear reactions and mixing: Finite streamtube ensemble formulation for multicomponent reaction systems with intra-streamtube dispersion. J. Contam. Hydrol. 2001, 47, 1−28.

(57) Tompson, A. F. B.; Bruton, C. J.; Pawloski, G. A.; Smith, D. K.; Bourcier, W. L.; Shumaker, D. E.; Kersting, A. B.; Carle, S. F.; Maxwell, R. M. On the evaluation of groundwater contamination from underground nuclear tests. Environ. Geol. 2002, 42, 235−247, DOI: 10.1007/s00254-001-0493-8. (58) Obi, E-OI.; Blunt, M. J. Streamline-based simulation of advective-dispersive solute transport. Adv. Water Resour. 2004, 27, 913−924. (59) Viswanathan, H. S.; Valocchi, A. J. Comparison of streamtube and three dimensional models of reactive transport in heterogeneous media. J. Hydraul. Research 2004, 42 (S1), 141−145, DOI: 10.1080/ 00221680409500057. (60) Obi, E-OI.; Blunt, M. J. Streamline-based simulation of carbon dioxide storage in a North Sea aquifer. Water Resour. Res. 2006, 42 DOI: 10.1029/2004WR003347. (61) Steefel, C. I.; Lasaga, A. C. A coupled model for transport of multiple chemical-species and kinetic precipitation dissolution reactions with application to reactive flow in single-phase hydrothermal systems. Am. J. Sci. 1994, 294 (5), 529−592. (62) Steefel, C. I.; Yabusaki, S. B. OS3D/GIMRT, Software for multicomponent-multidimensional reactive transport: User’s Manual and Programmer’s Guide, PNL-11166; Pacific Northwest National Laboratory: Richland, WA, 1996. (63) Steefel, C. I. Software for modeling multicomponent, multidimensional reactive transport. Lawrence Livermore National Laboratory, Livermore, Ca. UCRL-MA-143182. 2001. (64) EPA. Risk Assessment Guidance for Superfund (RAGS): Volume IHuman Health Evaluations Manual (HHEM) (Part A, Baseline Risk Assessment), Technical report EPA 540/1-89/002; U.S. EPA Office of Emergency and Remedial Response: Washington, DC, 1989. (65) Pollock, D. W. Semi-analytical computation of pathlines for finite-difference models. Ground Water 1988, 26 (6), 743−750. (66) Schafer-Perini, A. L.; Wilson, J. L. Efficient and accurate front tracking for two-dimensional groundwater flow models. Water Resour. Res. 1991, 27 (7), 1471−1485. (67) Zheng, L.; Apps, J. A.; Zhang, Y.; Xu, T.; Birkholzer, J. T. On mobilization of lead and arsenic in groundwater in response to CO2 leakage from deep geological storage. Chem. Geol. 2009, 268, 281− 297, DOI: 10.1016/j.chemgeo.2009.09.007. (68) Jennings, A. A. Critical chemical reactions rates for multicomponent groundwater contamination models. Water Resour. Res. 1987, 23 (9), 1775−1784. (69) Steefel, C. I.; Maher, K. Fluid-Rock Interaction: A Reactive Transport Approach; Reviews in Mineralogy & Geochemistry, Mineralogical Society of America, 2009; pp 485−532. (70) Brantley, S. L.; Mellott, N. P. Surface area and porosity of primary silicate minerals. Am. Mineral. 2000, 85 (11−12), 1767−1783. (71) Carroll, S. A.; O’Day, P. A.; Piechowski, M. Rock-water interactions controlling zinc, cadmium, and lead concentrations in surface waters and sediments, U.S. Tri-State mining district. 2. Geochemical Interpretation. Environ. Sci. Technol. 1998, 32, 956−965. (72) Carroll, S. A.; O’Day, P. A.; Esser, B.; Randall, S. Speciation and fate of trace metals in estuarine sediments under reduced and oxidized conditions, Seaplane Lagoon, Alameda Naval Air Station (USA). Geochem. Trans. 2002, 3, 81−101. (73) Navarre-Sitchler, A. K.; Maxwell, R. M.; Siirila, E. R.; Hammond, G. E.; Lichtner, P. R. Elucidating geochemical response of shallow heterogeneous aquifers to CO2 leakage using high-performance computing: Implications for monitoring of CO2 sequestration. Adv. in Water Resour. 2013, 53, 45−55, http://dx.doi.org/10.1016/j. advwatres.2012.10.005. (74) EPA. List of contaminants and there maximum contaminant level (MCLs), EPA 816-F-09-004, May 2009. http://water.epa.gov/ drink/contaminants/index.cfm#List.

I


Human Health Risk Assessment of CO2 Leakage into Overlying

Recommend Documents