Leakage of stored CO2 from a designated deep reservoir could contaminate overlying shallow potable aquifers by dissolution of arsenic-bearing minerals. To elucidate CO2 leakage-induced arsenic contamination, 2D multispecies reactive transport models were developed and CO2 leakage processes were simulated in the shallow groundwater aquifer. Throughout a series of numerical simulations, it was revealed that the movement of leaked CO2 was primarily governed by local flow fields within the shallow potable aquifer. The induced low-pH plume caused dissolution of aquifer minerals and sequentially increased permeabilities of the aquifer; in particular, the most drastic increase in permeability appeared at the rear margin of CO2 plume where two different types of groundwater mixed. The distribution of total arsenic (As) plume was similar to the one for the arsenopyrite dissolution. The breakthrough curve of As monitored at the municipal well was utilized to quantify the human health risk. In addition, sensitivity studies were conducted with different sorption rates of arsenic species, CO2 leakage rates, and horizontal permeability in the aquifer. In conclusion, the human health risk was influenced by the shape of As plume, which was, in turn, affected by the characteristics of CO2 plume behavior such as horizontal permeability and CO2 leakage rate.

1. Introduction

Carbon capture and storage are considered to be one of the mitigating strategies for reducing CO2 emissions to the atmosphere [13]. Among various carbon capture and storage technologies, CO2 can be injected into geologically stable formations, which typically have large storage capacities and are capped by low-permeability sealing formations. However, during CO2 injection activity, injection-induced pressure builds up within the storage formation [46]. If any unwanted pathways exist within the sealing formation, CO2 is able to migrate to the shallow aquifer through these pathways while experiencing phase change from supercritical to gaseous CO2 [710]. Especially in this study, the leakage of gaseous CO2 into a shallow confined aquifer was considered. Leaked gaseous CO2 dissolves in the potable groundwater and develops a low-pH plume [11, 12], which induces secondary contamination within the aquifer by enhancing the mobility of toxic heavy metals [1316].

Released toxic heavy metals are able to migrate with the ambient groundwater induced by natural hydraulic gradient or influenced by enforced gradient due to any active municipal wells nearby. If the contaminated groundwater produced from the municipal well is directly distributed without proper treatment (disinfection or chlorination processes) to residents who use it daily for the purpose of drinking, bathing, cleaning, or other household uses, these residents can be exposed to adverse carcinogenic health risks. The World Health Organization’s (WHO) International Agency for Research on Cancer classified various dissolved heavy metals (e.g., arsenic and lead) and other radioactive metals (e.g., uranium and cesium) as toxic substances hazardous to human health [17]. Among these carcinogenic heavy metals, this study focused on arsenic [18]. According to the US Agency for Toxic Substances and Disease Registry, who ranked hazardous substances based on their occurrence, toxicity, and potential for human exposure, arsenic was ranked the first in their Substance Priority List in 2013 and 2015 [19]. Moreover, arsenic contamination and its detrimental impacts have been reported by various countries such as China [20], Bangladesh [21], Vietnam [22], and India [23] in recent years.

Mobilization of arsenic in a shallow groundwater aquifer due to CO2 leakage has been investigated at one of the natural analog sites, Chimayó, New Mexico, where CO2-saturated brackish water was leaked into the shallow aquifer along the fault [7, 13, 24, 25]. At this site, decreased pH and the resulting mobilization of trace metals, including arsenic, were observed. Even if their adverse effects had been alleviated due to the high buffering capacity of the local groundwater aquifer, Keating et al. [7] reported significantly elevated trace metal concentrations at a number of local wells due to the influx of brackish waters. Later, both Keating et al. [13] and Viswanathan et al. [24] integrated the field dataset into a multiphase reactive transport model to understand the behavior of arsenic, since some wells in Chimayó exceeded the maximum contamination level (MCL). In addition to studies targeting natural CO2 release sites, several experiments have been conducted at field-scale CO2 injection sites to determine secondary contamination caused by the injected CO2 [11, 2628]. Decreased pH, increased concentrations, and subsequent changes in groundwater chemistry such as increased Fe2+, Mn2+, Mg2+, and Ca2+ concentrations were observed at both Frio-I brine pilot injection [11, 27] and ZERT field tests [12, 28].

In addition to field-oriented research, several studies have focused on the implementation of numerical simulations to evaluate geochemical behaviors associated with arsenic contamination. For example, Zhang et al. [29] and Xiao et al. [30] utilized the reactive transport simulation and investigated the complex chemical change induced by CO2 leaked into a shallow aquifer. In contrast to these studies which presented detailed geochemical interactions, Siirila et al. [31] simplified the geochemical processes by solving the advection-dispersion equation with a linear sorption. Without full assessment of geochemical reactions, they were able to account for the movement of toxic elements in complex 3D heterogeneous systems as well as for groundwater contamination-induced carcinogenic health risks with a probabilistic approach. Later, Navarre-Sitchler et al. [32] utilized PFLOTRAN and simulated the mobilization of lead in a complex heterogeneous system by assuming that released gaseous CO2 instantaneously dissolved in the groundwater aquifer. The proposed work in this study was built upon the framework of the previous studies mentioned above. Certain approaches accounted for the detailed geochemical behavior of toxic heavy metals (e.g., complexation, sorption, mineral dissolution, and precipitation) in the groundwater aquifer, whereas other approaches relied on either simplified chemical reactions (e.g., advection-dispersion equation) or multiphase fluid migration (e.g., gaseous CO2 leakage into the groundwater aquifer) while accounting for carcinogenic health risks or complex 3D heterogeneous systems. Therefore, the goal of this study was to integrate these two approaches and delineate the multiphase behavior of leaked gaseous CO2 to a shallow potable aquifer. Moreover, leaked CO2-induced geochemical changes, such as evolution of water quality and mobility of toxic trace metals (e.g., arsenic), were characterized by adapting the multispecies reactive transport model. Finally, the simulated concentration of arsenic species observed from the assigned municipal well was used to quantify the carcinogenic health risk for chronically exposed humans.

2. Behavior of Arsenic in Subsurface Environment

Arsenic in natural water typically originates from arsenic-bearing minerals that frequently possess sulfur, oxygen, and iron [33]. Generally, these naturally occurring arsenic-bearing minerals include arsenopyrite (FeAsS), realgar (AsS), enargite (CuAsS4), scorodite (FeAsO4·2H2O), and tennantite (Cu6[Cu4(Fe,Zn)2]As4S13) [29, 34, 35]. When these minerals dissolve, various forms of arsenic species such as , , , , , and can be released into groundwater. As seen in Figure 1, arsenic species can be present in several valence states (−3, 0, +3, and +5), but in natural groundwater, it is mostly found in oxyanions of trivalent arsenite (As(III)) or pentavalent arsenate (As(V)), depending on reducing or oxidizing conditions, respectively. Typically, the mobility and toxicity of As(III) are considered to be much higher than those of As(V) [3638]. Under the reducing condition, uncharged As(III) species, such as , is dominant below pH 9.2 (Figure 1). However, under the oxidizing condition, is dominant at pH < 6.9, while is dominant at higher pH (pH > 6.9). Moreover, and could exist under extremely acidic and alkaline conditions, respectively. In this study, when injected CO2 leaked into the shallow potable groundwater aquifer, the dissolution of leaked CO2 induced the acidification of the ambient groundwater to 4 < pH < 6 [11], where and were dominant under reducing and oxidizing conditions, respectively [33, 39]. Here, the shallow aquifer was considered to be under the reducing condition, and thus, was chosen as the most dominant species.

3. Work Flow: Numerical to Probabilistic Quantification

The following study comprised two major parts, namely, numerical prediction of CO2 transport from a leakage point to a municipal well (numerical simulation) and probabilistic quantification of health risks to humans who have been chronically exposed to a certain toxic heavy metal (health risk assessment), in sequence (Figure 2). Processes for both numerical simulation and probabilistic risk quantification were designed in three steps: “Data,” “Process,” and “Result.” During the “Data” step, input parameters were chosen selectively, and the sampling scheme was determined either deterministically or randomly. Subsequently, a series of calculations using either numerical or probabilistic approaches were conducted in the “Process” step. Finally, in the “Result” step, simulation outputs, such as CO2 plume distribution, maximum concentrations of selected toxic heavy metals, and human health risk, were analyzed. The link between “numerical simulation” and “health risk assessment” was the profile of arsenic concentration monitored at the municipal well; throughout complex spatial and temporal movements of the CO2 plume predicted from numerical simulations, a time series of arsenic concentration was observed at the municipal well, which was then used to assess carcinogenic human health risk.

4. Numerical Approach

The multiphase and multicomponent reactive transport simulator, TOUGHREACT, was used to simulate secondary contamination processes in a shallow confined aquifer induced by CO2 leakage [44] in conjunction with the ECO2N module, which was used to predict the fluid properties of H2O, NaCl, and gaseous CO2 [45]. For this study, the original thermodynamic database incorporated in TOUGHREACT was not adequate for delineating the complex chemical reactions associated with arsenic species. Therefore, the revised EQ3/6V7.2b database was specifically adopted to account for chemical reactions of arsenic and related chemical species [43]. In detail, the modified thermodynamic database adopted arsenite () as a primary species and incorporated associated aqueous complexes (e.g., , , H3AsO4(aq), HAsO2(aq), , and HAsS2(aq)), while accounting for their activity coefficients from the extended Debye-Hückel equation [46]. Moreover, solubility products of arsenic-bearing minerals such as arsenopyrite were included in this thermodynamic database.

Precipitation and dissolution of minerals were simulated kinetically by following the rate law, which was coupled with the equation representing the kinetic rate constant (see (1)). In this equation, the kinetic rate constant is dependent not only on temperature but also on the pH as shown below [47]:where  (J/mol) is the activation energy and  (mol/m2/s) is the rate constant at 25°C with superscripts nu, , and representing neutral, acid, and base mechanisms, respectively.  (J/mol/K) is the gas constant,  (K) is temperature, is the activity of dissolved species under acid or base conditions, and is the power term. is the specific reactive surface area (cm2/g), and is the kinetic mineral saturation ratio of the th mineral. Finally, both and are assumed to be unity.

5. Conceptual Model

5.1. Model Description

The model designed for this study delineates the CO2 leakage process throughout undetected and unexpected pathways in the sealing formation. In order to understand such processes, the 2D cross-sectional potable aquifer model was designed as shown in Figure 3. The potable aquifer was assumed to be relatively deep (300 m with a reducing condition), where the municipal or high-capacity wells produced a large amount of groundwater, which was then distributed for residential purposes [48, 49]. In addition, the width and thickness of the potable aquifer were assigned as 200 m and 40 m, respectively, with the size of grid blocks of 2 m; the total number of grid blocks was 2,000 (100 × 20). Assuming that the aquifer is located at a depth of 300 m, the initial formation pressure and temperature were assigned as 3 MPa and 25°C, respectively. Top and bottom boundaries of the model were assigned as no-flow conditions assuming that the upper and lower confining formations act as nearly impermeable sealing units. Lateral boundaries (purple-colored grid blocks) were assigned as the Dirichlet condition, where the left and right pressure were 3.2 MPa and 2.8 MPa, respectively. The difference of pressure in lateral boundaries in addition to aquifer properties ( m2) results in 1.94 cm/day of the ambient groundwater flow in the aquifer. The leakage point for gaseous CO2 was located at a 50 m distance from the left boundary where CO2 was leaked at a rate of 0.05 kg/s (Figure 3). The CO2 leakage was maintained for 1 year. Furthermore, the municipal well with a continuous pumping rate of 0.5 kg/s was located 100 m away from the CO2 leakage point and penetrated a depth of 30 m from the upper confining seal. Finally, the simulation was conducted for 100 years.

5.2. Physical and Chemical Parameters

The aquifer was assumed to be homogeneous with its horizontal permeability () and porosity () to be 10–13 m2 and 0.2, respectively (Table 1); the anisotropy ratio (/) was assumed to be 0.1. In addition, delineation of gaseous CO2 transport through the groundwater-saturated aquifer required constitutive relations such as relative permeability and capillary pressure. In this study, Van Genuchten’s functions were adapted, and the relevant parameters are shown in Table 1 [40].

The mineralogical composition of the designated aquifer materials was assumed to be sandstone, which is the typical hosting formation for geological CO2 sequestration [50, 51] and natural analog CO2 leakage sites [52, 53]. As an example, Kampman et al. [54] analyzed the mineralogical assemblage of Navajo Sandstone, which was considered to be the primary sourcing aquifer for CO2-charged brine [55]. Chemical analyses of Navajo Sandstone fluids collected either from adjacent geysers or from springs revealed elevated concentrations of both arsenic and other toxic heavy metals [54, 5658]. Due to this reason, mineralogical composition in this study was adopted from the composition of Navajo Sandstone while assuming that arsenopyrite is the primary source of arsenic species in the ambient groundwater. A total of 7 primary minerals, that is, quartz, kaolinite, illite, K-feldspar, arsenopyrite, calcite, and magnesite, were chosen, and another 5 secondary minerals (chlorite, dolomite, goethite, oligoclase, and smectite-Ca) were expected to be precipitated (Table 2).

Mineral volume fraction of each primary mineral was chosen from a similar quantity of Navajo Sandstone with the addition of 1% arsenopyrite; quartz was predominant (81%), and kaolinite accounted for the second-largest amount (16%). Due to the small proportion of carbonate minerals, fluid chemistry possessed less buffering capacity against CO2 dissolution. With the chosen mineral composition, the batch reaction was conducted to determine a list of primary species and their initial concentrations. As a result, , Ca2+, Cl, Fe2+, H+, H2O, H3AsO3(aq), , , HS, K+, Mg2+, Na+, O2(aq), SiO2(aq), and were selected as primary species, and associated aqueous complexes were chosen as secondary species (Table 3).

As described by Keating et al. [7] and Zheng et al. [43], the concentration of arsenic species in the aquifer media was significantly controlled by both adsorption and desorption processes. Due to this reason, many researchers experimentally measured sorption values () for arsenic species under various conditions (e.g., different minerals and pH) and revealed that of As(III) and As(V) varied from 1.47 L/kg to 275 L/kg within geologic media including subsurface soil and aquifers [39, 59, 60]. In this study, a shallow sandstone aquifer was targeted, composed of over 90% quartz and kaolinite (Table 2). Previous experimental measurements revealed that for quartz and kaolinite was measured at 2 and 19 L/kg, respectively [39]. Therefore, the value was assigned as 10 L/kg.

The multiple parameters required to address the kinetic rates for mineral reactions following (1) are listed in Table 4. In addition, the calculation of reactive surface area for minerals was followed by both Xu et al. [27] and Sonnenthal et al. [61], who assumed that a mineral is a cubic array of truncated spheres, in which the radius of the sphere is assumed to be 0.001 m. In this study, the surface roughness-based area predicted from the spherical radius was reduced by two orders to reasonably represent the reactive surface area. Typically, chemical reactions only occur at selected sites on the mineral surface, and furthermore, only a small fraction of mineral surface is involved in this reaction due to grain coating and armoring. Therefore, the reactive surface areas of most silicate and carbonate minerals were chosen to be approximately 10 cm2/g, similar to those chosen by Knauss et al. [62] and Zerai et al. [63]. Finally, reactive surface areas of clay minerals such as kaolinite, illite, chlorite, and smectite were selected for larger values because of their smaller grain sizes [43].

Porosity was associated with changes in volume due to mineral dissolution or precipitation [64]. The porosity in this study was simulated by the following equation:where is the number of minerals and and are the volume fractions of th mineral in the rock and nonreactive rock, respectively. Finally, permeability change was calculated using the porosity changes with cubic law [65]:where and are permeability and porosity, respectively, with subscript representing the initial value.

Even if the proposed approach was capable of simulating water-rock interaction, it would still have some limitations. Firstly, in this study, the sorption effect of arsenic onto the surface of clay minerals was simulated with a linear approach instead of surface complexation. Although the linear approach had limitations in terms of delineating chemical heterogeneities on both temporal and spatial scales, this approach was effective for application to a large-scale simulation by reducing computation [6668]. Secondly, arsenopyrite is a solid solution of pyrite (FeS2), in which the ratio of arsenic and sulfur varies depending on their mole fractions. However, in this study, for simplification in predicting the thermodynamic properties, it was assumed that 1 mole fraction of arsenic replaced the sulfur (FeAsS). Finally, it was assumed that arsenopyrite oxidatively dissolves in the presence of common geologic oxidants such as dissolved O2, , and Fe3+. Under acidic conditions, Fe3+ quickly oxidizes, and sequential dissolution of arsenopyrite releases arsenic species ( and ) by following (4). Therefore, the stoichiometric reaction for arsenopyrite dissolution, which was used in this numerical simulation, involved the reduction of Fe3+ to Fe2+ [43].

5.3. Probabilistic Health Risk Approach

Carcinogenic effects on humans who have been chronically exposed to arsenic species through multiple pathways were probabilistically quantified based on individual exposure rate and toxicity, which were suggested by the “Guidelines for Carcinogenic Risk Assessment” of the US Environmental Protection Agency (EPA) [69], Siirila et al. [31], and the EPA Superfund Risk Assessment [70]. In general, three different uptake pathways, namely, ingestion, dermal sorption, and inhalation, were considered. Among these, ingestion and dermal sorption were considered to be major pathways, because humans are often exposed to risk by drinking dissolved toxic species in tap water or showering. Uptake through inhalation was not considered in this study, as it is unlikely that the concentration of vaporized trace metals was high enough to cause carcinogenic effects indoors under normal conditions. Therefore, following Siirila et al. [31], only two exposure pathways (ingestion and dermal sorption) were considered for quantifying exposure rate and toxicity of arsenic.

Arsenic toxicity was predicted from the product of cancer potency factor (CPF) (kg·day/mg) and average daily dose (ADD) (mg/kg·day). Typically, the CPF is different at the individual pathway even for the same toxic metal; in this study, the CPF for ingestion and dermal sorption of arsenic were assumed to be 1.5 and 1.6 kg·day/mg by following IRIS [71] and EPA [72], respectively. Subsequently, ADD was used to assess the individual exposure rate with the following equation:where is the individual intake per body weight (L/kg/day), AT is the average lifetime, which was assumed to be 25,550 days (70 years), and EF is the standard exposure frequency during 1 year, which was assumed to be 350 (days/year) [73]. The most significant term in (6) is , which is the maximum average of arsenic concentration (mg/L) monitored at the municipal well during the exposure duration (ED) (years). In this study, the ED was assumed to be 30 years following the EPA guideline [74]. Accordingly, , calculated from concentration profiles () of arsenic species at the municipal well in the numerical simulation can be represented as follows:

Specific to ADD, exposure pathways through both ingestion and dermal sorption are defined as follows:where is the skin surface area per body weight (m2/kg), is the dermal permeability coefficient for arsenic in water, designated as 1.0 × 10−5 m/hour, is the fraction of skin in contact with water (−), is exposure time of shower per day (hours/day), and CF is the unit conversion factor (0.001 L/m3). Finally, the total risk is the summation of individual risks representing each exposure pathway as shown below:In this study, the calculation of the total risk was repeated 100,000 times to account for variability, which stands for probabilistic quantification in human health risk. Except for (the concentration of arsenic species obtained from the result of numerical simulation), individual parameters considered in the quantification of total risk were randomly sampled within each intrinsic distribution (Table 5) [31, 75]. Finally, the calculated total health risks were plotted as a cumulative density function, which enables the estimation of the probability exceeding the risk level of concern (10−4) [76, 77].

6. Model Scenarios

Case 1 (base case) was designed to delineate CO2 leakage processes associated water-rock interactions and secondary contamination caused by arsenic species (Table 6). Subsequently, sensitivity studies were conducted with different , CO2 leakage rates, and horizontal permeability of the aquifer (). First, the degree of sorption intensity was evaluated by varying between 25, 50, and 100 L/kg (Cases 2–4). Second, in Cases 5–7, the effect of CO2 leakage rate was evaluated by varying its rate between 0.020, 0.025, and 0.030 kg/s. Different CO2 leakage rates could induce the development of CO2 plumes with different sizes. For example, as the size of a CO2 plume increases, larger areas expect to experience water-rock interactions and more dissolution of arsenic species. Consequently, the municipal well captures dissolved arsenic species more when the size of the CO2 plume is greater, which eventually increases the carcinogenic health risk on humans. Finally, in Cases 8–10, varied ( = 0.2 × 10−13, 1.0 × 10−13, and 5.0 × 10−13 m2) while maintaining the vertical permeability (1.0 × 10−14 m2); increased accelerates the horizontal velocity of ambient groundwater while reducing buoyancy forces on the CO2 plume.

7. Results and Discussion

7.1. Base Case
7.1.1. Migration of Leaked CO2 Plume within the Shallow Potable Aquifer

Figures 4(a)4(e) represent the evolution of the leaked CO2 plume at designated times of 120, 240, 360, 480, and 600 days. The mass centers of the CO2 plume, shown as red, black, and yellow circles, were calculated and plotted every 120 days; black circles represent the present time of the mass center, and red and yellow circles represent past and future times, respectively. From the leakage point, CO2 continuously leaked at a rate of 0.05 kg/s only until 365 days. During this period, three flow systems induced by ambient groundwater, CO2 leakage, and pumping activity interacted with each other (Figure 4(c)), which developed two mixing zones for geochemically different types of groundwater at both the front and the rear margins of the CO2 plume. The ambient groundwater flow (orange arrows) was developed from the left to the right boundaries at an approximate rate of 1.94 cm/day. Additionally, the CO2 plume gradually expanded from the leakage point where CO2 saturation remained at 0.3; the rate of CO2 flux was approximately 2.1 × 10−5 kg/(s·m2) (black arrows) adjacent to the leakage point. At the rear margin of the CO2 plume, two chemically different types of groundwater (ambient groundwater and CO2-dissolved groundwater) flowed in opposite directions, inducing the development of a vigorous geochemical mixing zone. Subsequently, the CO2 plume migrated together with the ambient groundwater until it was captured by the municipal well (Figure 4(c)).

After CO2 leakage had stopped at 365 days, the size of the CO2 plume gradually decreased due to both dissolution to ambient groundwater and extraction from the municipal well (Figures 4(d) and 4(e)); CO2 solubility predicted by Duan and Sun [78] was 0.83 mol/kg water in this aquifer (3 MPa and 25°C). After movable CO2 was captured by the municipal well, residually trapped CO2 governed by irreducible CO2 saturation (shown in Table 1) remained until complete dissolution to the ambient groundwater; residually trapped CO2 eventually vanished approximately 5 years after CO2 leakage had stopped.

The calculated mass center was located close to the plume center while the CO2 plume migrated toward the municipal well. The migration rate of the CO2 plume estimated from the mass center locations was approximately 12.5 cm/day until 480 days (Figures 4(a)4(d)), indicating that the migration rate of the CO2 plume was faster than the ambient groundwater flow (1.94 cm/day). This is because the pumping activity developed an additional head gradient, which was greater than ambient groundwater flow. After 480 days, the mass center slightly moved back until 600 days (Figure 4(e)), implying that all movable CO2, which was weighted CO2 mass at the plume front, was pumped out. After 600 days, continuous CO2 dissolution occurred at the plume rear, and the pumping activity at the front induced the movement of the mass center slowly.

7.1.2. Induced Geochemical Reactions

The dissolution of gaseous CO2 into ambient groundwater increased concentration from 0.2 to 1 mol/L while decreasing the pH from 8.0 to 5.5 within the CO2 plume (Figures 5(a) and 5(d)). Concurrently, both carbonate and silicate minerals were either dissolved or precipitated. In particular, dramatic changes in minerals were localized at both rear and front margins of the CO2 plume where the vigorous advective mixing of chemically different types of groundwater occurred.

Dissolutions in calcite (CaCO3) and magnesite (MgCO3) were distinct (Figures 5(b) and 5(c)). Such dissolution released into the ambient groundwater and subsequently induced positive feedback to lower pH. Calculation of the saturation index (SI) using the initial concentration of ambient groundwater indicated that calcite () was more saturated than magnesite (); the initial concentrations of Ca2+ and Mg2+ were 3.3 × 10−3 mol/L and 4.2 × 10–12 mol/L, respectively, in the ambient groundwater. Due to this reason, when CO2 was leaked, dissolution of magnesite (−6 mol/m3) was greater than that of calcite (−4 mol/m3) within the CO2 plume. However, even if overall magnesite dissolution was greater than that of calcite, the greatest change in mineral dissolution appeared in calcite (−8 mol/m3), focusing at the rear margin of the CO2 plume (Figure 5(b)). This implies that the localized dissolution of calcite was primarily induced by the mixing of two chemically different types of groundwater, such as the ambient and CO2-dissolved groundwater. Initially, the Ca2+ concentration in the ambient groundwater was 3.3 × 10−3 mol/L (Figure 5(e)). Inside the CO2 plume, calcite was dissolved and increased Ca2+ concentration 10-fold to 2.5 × 10−2 mol/L. At the rear margin of the CO2 plume, Ca2+ concentration was increased even more (to 3.1 × 10−2 mol/L). Overall, the distribution of Ca2+ was similar to that of calcite (Figures 5(b) and 5(e)), and distributions of both Mg2+ and were similar to that of magnesite (Figures 5(c), 5(d), and 5(f)).

Patterns of dissolution and precipitation in silicate minerals were more complex than those of carbonate minerals (Figures 6(a)–6(c)). CO2 leakage primarily induced the dissolution of K-feldspar (KAlSi3O8) (Figure 6(a)); the greatest dissolution (3.5 × 10−2 mol/m3) occurred at the rear of the CO2 plume, and the degree of dissolution gradually decreased as the plume approached the municipal well. Dissolution of K-feldspar increased concentrations of K+, SiO2(aq), and in the groundwater (Figures 6(d)–6(f)). Nevertheless, distributions of such species did not imitate the dissolution pattern of K-feldspar. Rather, SiO2(aq) and K+ showed the highest concentrations at the rear and front margins of the CO2 plume with values of 1.8 × 10−4 mol/L and 3.5 × 10−4 mol/L, respectively (Figures 6(d) and 6(e)), but the decrease in concentration occurred uniformly throughout the CO2 plume (Figure 6(f)). The discrepancy in patterns between K-feldspar and other dissolved species was presumably caused by a combination of both dissolution and precipitation among various silicate minerals such as illite, kaolinite, and chlorite as described below.

For the distribution of illite (K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2), a small amount (3.3 × 10−2 mol/m3) was precipitated throughout the CO2 plume (Figure 6(b)). However, at the rear margin and immediately adjacent to the municipal well, a relatively large degree of illite dissolution was predicted (−9.1 × 10−3 and −7.4 × 10−2 mol/m3, resp.). In contrast to illite, a small amount (−2.9 × 10−2 mol/m3) of kaolinite was dissolved within the plume, and reversely, a small amount of precipitation was predicted at the rear and front margins (1.9 × 10−2 and 4.8 × 10−2 mol/m3, resp.) (Figure 6(c)). Presumably, the dissolution and precipitation of illite and kaolinite would influence the distribution of SiO2(aq) and K+ in addition to K-feldspar dissolution (Figures 6(d) and 6(e)). Finally, the behavior of individual silicate mineral influenced the distribution of , the concentration of which within the plume was lower than that outside (Figure 6(f)). Overall, CO2 leakage induced dissolution or precipitation of both carbonate and silicate minerals and ultimately changed both the porosity and the permeability of the shallow aquifer. The dissolution of carbonate minerals primarily caused an increase in permeability; permeability increased to 0.43% () within the CO2 plume, and the most drastic increase (0.77%) occurred at the rear margin of the CO2 plume (Figure 7).

Dissolution of arsenopyrite, which was the primary reaction for predicting carcinogenic health risk, occurred only within the CO2 plume with a dissolved amount of 4.84 × 10−5 mol/m3 (Figure 8(a)). Similar to that of carbonate minerals, the greatest amount of arsenopyrite dissolution (8.57 × 10−5 mol/m3) occurred at the rear margin of the plume due to the vigorous mixing of two chemically different types of groundwater. Following (4), oxidative dissolution of arsenopyrite consumed 0.75 moles of O2(aq) and 1 mole of H+ while increasing the concentrations of total arsenic (As), Fe2+, and HS (Figures 8(b)8(e)). In this study, As represents the summation of primary species such as arsenite ((aq)), which is the by-product of arsenopyrite dissolution, as well as other arsenic species such as , H3AsO4(aq), , , , HAsO2(aq), and HAsS2(aq). Distribution of As concentration mimics that of arsenopyrite (Figures 8(a) and 8(b)); generally, As concentration within the CO2 plume was greater than that outside, while the rear margin revealed the highest concentration. However, the other associated species such as Fe2+, HS, and O2(aq) revealed relatively uniform distribution (Figures 8(c)8(e)). The difference between arsenopyrite-produced species such as As and other associated species (e.g., Fe2+, HS, and O2(aq)) presumably occurred due to the sorption effect, which was accounted with the linear approach. was designated for arsenic species only, and thus, as shown in Figure 8(b), enrichment of As concentration occurred at the rear of the CO2 plume. In summary, the greatest concentration of As was 4.9 × 10−7 mol/L at the rear margin of the plume and the average concentration of As within the plume was 2.9 mol/L (Figure 8(b)). The average concentrations of other species such as Fe2+ and HS were 3.9 × 10−7 mol/L and 3.6 × 10−7 mol/L, respectively (Figures 8(c) and 8(d)).

7.1.3. Health Risk Assessment of Carcinogenic Effect

To account for carcinogenic health risk, selected species such as pH, As, and arsenite ((aq)) concentrations were monitored at the municipal well for 100 years (Figure 9(a)). Dramatic changes in gaseous CO2 saturation (), pH, and mass fraction of CO2 dissolved in groundwater () predicted during 10 years were magnified at the small window (Figure 9(a)). Depending on the profiles of both pH and dissolved arsenic species, two stages (Stages I and II) were characterized. During Stage I (0–6.3 years), the leaked CO2 plume, which existed as either gaseous CO2 (, black dotted line) or dissolved CO2 (, purple dotted line), arrived at the municipal well approximately after 360 days. Once the CO2 plume arrived at the well, immediate reduction of pH from 8.8 to 5.2 was observed, while As and (aq) concentrations sharply increased to 2.93 × 10−7 and 1.32 × 10−7 mol/L, respectively, exceeding the maximum contaminant level (MCL = 1.33 × 10−7 mol/L, red dotted line) [79]. After 2.1 years, decreased to 0 at the municipal well, implying that all movable gaseous CO2 was pumped out. Even after all gaseous CO2 had been diminished due to the pumping activity, residually trapped CO2 remained within the pores while dissolving into the groundwater. Due to the dissolution of residually trapped CO2, the mass fraction of dissolved CO2 () was invariant at 0.035 until 5 years; in this shallow aquifer, flow caused by both ambient fresh groundwater and pumping activity accelerated CO2 dissolution. Therefore, complete dissolution of residually trapped CO2 appeared at 6.3 years when became 0 and pH returned to 8.2. In addition, concentrations of As and (aq) reached 3.77 × 10−7 and 1.70 × 10−7 mol/L, respectively. Immediately beginning with Stage II (6.3–100 years), elevated pH (8.2) of groundwater inhibited the dissolution of arsenopyrite, the reaction of which requires the consumption of H+ (see (4)). Due to decreased arsenopyrite dissolution, concentrations of both As and (aq) were stabilized and gradually decreased due to sorption on aquifer media. The pH continuously decreased until 90 years. However, concentrations of As and (aq) reached the background level (6.60 × 10−8 and 3.01 × 10−8 mol/L, resp.) at 62 years.

Based on the simulated profiles of both As and (aq), , which is the maximum average concentration calculated from (7), was predicted, and the carcinogenic health risk was quantified following the method described in Section 5.3. Figure 9(b) shows histograms representing calculated frequencies of carcinogenic risk for both As (blue bar) and (aq) (green bar). The carcinogenic risk predicted from (aq) profile, which revealed relatively low concentrations, showed a mean, median, and standard deviation of 4.00 × 10−4, 4.30 × 10−4, and 1.69 × 10−4, respectively. For risk predicted from As concentration profile, the mean, median, and standard deviation were 8.94 × 10−4, 9.60 × 10−4, and 3.77 × 10−4, respectively. Cumulative density functions (blue and green lines) were also plotted together with the risk level of concern (10−4) [77, 80]. From the cumulative density functions, the risk, which exceeds the risk level of concern, can be considered to have carcinogenic potential after chronic exposure to arsenic-contaminated groundwater. As shown, the risk level predicted from both As and (aq) exceeded the risk level (red dotted line) of concern.

7.2. Sensitivity Studies
7.2.1. Effect of Sorption Intensity () (Cases 2–4)

In this sensitivity study, the intensity of the sorption effect ( of 25, 50, and 100 L/kg) on As was evaluated while gaseous CO2 was leaked into the shallow aquifer (Table 6). As shown in Figure 4, once gaseous CO2 was leaked from unidentifiable pathways, it migrated with the ambient groundwater. Here, characteristics of gaseous CO2 plume such as its size, shape, migrating velocity, and gas saturation () were influenced by multiphase parameters such as capillary pressure and relative permeability (Table 1); the role of in gaseous CO2 transport was minimal because the sorption typically accounted for the movement of dissolved species within geologic media. Due to this reason, regardless of the variation in , the distribution and behavior of gaseous CO2 plume were not affected, and therefore the CO2 plume remained essentially the same in all cases as shown in Figure 4. However, within the gaseous CO2 plume, various geochemical processes occurred, including reductions in pH, dissolution of arsenopyrite, and interactions between dissolved species. In particular, variation in was anticipated to affect the behavior of dissolved arsenic species after all gaseous CO2 was pumped out or dissolved into the groundwater.

In Stage I, the gaseous CO2 plume arrived at the municipal well after 350 days (0.95 years) for all cases (Figure 10(a)). The elevated concentration (2.9 × 10−7 mol/L) of As occurred simultaneously in all cases because the source of arsenic species was the dissolution of arsenopyrite, which resulted from CO2 dissolution. Despite the variation in , the arrival times of As at the municipal well were the same because the sorption did not affect the migration of the gaseous CO2 plume. Therefore, the As profiles evolved similarly until approximately 30 years (the middle of Stage II) when its concentration reached a maximum (3.9 × 10−7 mol/L). After the arrival of the As peak, the differences between As profiles were amplified until the As concentration reached the background level (6.7 × 10−8 mol/L). In detail, differences in the slopes of As profiles were small immediately after the arrival of the As peak (30–40 years), but the discrepancy was amplified from 40 years while the slopes for As profiles sharply dropped. Differences in As profiles were attributed to the degree of , which determined the amount of arsenic adsorbed to aquifer media, especially at the rear margin of the CO2 plume; at this location, the highest As concentration occurred due to the mixing of two chemically different types of groundwater as shown in Figure 8(b).

In these simulations, with increasing , more arsenic was adsorbed to the aquifer media, and the migration of As was therefore retarded. In other words, stronger retardation caused As concentration to be maintained higher and longer in the aquifer, and therefore the arrival of As concentration at the background level was delayed. For example, in Case 2 ( = 25 L/kg), As concentration revealed the earliest recovery (69.3 years) at the municipal well (Figure 10(a)). As increased to 50 and 100 L/kg, the recovery time was delayed to 72.1 and 73.9 years, respectively.

Figure 10(b) represents the predicted probabilistic health risk for Cases 2–4. While calculating the health risk, , the peak calculated from the moving average of As concentration by adopting the designated interval of ED (30 years), influenced the health risk most significantly based on (8), (9), and (10). Since As profiles revealed similar patterns with the same peak values while the only difference being the recovery time, the calculated values for Cases 2–4 were almost the same (3.9 × 10−7 mol/L). Consequently, predicted carcinogenic health risk for humans was almost identical to variation in . This result implies that variation in was a less influential parameter for assessing health risk for arsenic species. This was because variation in did not affect the characteristics of the CO2 plume such as size, shape, and migration velocity, which determined the dissolved amount or reaction rate for arsenic sources such as arsenopyrite. Due to this reason, additional simulations were conducted and described in the following section after varying parameters (CO2 leakage rate and aquifer permeability), which directly affected the size of the CO2 plume and its migration rate.

7.2.2. Effect of CO2 Leakage Rate () (Cases 5–7)

In Cases 5–7, the effect of CO2 leakage rate () on the quantification of carcinogenic health risk was evaluated. Figures 11(a)–11(c) show the distribution of the gaseous CO2 plume after CO2 leakage was stopped at 1 year, and Figures 11(d)–11(f) present the distribution of As concentration after 20 years. As increased from 0.020 to 0.025 to 0.030 kg/s, the CO2 plume approached the municipal well more closely; the calculated mass center of the CO2 plume, plotted as black circles at a 120-day interval, moved rapidly as increased. In addition, adjacent to the CO2 leakage point, CO2 saturation was elevated from 0.27 to 0.30, implying occurrence of active dissolution of gaseous CO2. As described before, although the movable gaseous CO2 plume was pumped out after approximately 2 years, residually trapped CO2 still remained in the pores, concurrently dissolving into the ambient groundwater until 8 years. Increases in dissolved CO2 concentration due to dissolution of residually trapped CO2 and resulting low-pH groundwater influenced the distribution of As concentration even after all gaseous CO2 was pumped out (Figures 11(d)–11(f)). For example, the intensity of governed the size of the CO2 plume where the active dissolution of arsenopyrite, the source of As, primarily occurred. Therefore, as increased from 0.020 to 0.025 to 0.030 kg/s, the size of As plume expanded at 20 years. Additionally, the effects of dispersion, diffusion, and sorption were amplified while the As plume migrated in the ambient groundwater.

The As profiles monitored at the municipal well revealed the drastic increase in As concentration immediately after the arrival of the As plume at 1.31, 1.07, and 0.98 years for Cases 5, 6, and 7, respectively (Figure 12(a)). The arrival time of the As plume coincided with that of the gaseous CO2 plume, implying that multiphase migration of CO2 governed the movement of dissolved As. The size of the As peaks, all of which were over the MCL, increased with (1.56 × 10−7, 2.75 × 10−7, and 4.07 × 10−7 mol/L for Cases 5, 6, and 7, resp.). In addition, its peak was maintained for a longer time with increased (13.5, 16.6, and 23.0 years for Cases 5, 6, and 7, resp.). Figure 12(b) shows the calculated probabilistic health risk. While calculating the health risk using (9) and (10), an important parameter was , which was the maximum average of arsenic concentration monitored at the municipal well (see (7)). The size of varied with as shown in Figure 12(a), which directly affected the risk prediction. Therefore, as the As concentration increased, the health risk for humans increased accordingly. In particular, the medians of Cases 5, 6, and 7, located at the half percentile in variability, were 3.39 × 10−4, 6.53 × 10−4, and 9.56 × 10−4, respectively, all of which exceeded the risk level of concern (10−4).

7.2.3. Effect of Horizontal Permeability () (Cases 8–10)

From two previous sensitivity studies, it was revealed that the driving force of the As plume was essentially the movement of the gaseous CO2 plume; depending on the size of the CO2 plume, the amount of dissolved As was determined. In this study, the magnitude of horizontal permeability (), which altered the velocity of ambient groundwater flow, varied between 0.2 × 10−13, 1.0 × 10−13, and 5.0 × 10–13 m2 (Cases 8–10, resp.). As increased, the ambient groundwater rate increased from 0.39 to 1.94 to 9.71 cm/day, which induced a change in the shape of the CO2 plume from oval to flat (Figures 13(a)–13(c)).

Typically, the shape and movement of the gaseous CO2 plume within the groundwater aquifer are governed by the balance between buoyancy and viscous forces due to the density contrast of these two fluids. Even the density contrast is amplified as CO2 leaks into the shallower aquifer, because CO2 density radically decreases while approaching the surface. The magnitude of the buoyancy number () reflects the change in CO2 plume shape, where is gravitational acceleration, and are the length and height of the model domain, respectively, is the ambient groundwater flow rate, and is dynamic viscosity of CO2 (1.47 × 10−5 Pa·s). Finally, is the density difference between groundwater (1,000 kg/m3) and CO2 (1.842 kg/m3) [81]. The calculated of Cases 8, 9, and 10 was 740.30, 148.06, and 29.61, respectively, implying that the buoyancy force acting on CO2 plume increased as decreased.

Even if the velocity of the CO2 plume was the lowest in Case 8 (or the largest ), the CO2 plume swept the largest area, covering the entire thickness of the aquifer; the calculated mass center moved to the middle height of the aquifer (Figure 13(a)). Due to large coverage of the CO2 plume, the size of the corresponding As plume was the largest in Case 8 (Figure 13(d)). In contrast, Case 10 with high accelerated the advective groundwater flow and dispersion, resulting in the flattened shape of the CO2 plume sinking to the bottom of the aquifer (Figure 13(c)). Due to the poor sweeping efficiency of the CO2 plume, the As plume only migrated beneath the municipal well.

The As concentrations at the municipal well for Cases 8–10 are plotted in Figure 14(a). The arrival time of As concentration was the latest (4.95 years) in Case 8 due to a low . However, due to large coverage of the As plume, the As concentration at the municipal well continuously increased to 3.77 × 10−7 mol/L until 190 years. In contrast, for Case 10, the arrival time of the As plume was the shortest (0.95 years), reaching a peak of 1.73 × 10−7 mol/L. However, due to dominance of high advective flow, the As plume was flattened below the municipal well. Therefore, the amount of As captured from the municipal well decreased soon, reaching the background level (0.67 × 10−7 mol/L) only after 25 years. These results imply that the location of the wellbore (e.g., fully or partially penetrating well, the location of screen interval) and the size of the capture zone (e.g., the pumping capacity) are important characteristics for governing As concentration at the well.

The calculated health risk for each case is plotted in Figure 14(b). As expected, Case 10, which showed the smallest breakthrough of the As concentration (e.g., the smallest ), revealed the lowest risk. Both Cases 8 and 9 showed almost equivalent high-risk prediction even if the profile of the As concentration appeared differently at the municipal well (Case 8: ~80 years; Case 9: ~180 years). While accounting for the risk assessment, the exposure duration (ED) was chosen to be 30 years in this work [74]. Following (7), calculated from the As concentration profile is typically dependent on the duration of ED [31, 82, 83]. For example, can decrease when the breakthrough of As concentration is shorter than the ED. However, when the breakthrough of the As concentration is sufficiently longer than the chosen ED, does not change. Similarly, in both Cases 8 and 9, the profiles of As concentration at the municipal well were sufficiently longer than the chosen ED (Figure 14(a)). Due to this reason, for these two cases were similar to each other, and subsequently the predicted risk levels only showed a slight difference.

8. Conclusion

Even if there is no direct evidence showing that the stored CO2 has leaked to the shallow aquifers from any major CO2 injection demonstration sites [50, 84, 85], there exist a few natural analog sites indicating that naturally stored CO2 has leaked through preexisting fault systems [25, 8688]. These natural sites where CO2 leakage is primarily driven by geothermal or tectonic activities are typically less populated with lack of concern in potable groundwater resources. However, as the number of CO2 injection demonstration activities is growing, the secondary contamination of leaked CO2 to the shallow potable aquifer becomes an important issue. In this study, with the presence of arsenic-bearing minerals in the aquifer, it is suggested that water-rock interactions induced by CO2 leakage could mobilize arsenic species to the shallow potable aquifer. Therefore, without proper treatments, any residents who continuously utilize these groundwater resources could have great probability of exposure to cancer-related diseases.

Throughout this study, we developed a 2D confined aquifer model where gaseous CO2 was leaked, and a nearby municipal well concurrently pumped out both leaked CO2 and groundwater. Immediately after a leaked CO2 plume arrived at the municipal well, concentrations of As species increased together, indicating occurrence of active arsenopyrite dissolution in the aquifer media. Subsequently, by analyzing As concentration from the municipal well, the carcinogenic health risk was quantified. The processes characterizing the movement of gaseous CO2 and associated CO2-water-rock reactions were simulated with the multiphase, multispecies reactive transport model, and subsequent carcinogenic health risks were predicted with probabilistic approach. The simulation results revealed that the movement of leaked CO2 plume was governed by local flow fields within the shallow potable aquifer; three driving forces, that is, ambient groundwater flow, CO2 leakage-driven flow, and pumping-driven flow, were characterized. This complex flow field governed chemical reactions, resulting in the most drastic increase (0.77%) in permeability occurring at the rear margin of the CO2 plume where the vigorous mixing between the ambient groundwater and CO2-dissolved fluid accelerated dissolution of the minerals. Additionally, sensitivity studies were conducted while varying the sorption intensity, leakage rate of CO2, and horizontal permeability.

In summary, key factors that exacerbate the secondary contamination of arsenic species at the municipal well were physical characteristics of CO2 plume such as shape, size, and migration velocity of CO2 plume; these physical characteristics govern the area where arsenopyrite dissolves, subsequently affecting As concentration. Furthermore, the size of capture zone (e.g., screen interval, pump capacity) also governed the As concentration in the municipal well. Therefore, when the secondary contamination occurs within the shallow potable aquifer, it is suggested that the aquifer characteristics as well as the amount of leaked CO2 and its plume size should be evaluated to develop a proper remediation protocol. Additionally, the prompt shutdown of any nearby municipal wells could minimize any potential hazards.


A primitive version of this research was presented at the 2016 Fall General Assembly of the American Geophysical Union.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The research was partially supported by multiple programs including R&D Project on Environmental Management of Geologic CO2 Storage offered by the Korea Environmental Industry & Technology Institute (Project no. 201700180004), the R&D Convergence Program of National Research Council of Science & Technology in Korea (Grant no. CAP-15-07-KICT), and, finally, the Yonsei University Future-Leading Research Initiative of 2017-22-0044.