# CO2 Leakage-Induced Contamination in Shallow Potable Aquifer and Associated Health Risk Assessment

Authors: Chan Yeong Kim; Weon Shik Han; Eungyu Park; Jina Jeong; Tianfu Xu
Journal: Geofluids
 (2018)

---

## Abstract

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.


---

## Body

## 1. Introduction

Carbon capture and storage are considered to be one of the mitigating strategies for reducing CO2 emissions to the atmosphere [1–3]. 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 [4–6]. 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 [7–10]. 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 [13–16].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, 26–28]. Decreased pH, increased HCO3- 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 H3AsO40, H2AsO4-, HAsO42-, AsO43-, H3AsO30, and H2AsO3- 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) [36–38]. Under the reducing condition, uncharged As(III) species, such as H3AsO30, is dominant below pH 9.2 (Figure 1). However, under the oxidizing condition, H2AsO4- is dominant at pH < 6.9, while HAsO42- is dominant at higher pH (pH > 6.9). Moreover, H3AsO40 and AsO43- 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 H3AsO30 and H2AsO4- were dominant under reducing and oxidizing conditions, respectively [33, 39]. Here, the shallow aquifer was considered to be under the reducing condition, and thus, H3AsO30 was chosen as the most dominant species.Figure 1
Eh-pH diagram for arsenic species under H-O-As system at 1 bar and 25°C.


## 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.Figure 2
Flow chart explaining the sequence of both numerical simulation and probabilistic health risk assessment. Numerical simulation composed of multiphase and multispecies reactive transport modeling required both physical and chemical parameters within the shallow aquifer. Subsequently, the output of numerical simulation was utilized for calculating probabilistic health risk for carcinogenic effect.


## 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 (H3AsO30) as a primary species and incorporated associated aqueous complexes (e.g., H2AsO3-, H2AsO4-, H3AsO4(aq), HAsO2(aq), HAsO42-, 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]:(1)rn=±k25nuexp⁡-EanuR1T-1298.15+k25Hexp⁡-EaHR1T-1298.15aHnH+k25OHexp⁡-EaOHR1T-1298.15aOHnOHAn1-Ωnθη,where Ea (J/mol) is the activation energy and k25 (mol/m2/s) is the rate constant at 25°C with superscripts nu, H, and OH representing neutral, acid, and base mechanisms, respectively. R (J/mol/K) is the gas constant, T (K) is temperature, a is the activity of dissolved species under acid or base conditions, and n is the power term. An is the specific reactive surface area (cm2/g), and Ωn is the kinetic mineral saturation ratio of the nth 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 (kh=10-13 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.Figure 3
Conceptual model representing the shallow potable aquifer with a length and thickness of 200 m and 40 m, respectively. Blue arrows represent ambient groundwater flow with the rate of 2 cm/day. The grid blocks with the purple color stand for constant pressure boundaries; left and right boundary are 3.2 MPa and 2.8 MPa, respectively. Both upper and lower boundaries are assumed to be the no-flow boundary. CO2 leakage point (red) is located at 50 m away from the left boundary. In addition, the municipal well (green) located at 50 m away from the right boundary penetrates 30 m deep from the upper seal while the pumping activity was only conducted through the screened section.


### 5.2. Physical and Chemical Parameters

The aquifer was assumed to be homogeneous with its horizontal permeability (kh) and porosity (ϕ) to be 10–13 m2 and 0.2, respectively (Table 1); the anisotropy ratio (kv/kh) 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].Table 1
Input parameters for shallow groundwater aquifer materials. Relative permeability and capillary pressure follow equations developed by van Genuchten [40].
Aquifer characteristic  Permeability Horizontal 10−13 Vertical 10−14 Porosity 0.2 Density (kg/m3) 2600 Tortuosity 0.3 Compressibility 10−8 Diffusivity (m2/s) 10−9  Relative permeability  krl=S-1-1-S-1/mm2 S - = S l - S l r 1 - S l r - S g r krg=1-S-2(1-S-2) S l r: irreducible water saturation 0.20 S g r: irreducible gas saturation 0.05 m: exponent 0.457  Capillary pressure  pcap=-P0S∗-1/m-11-m S ∗ = S l - S l r 1 - S l r S l r: irreducible water saturation 0.20 m: exponent 0.457 P 0: strength coefficient 19.59 kPaThe 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, 56–58]. 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).Table 2
Mineral volume fractions for shallow groundwater aquifer materials. A total of 7 primary and 5 secondary minerals were selected.
Primary mineral Vol. fraction  Quartz SiO2 8.10 × 10 - 1 Kaolinite Al2Si2O5(OH)4 1.58 × 10 - 1 Illite K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2 1.28 × 10 - 2 K-feldspar KAlSi3O8 9.43 × 10 - 3 Arsenopyrite FeAsS 9.00 × 10 - 3 Calcite CaCO3 3.38 × 10 - 4 Magnesite MgCO3 2.12 × 10 - 4  Secondary mineral Vol. fraction  Chlorite Mg2.5Fe2.5Al2Si3O10(OH)8 0.00 Dolomite (CaMg(CO3)2) 0.00 Goethite FeO(OH) 0.00 Oligoclase CaNa4Al6Si14O40 0.00 Smectite-Ca Ca0.145Mg0.26Al1.77Si3.97O10(OH) 0.00Mineral 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, AlO2-, Ca2+, Cl−, Fe2+, H+, H2O, H3AsO3(aq), HCO3-, HSO42-, HS−, K+, Mg2+, Na+, O2(aq), SiO2(aq), and SO42- were selected as primary species, and associated aqueous complexes were chosen as secondary species (Table 3).Table 3
The 16 primary species and 39 secondary species are selected after the batch-reaction simulation.
Primary species (16)  A l O 2 - H+ H S O 4 2 - Na+ Ca2+ H2O HS- O2(aq) Cl− H3AsO3(aq) K+ SiO2(aq) Fe2+ H C O 3 - Mg2+ S O 4 2 -  Secondary species (39)  A l O H 2 + CO2(aq) H2S(aq) M g H C O 3 + Al(OH)3(aq) C O 3 2 - H3AsO4(aq) MgSO4(aq) Al3+ Fe3+ H 3 S i O 4 - NaCl(aq) AlOH2+ FeCl+ HAlO2(aq) N a C O 3 - CaCl2(aq) F e C l 4 2 - HAsO4(aq) NaHCO3(aq) CaCO3(aq) F e C O 3 + H A s O 4 2 - NaOH(aq) C a H C O 3 + F e H C O 3 + HAsS2(aq) N a S O 4 - CaOH+ H2(aq) HS− OH− CaSO4(aq) H 2 A s O 3 - H S O 3 - SO2(aq) CH4(aq) H 2 A s O 4 - MgCl+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 (Kd) for arsenic species under various conditions (e.g., different minerals and pH) and revealed that Kd 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 Kd for quartz and kaolinite was measured at 2 and 19 L/kg, respectively [39]. Therefore, the Kd 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].Table 4
Kinetic data for selected primary and secondary minerals. Only calcite was treated similar to the equilibrium reaction, and the rest of the minerals were kinetically treated with consideration of neutral, acid, and base mechanisms.
Mineral A n (cm2/g) Neutral mechanism Acid mechanism Base mechanism References k 25 (mol/m2·s) E a (kJ/mol) k 25 (mol/m2·s) E a (kJ/mol) n (H+) k 25 (mol/m2·s) E a (kJ/mol) n (H+) Quartz 9.8 1.02 × 10 - 14 87.7       [41] Kaolinite 1.95 × 10 5 6.91 × 10 - 14 22.2 4.89 × 10 - 12 65.9 0.78 8.91 × 10 - 18 17.9 −0.47 [42] Illite 6.68 × 10 5 1.66 × 10 - 13 35.0 1.05 × 10 - 14 23.6 0.34 3.02 × 10 - 17 85.9 −0.40 Set to smectites K-feldspar 9.8 3.89 × 10 - 13 38.0 8.71 × 10 - 11 51.7 0.50 6.31 × 10 - 12 94.1 −0.82 [42] Arsenopyrite 12.9 2.52 × 10 - 12 62.8       [43] Calcite Equilibrium reaction Magnesite 9.8 4.57 × 10 - 10 23.5 4.17 × 10 - 7 14.4 1.00    [42]  Chlorite 9.8 3.02 × 10 - 13 88.0 7.76 × 10 - 12 88.0 0.50    [42] Dolomite 12.9 2.52 × 10 - 12 52.2 2.34 × 10 - 7 43.5 0.50    [27] Goethite 9.8 2.52 × 10 - 12 62.8       [43] Oligoclase 9.8 1.44 × 10 - 12 69.8 2.13 × 10 - 10 65.0 0.46    [42] Smectite-Ca 5.64 × 10 5 1.66 × 10 - 13 35.0 1.05 × 10 - 14 12.6 0.34 3.02 × 10 - 17 58.9 −0.40 [42]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:(2)ϕ=1-∑m=1nmfrm-fru,where nm is the number of minerals and frm and fru are the volume fractions of mth mineral in the rock and nonreactive rock, respectively. Finally, permeability change was calculated using the porosity changes with cubic law [65]:(3)k=kiϕϕi3,where k and ϕ are permeability and porosity, respectively, with subscript i 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 linearKd approach instead of surface complexation. Although the linear Kd 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 [66–68]. 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, NO3-, and Fe3+. Under acidic conditions, Fe3+ quickly oxidizes, and sequential dissolution of arsenopyrite releases arsenic species (H3AsO40 and H3AsO30) 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].(4)FeAsS+1.50H2O+0.75O2aq+H+=Fe2++HS-+H3AsO3aq


### 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.(5)Risk=e1-ADDi×CPFi≈ADDi×CPFi.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:(6)ADD=C-INBWED×EFAT,where IN/BW 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 C-, 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, C-, calculated from concentration profiles (C(t)) of arsenic species at the municipal well in the numerical simulation can be represented as follows:(7)C-=max⁡∑tt+EDCt-∞+∞.Specific to ADD, exposure pathways through both ingestion and dermal sorption are defined as follows:(8)ADDingestion=C-IRBWED×EFAT,(9)ADDdermal=C-SABWED×EFATKpfskinEDshowerCF,where SA/BW is the skin surface area per body weight (m2/kg), Kp is the dermal permeability coefficient for arsenic in water, designated as 1.0 × 10−5 m/hour, fskin is the fraction of skin in contact with water (−), EDshower 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:(10)Total  Risk=Riskingestion+Riskdermal=ADDingestion×CPFingestion+ADDdermal×CPFdermal.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 C- (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].Table 5
Parameters for health risk assessment. Each parameter is randomly sampled within the designated value range and distributions to represent individual characteristics of potential victims.
Parameter Symbol Unit Distribution Values Exposure duration ED [yr] Constant 30 Exposure frequency EF [d/yr] Constant 350 Averaging time AT [d] Constant 25,550 Ingestion rate per unit body weight IR/BW [L·kg/d] Lognormal ( 3.3 × 10 - 2 , 1.3 × 10 - 2 ) Skin surface area per unit body weight SA/BW [m2/kg] Lognormal ( 2.7 × 10 - 2 , 2.5 × 10 - 3 ) Fraction skin in contact with water f s k i n [-] Uniform ( 4.0 × 10 - 1 , 9.0 × 10 - 1 ) Shower exposure duration E D s h o w e r [h/d] Lognormal ( 1.3 × 10 - 1 , 9.0 × 10 - 2 ) Unit conversion factor C F [L/m3] Constant 1.0 × 10 - 3 Cancer potency factor C P F i n g e s t i o n [kg·day/mg] Constant 1.5 C P F d e r m a l 1.6


## 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 (kh=10-13 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.Figure 3
Conceptual model representing the shallow potable aquifer with a length and thickness of 200 m and 40 m, respectively. Blue arrows represent ambient groundwater flow with the rate of 2 cm/day. The grid blocks with the purple color stand for constant pressure boundaries; left and right boundary are 3.2 MPa and 2.8 MPa, respectively. Both upper and lower boundaries are assumed to be the no-flow boundary. CO2 leakage point (red) is located at 50 m away from the left boundary. In addition, the municipal well (green) located at 50 m away from the right boundary penetrates 30 m deep from the upper seal while the pumping activity was only conducted through the screened section.


## 5.2. Physical and Chemical Parameters

The aquifer was assumed to be homogeneous with its horizontal permeability (kh) and porosity (ϕ) to be 10–13 m2 and 0.2, respectively (Table 1); the anisotropy ratio (kv/kh) 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].Table 1
Input parameters for shallow groundwater aquifer materials. Relative permeability and capillary pressure follow equations developed by van Genuchten [40].
Aquifer characteristic  Permeability Horizontal 10−13 Vertical 10−14 Porosity 0.2 Density (kg/m3) 2600 Tortuosity 0.3 Compressibility 10−8 Diffusivity (m2/s) 10−9  Relative permeability  krl=S-1-1-S-1/mm2 S - = S l - S l r 1 - S l r - S g r krg=1-S-2(1-S-2) S l r: irreducible water saturation 0.20 S g r: irreducible gas saturation 0.05 m: exponent 0.457  Capillary pressure  pcap=-P0S∗-1/m-11-m S ∗ = S l - S l r 1 - S l r S l r: irreducible water saturation 0.20 m: exponent 0.457 P 0: strength coefficient 19.59 kPaThe 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, 56–58]. 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).Table 2
Mineral volume fractions for shallow groundwater aquifer materials. A total of 7 primary and 5 secondary minerals were selected.
Primary mineral Vol. fraction  Quartz SiO2 8.10 × 10 - 1 Kaolinite Al2Si2O5(OH)4 1.58 × 10 - 1 Illite K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2 1.28 × 10 - 2 K-feldspar KAlSi3O8 9.43 × 10 - 3 Arsenopyrite FeAsS 9.00 × 10 - 3 Calcite CaCO3 3.38 × 10 - 4 Magnesite MgCO3 2.12 × 10 - 4  Secondary mineral Vol. fraction  Chlorite Mg2.5Fe2.5Al2Si3O10(OH)8 0.00 Dolomite (CaMg(CO3)2) 0.00 Goethite FeO(OH) 0.00 Oligoclase CaNa4Al6Si14O40 0.00 Smectite-Ca Ca0.145Mg0.26Al1.77Si3.97O10(OH) 0.00Mineral 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, AlO2-, Ca2+, Cl−, Fe2+, H+, H2O, H3AsO3(aq), HCO3-, HSO42-, HS−, K+, Mg2+, Na+, O2(aq), SiO2(aq), and SO42- were selected as primary species, and associated aqueous complexes were chosen as secondary species (Table 3).Table 3
The 16 primary species and 39 secondary species are selected after the batch-reaction simulation.
Primary species (16)  A l O 2 - H+ H S O 4 2 - Na+ Ca2+ H2O HS- O2(aq) Cl− H3AsO3(aq) K+ SiO2(aq) Fe2+ H C O 3 - Mg2+ S O 4 2 -  Secondary species (39)  A l O H 2 + CO2(aq) H2S(aq) M g H C O 3 + Al(OH)3(aq) C O 3 2 - H3AsO4(aq) MgSO4(aq) Al3+ Fe3+ H 3 S i O 4 - NaCl(aq) AlOH2+ FeCl+ HAlO2(aq) N a C O 3 - CaCl2(aq) F e C l 4 2 - HAsO4(aq) NaHCO3(aq) CaCO3(aq) F e C O 3 + H A s O 4 2 - NaOH(aq) C a H C O 3 + F e H C O 3 + HAsS2(aq) N a S O 4 - CaOH+ H2(aq) HS− OH− CaSO4(aq) H 2 A s O 3 - H S O 3 - SO2(aq) CH4(aq) H 2 A s O 4 - MgCl+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 (Kd) for arsenic species under various conditions (e.g., different minerals and pH) and revealed that Kd 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 Kd for quartz and kaolinite was measured at 2 and 19 L/kg, respectively [39]. Therefore, the Kd 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].Table 4
Kinetic data for selected primary and secondary minerals. Only calcite was treated similar to the equilibrium reaction, and the rest of the minerals were kinetically treated with consideration of neutral, acid, and base mechanisms.
Mineral A n (cm2/g) Neutral mechanism Acid mechanism Base mechanism References k 25 (mol/m2·s) E a (kJ/mol) k 25 (mol/m2·s) E a (kJ/mol) n (H+) k 25 (mol/m2·s) E a (kJ/mol) n (H+) Quartz 9.8 1.02 × 10 - 14 87.7       [41] Kaolinite 1.95 × 10 5 6.91 × 10 - 14 22.2 4.89 × 10 - 12 65.9 0.78 8.91 × 10 - 18 17.9 −0.47 [42] Illite 6.68 × 10 5 1.66 × 10 - 13 35.0 1.05 × 10 - 14 23.6 0.34 3.02 × 10 - 17 85.9 −0.40 Set to smectites K-feldspar 9.8 3.89 × 10 - 13 38.0 8.71 × 10 - 11 51.7 0.50 6.31 × 10 - 12 94.1 −0.82 [42] Arsenopyrite 12.9 2.52 × 10 - 12 62.8       [43] Calcite Equilibrium reaction Magnesite 9.8 4.57 × 10 - 10 23.5 4.17 × 10 - 7 14.4 1.00    [42]  Chlorite 9.8 3.02 × 10 - 13 88.0 7.76 × 10 - 12 88.0 0.50    [42] Dolomite 12.9 2.52 × 10 - 12 52.2 2.34 × 10 - 7 43.5 0.50    [27] Goethite 9.8 2.52 × 10 - 12 62.8       [43] Oligoclase 9.8 1.44 × 10 - 12 69.8 2.13 × 10 - 10 65.0 0.46    [42] Smectite-Ca 5.64 × 10 5 1.66 × 10 - 13 35.0 1.05 × 10 - 14 12.6 0.34 3.02 × 10 - 17 58.9 −0.40 [42]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:(2)ϕ=1-∑m=1nmfrm-fru,where nm is the number of minerals and frm and fru are the volume fractions of mth mineral in the rock and nonreactive rock, respectively. Finally, permeability change was calculated using the porosity changes with cubic law [65]:(3)k=kiϕϕi3,where k and ϕ are permeability and porosity, respectively, with subscript i 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 linearKd approach instead of surface complexation. Although the linear Kd 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 [66–68]. 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, NO3-, and Fe3+. Under acidic conditions, Fe3+ quickly oxidizes, and sequential dissolution of arsenopyrite releases arsenic species (H3AsO40 and H3AsO30) 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].(4)FeAsS+1.50H2O+0.75O2aq+H+=Fe2++HS-+H3AsO3aq


## 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.(5)Risk=e1-ADDi×CPFi≈ADDi×CPFi.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:(6)ADD=C-INBWED×EFAT,where IN/BW 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 C-, 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, C-, calculated from concentration profiles (C(t)) of arsenic species at the municipal well in the numerical simulation can be represented as follows:(7)C-=max⁡∑tt+EDCt-∞+∞.Specific to ADD, exposure pathways through both ingestion and dermal sorption are defined as follows:(8)ADDingestion=C-IRBWED×EFAT,(9)ADDdermal=C-SABWED×EFATKpfskinEDshowerCF,where SA/BW is the skin surface area per body weight (m2/kg), Kp is the dermal permeability coefficient for arsenic in water, designated as 1.0 × 10−5 m/hour, fskin is the fraction of skin in contact with water (−), EDshower 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:(10)Total  Risk=Riskingestion+Riskdermal=ADDingestion×CPFingestion+ADDdermal×CPFdermal.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 C- (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].Table 5
Parameters for health risk assessment. Each parameter is randomly sampled within the designated value range and distributions to represent individual characteristics of potential victims.
Parameter Symbol Unit Distribution Values Exposure duration ED [yr] Constant 30 Exposure frequency EF [d/yr] Constant 350 Averaging time AT [d] Constant 25,550 Ingestion rate per unit body weight IR/BW [L·kg/d] Lognormal ( 3.3 × 10 - 2 , 1.3 × 10 - 2 ) Skin surface area per unit body weight SA/BW [m2/kg] Lognormal ( 2.7 × 10 - 2 , 2.5 × 10 - 3 ) Fraction skin in contact with water f s k i n [-] Uniform ( 4.0 × 10 - 1 , 9.0 × 10 - 1 ) Shower exposure duration E D s h o w e r [h/d] Lognormal ( 1.3 × 10 - 1 , 9.0 × 10 - 2 ) Unit conversion factor C F [L/m3] Constant 1.0 × 10 - 3 Cancer potency factor C P F i n g e s t i o n [kg·day/mg] Constant 1.5 C P F d e r m a l 1.6


## 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 Kd, CO2 leakage rates, and horizontal permeability of the aquifer (kh). First, the degree of sorption intensity was evaluated by varying Kd 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, kh varied (kh = 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 kh accelerates the horizontal velocity of ambient groundwater while reducing buoyancy forces on the CO2 plume.Table 6
Different parameters for sensitivity studies. Three groups except Case 1 (base case) are designed for sorption rate (Kd) (Cases 2–4), gaseous CO2 leakage rate (QL) (Cases 5–7), and horizontal permeability (kh) (Cases 8–10), respectively.
Case name Sorption rate(Kd) CO2 leakage rate (QL) Horizontal permeability(kh) Case 1 (base case) 10 0.050 kg/s 1.0 × 10 - 13 m2  Case 2 25 0.050 kg/s 1.0 × 10 - 13 m2 Case 3 50 0.050 kg/s 1.0 × 10 - 13 m2 Case 4 100 0.050 kg/s 1.0 × 10 - 13 m2  Case 5 50 0.020 kg/s 1.0 × 10 - 13 m2 Case 6 50 0.025 kg/s 1.0 × 10 - 13 m2 Case 7 50 0.030 kg/s 1.0 × 10 - 13 m2  Case 8 50 0.050 kg/s 0.2 × 10 - 13 m2 Case 9 50 0.050 kg/s 1.0 × 10 - 13 m2 Case 10 50 0.050 kg/s 5.0 × 10 - 13 m2


## 7. Results and Discussion

### 7.1. Base Case

#### 7.1.1. Migration of Leaked CO2 Plume within the Shallow Potable Aquifer

Figures4(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)).Figure 4
Saturation of leaked CO2 at different times (120, 240, 360, 480, and 600 days). Red, black, and yellow circles with labels indicate the locations of mass center for CO2 plume. Red circles in (d) and (e) stand for the evolution of the mass center in multiple times; yellow circles in (e) express the mass center from 600 to 1,080 days.
(a)
120 days (b)
240 days (c)
360 days (d)
480 days (e)
600 daysAfter 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 HCO3- 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.Figure 5
Distribution of pH, carbonate minerals (calcite and magnesite), and associated dissolved species (HCO3-, Ca2+, and Mg2+) after 365 days.Dissolutions in calcite (CaCO3) and magnesite (MgCO3) were distinct (Figures 5(b) and 5(c)). Such dissolution released HCO3- 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 (SIcalcite=-3.1) was more saturated than magnesite (SImagnesite=-12.4); 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 HCO3- 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 (Figures6(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 AlO2- 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 AlO2- 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.Figure 6
Distribution of silicate minerals (K-feldspar, illite, chlorite, and kaolinite) and associated dissolved species (SiO2(aq), K+, and AlO2) after 365 days.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 AlO2-, 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% (k′%=(kchanged-kinitial)/kinitial×100) within the CO2 plume, and the most drastic increase (0.77%) occurred at the rear margin of the CO2 plume (Figure 7).Figure 7
Distribution of altered permeability (%) after 365 days. The changes in permeability were predicted from the following equation (k′%=(kchanged-kinitial)/kinitial×100) because its variation was too small. Here, kinitial and kchanged indicate the initial permeability and altered permeability at designated time, respectively.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 (H3AsO30(aq)), which is the by-product of arsenopyrite dissolution, as well as other arsenic species such as H2AsO3-, H3AsO4(aq), H2AsO4-, H2AsO42-, HAsO42-, 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 Kd approach. Kd 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−7mol/L and 3.6 × 10−7 mol/L, respectively (Figures 8(c) and 8(d)).Figure 8
Distribution of arsenopyrite dissolution and associated dissolved species (∑As, Fe2+, HS−, and O2(aq)) after 365 days.
(a)
Arsenopyrite (FeAsS) [mol/m3] (b)
∑As [mol/L] (c)
Fe2+ [mol/L] (d)
HS− [mol/L] (e)
O2(aq) [mol/L]


#### 7.1.3. Health Risk Assessment of Carcinogenic Effect

To account for carcinogenic health risk, selected species such as pH,∑As, and arsenite (H3AsO30(aq)) concentrations were monitored at the municipal well for 100 years (Figure 9(a)). Dramatic changes in gaseous CO2 saturation (Sg), pH, and mass fraction of CO2 dissolved in groundwater (XCO2) 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 (Sg, black dotted line) or dissolved CO2 (XCO2, 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 H3AsO30(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, Sg 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 (XCO2) 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 XCO2 became 0 and pH returned to 8.2. In addition, concentrations of ∑As and H3AsO30(aq) reached 3.77 × 10−7 and 1.70 × 10−7mol/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 H3AsO30(aq) were stabilized and gradually decreased due to sorption on aquifer media. The pH continuously decreased until 90 years. However, concentrations of ∑As and H3AsO30(aq) reached the background level (6.60 × 10−8 and 3.01 × 10−8 mol/L, resp.) at 62 years.Figure 9
(a) Changes in pH, total arsenic (∑As), and arsenite (H3AsO30) concentrations monitored at the municipal well. The maximum contaminant level (MCL) of arsenic species is 1.33 × 10−7 mol/L. Gaseous CO2 saturation (Sg) and pH for 10 years were also plotted at the magnified window. Two stages were characterized based on the behavior of CO2 plume. (b) Frequency of predicted health risk accounting for both ∑As and H3AsO30 and cumulative density functions representing probabilistic health risk. The bin size of both histograms was designated to be 10−4, and thus the number of bins was 10 and 20, respectively. Finally, the threshold of risk or the risk level of concern was 10−4 (red dotted line).
(a) (b)Based on the simulated profiles of both∑As and H3AsO30(aq), C-, 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 H3AsO30(aq) (green bar). The carcinogenic risk predicted from H3AsO30(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 H3AsO30(aq) exceeded the risk level (red dotted line) of concern.


### 7.2. Sensitivity Studies

#### 7.2.1. Effect of Sorption Intensity (Kd) (Cases 2–4)

In this sensitivity study, the intensity of the sorption effect (Kd 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 (Sg) were influenced by multiphase parameters such as capillary pressure and relative permeability (Table 1); the role of Kd 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 Kd, 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 Kd 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 Kd, 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 Kd, 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).Figure 10
(a)∑As profiles with different Kd of 25, 50, and 100 L/kg−1 (Cases 2–4) and (b) probabilistic health risk plotted as cumulative density functions. The red dotted line represents the risk level of concern (10−4).
(a) (b)In these simulations, with increasingKd, 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 (Kd = 25 L/kg), ∑As concentration revealed the earliest recovery (69.3 years) at the municipal well (Figure 10(a)). As Kd increased to 50 and 100 L/kg, the recovery time was delayed to 72.1 and 73.9 years, respectively.Figure10(b) represents the predicted probabilistic health risk for Cases 2–4. While calculating the health risk, C-, 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 C- 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 Kd. This result implies that variation in Kd was a less influential parameter for assessing health risk for arsenic species. This was because variation in Kd 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 (QL) (Cases 5–7)

In Cases 5–7, the effect of CO2 leakage rate (QL) 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 QL 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 QL 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 QL governed the size of the CO2 plume where the active dissolution of arsenopyrite, the source of ∑As, primarily occurred. Therefore, as QL 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.Figure 11
Saturation of leaked CO2 at different leakage rates after 1 year (365 days). (a) Case 5, (b) Case 6, and (c) Case 7. Black circles with labels indicate the locations of mass center for CO2 plume. Distribution of ∑As concentration induced by leaked CO2 at different leakage intensities (QL) after 20 years. (d) Case 5, (e) Case 6, and (f) Case 7.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 QL (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 QL (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 C-, which was the maximum average of arsenic concentration monitored at the municipal well (see (7)). The size of QL varied with C- 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).Figure 12
(a)∑As profiles with different leakage intensities (QL) of 0.01, 0.02, and 0.03 kg/s (Cases 5–7) and (b) probabilistic health risk for Cases 5–7 plotted as cumulative density functions (CDF). The red dotted line represents the risk level of concern (10−4).
(a) (b)


#### 7.2.3. Effect of Horizontal Permeability (kh) (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 (kh), 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 kh 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)).Figure 13
Saturation of leaked CO2 at different horizontal permeabilities (kh) after 1 year (365 days). (a) Case 8, (b) Case 9, and Case 10. Black circles with labels indicate the locations of mass center for CO2 plume. Distribution of ∑As Concentration induced by leaked CO2 at different horizontal permeabilities (kh) after 20 years: Cases 8, 9, and 10.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 (Nb=kvLΔρg/Hvμ) reflects the change in CO2 plume shape, where g is gravitational acceleration, L and H are the length and height of the model domain, respectively, v is the ambient groundwater flow rate, and μ is dynamic viscosity of CO2 (1.47 × 10−5Pa·s). Finally, Δρ is the density difference between groundwater (1,000 kg/m3) and CO2 (1.842 kg/m3) [81]. The calculated Nb 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 kh decreased.Even if the velocity of the CO2 plume was the lowest in Case 8 (or the largest Nb), 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 kh 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 kh. 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.Figure 14
(a)∑As profiles with different horizontal permeabilities (kh) of 0.2 × 10−13, 1.0× 10−13, and 5.0× 10–13 m2 (Cases 8–10) and (b) probabilistic health risk for Cases 8–10 plotted as cumulative density functions (CDF). The red dotted line represents the risk level of concern (10−4).
(a) (b)The calculated health risk for each case is plotted in Figure14(b). As expected, Case 10, which showed the smallest breakthrough of the ∑As concentration (e.g., the smallest C-), 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), C- calculated from the ∑As concentration profile is typically dependent on the duration of ED [31, 82, 83]. For example, C- 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, C- 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, C- for these two cases were similar to each other, and subsequently the predicted risk levels only showed a slight difference.


## 7.1. Base Case

### 7.1.1. Migration of Leaked CO2 Plume within the Shallow Potable Aquifer

Figures4(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)).Figure 4
Saturation of leaked CO2 at different times (120, 240, 360, 480, and 600 days). Red, black, and yellow circles with labels indicate the locations of mass center for CO2 plume. Red circles in (d) and (e) stand for the evolution of the mass center in multiple times; yellow circles in (e) express the mass center from 600 to 1,080 days.
(a)
120 days (b)
240 days (c)
360 days (d)
480 days (e)
600 daysAfter 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 HCO3- 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.Figure 5
Distribution of pH, carbonate minerals (calcite and magnesite), and associated dissolved species (HCO3-, Ca2+, and Mg2+) after 365 days.Dissolutions in calcite (CaCO3) and magnesite (MgCO3) were distinct (Figures 5(b) and 5(c)). Such dissolution released HCO3- 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 (SIcalcite=-3.1) was more saturated than magnesite (SImagnesite=-12.4); 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 HCO3- 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 (Figures6(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 AlO2- 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 AlO2- 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.Figure 6
Distribution of silicate minerals (K-feldspar, illite, chlorite, and kaolinite) and associated dissolved species (SiO2(aq), K+, and AlO2) after 365 days.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 AlO2-, 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% (k′%=(kchanged-kinitial)/kinitial×100) within the CO2 plume, and the most drastic increase (0.77%) occurred at the rear margin of the CO2 plume (Figure 7).Figure 7
Distribution of altered permeability (%) after 365 days. The changes in permeability were predicted from the following equation (k′%=(kchanged-kinitial)/kinitial×100) because its variation was too small. Here, kinitial and kchanged indicate the initial permeability and altered permeability at designated time, respectively.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 (H3AsO30(aq)), which is the by-product of arsenopyrite dissolution, as well as other arsenic species such as H2AsO3-, H3AsO4(aq), H2AsO4-, H2AsO42-, HAsO42-, 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 Kd approach. Kd 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−7mol/L and 3.6 × 10−7 mol/L, respectively (Figures 8(c) and 8(d)).Figure 8
Distribution of arsenopyrite dissolution and associated dissolved species (∑As, Fe2+, HS−, and O2(aq)) after 365 days.
(a)
Arsenopyrite (FeAsS) [mol/m3] (b)
∑As [mol/L] (c)
Fe2+ [mol/L] (d)
HS− [mol/L] (e)
O2(aq) [mol/L]


### 7.1.3. Health Risk Assessment of Carcinogenic Effect

To account for carcinogenic health risk, selected species such as pH,∑As, and arsenite (H3AsO30(aq)) concentrations were monitored at the municipal well for 100 years (Figure 9(a)). Dramatic changes in gaseous CO2 saturation (Sg), pH, and mass fraction of CO2 dissolved in groundwater (XCO2) 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 (Sg, black dotted line) or dissolved CO2 (XCO2, 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 H3AsO30(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, Sg 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 (XCO2) 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 XCO2 became 0 and pH returned to 8.2. In addition, concentrations of ∑As and H3AsO30(aq) reached 3.77 × 10−7 and 1.70 × 10−7mol/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 H3AsO30(aq) were stabilized and gradually decreased due to sorption on aquifer media. The pH continuously decreased until 90 years. However, concentrations of ∑As and H3AsO30(aq) reached the background level (6.60 × 10−8 and 3.01 × 10−8 mol/L, resp.) at 62 years.Figure 9
(a) Changes in pH, total arsenic (∑As), and arsenite (H3AsO30) concentrations monitored at the municipal well. The maximum contaminant level (MCL) of arsenic species is 1.33 × 10−7 mol/L. Gaseous CO2 saturation (Sg) and pH for 10 years were also plotted at the magnified window. Two stages were characterized based on the behavior of CO2 plume. (b) Frequency of predicted health risk accounting for both ∑As and H3AsO30 and cumulative density functions representing probabilistic health risk. The bin size of both histograms was designated to be 10−4, and thus the number of bins was 10 and 20, respectively. Finally, the threshold of risk or the risk level of concern was 10−4 (red dotted line).
(a) (b)Based on the simulated profiles of both∑As and H3AsO30(aq), C-, 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 H3AsO30(aq) (green bar). The carcinogenic risk predicted from H3AsO30(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 H3AsO30(aq) exceeded the risk level (red dotted line) of concern.


## 7.1.1. Migration of Leaked CO2 Plume within the Shallow Potable Aquifer

Figures4(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)).Figure 4
Saturation of leaked CO2 at different times (120, 240, 360, 480, and 600 days). Red, black, and yellow circles with labels indicate the locations of mass center for CO2 plume. Red circles in (d) and (e) stand for the evolution of the mass center in multiple times; yellow circles in (e) express the mass center from 600 to 1,080 days.
(a)
120 days (b)
240 days (c)
360 days (d)
480 days (e)
600 daysAfter 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 HCO3- 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.Figure 5
Distribution of pH, carbonate minerals (calcite and magnesite), and associated dissolved species (HCO3-, Ca2+, and Mg2+) after 365 days.Dissolutions in calcite (CaCO3) and magnesite (MgCO3) were distinct (Figures 5(b) and 5(c)). Such dissolution released HCO3- 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 (SIcalcite=-3.1) was more saturated than magnesite (SImagnesite=-12.4); 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 HCO3- 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 (Figures6(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 AlO2- 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 AlO2- 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.Figure 6
Distribution of silicate minerals (K-feldspar, illite, chlorite, and kaolinite) and associated dissolved species (SiO2(aq), K+, and AlO2) after 365 days.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 AlO2-, 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% (k′%=(kchanged-kinitial)/kinitial×100) within the CO2 plume, and the most drastic increase (0.77%) occurred at the rear margin of the CO2 plume (Figure 7).Figure 7
Distribution of altered permeability (%) after 365 days. The changes in permeability were predicted from the following equation (k′%=(kchanged-kinitial)/kinitial×100) because its variation was too small. Here, kinitial and kchanged indicate the initial permeability and altered permeability at designated time, respectively.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 (H3AsO30(aq)), which is the by-product of arsenopyrite dissolution, as well as other arsenic species such as H2AsO3-, H3AsO4(aq), H2AsO4-, H2AsO42-, HAsO42-, 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 Kd approach. Kd 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−7mol/L and 3.6 × 10−7 mol/L, respectively (Figures 8(c) and 8(d)).Figure 8
Distribution of arsenopyrite dissolution and associated dissolved species (∑As, Fe2+, HS−, and O2(aq)) after 365 days.
(a)
Arsenopyrite (FeAsS) [mol/m3] (b)
∑As [mol/L] (c)
Fe2+ [mol/L] (d)
HS− [mol/L] (e)
O2(aq) [mol/L]


## 7.1.3. Health Risk Assessment of Carcinogenic Effect

To account for carcinogenic health risk, selected species such as pH,∑As, and arsenite (H3AsO30(aq)) concentrations were monitored at the municipal well for 100 years (Figure 9(a)). Dramatic changes in gaseous CO2 saturation (Sg), pH, and mass fraction of CO2 dissolved in groundwater (XCO2) 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 (Sg, black dotted line) or dissolved CO2 (XCO2, 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 H3AsO30(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, Sg 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 (XCO2) 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 XCO2 became 0 and pH returned to 8.2. In addition, concentrations of ∑As and H3AsO30(aq) reached 3.77 × 10−7 and 1.70 × 10−7mol/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 H3AsO30(aq) were stabilized and gradually decreased due to sorption on aquifer media. The pH continuously decreased until 90 years. However, concentrations of ∑As and H3AsO30(aq) reached the background level (6.60 × 10−8 and 3.01 × 10−8 mol/L, resp.) at 62 years.Figure 9
(a) Changes in pH, total arsenic (∑As), and arsenite (H3AsO30) concentrations monitored at the municipal well. The maximum contaminant level (MCL) of arsenic species is 1.33 × 10−7 mol/L. Gaseous CO2 saturation (Sg) and pH for 10 years were also plotted at the magnified window. Two stages were characterized based on the behavior of CO2 plume. (b) Frequency of predicted health risk accounting for both ∑As and H3AsO30 and cumulative density functions representing probabilistic health risk. The bin size of both histograms was designated to be 10−4, and thus the number of bins was 10 and 20, respectively. Finally, the threshold of risk or the risk level of concern was 10−4 (red dotted line).
(a) (b)Based on the simulated profiles of both∑As and H3AsO30(aq), C-, 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 H3AsO30(aq) (green bar). The carcinogenic risk predicted from H3AsO30(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 H3AsO30(aq) exceeded the risk level (red dotted line) of concern.


## 7.2. Sensitivity Studies

### 7.2.1. Effect of Sorption Intensity (Kd) (Cases 2–4)

In this sensitivity study, the intensity of the sorption effect (Kd 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 (Sg) were influenced by multiphase parameters such as capillary pressure and relative permeability (Table 1); the role of Kd 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 Kd, 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 Kd 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 Kd, 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 Kd, 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).Figure 10
(a)∑As profiles with different Kd of 25, 50, and 100 L/kg−1 (Cases 2–4) and (b) probabilistic health risk plotted as cumulative density functions. The red dotted line represents the risk level of concern (10−4).
(a) (b)In these simulations, with increasingKd, 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 (Kd = 25 L/kg), ∑As concentration revealed the earliest recovery (69.3 years) at the municipal well (Figure 10(a)). As Kd increased to 50 and 100 L/kg, the recovery time was delayed to 72.1 and 73.9 years, respectively.Figure10(b) represents the predicted probabilistic health risk for Cases 2–4. While calculating the health risk, C-, 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 C- 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 Kd. This result implies that variation in Kd was a less influential parameter for assessing health risk for arsenic species. This was because variation in Kd 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 (QL) (Cases 5–7)

In Cases 5–7, the effect of CO2 leakage rate (QL) 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 QL 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 QL 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 QL governed the size of the CO2 plume where the active dissolution of arsenopyrite, the source of ∑As, primarily occurred. Therefore, as QL 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.Figure 11
Saturation of leaked CO2 at different leakage rates after 1 year (365 days). (a) Case 5, (b) Case 6, and (c) Case 7. Black circles with labels indicate the locations of mass center for CO2 plume. Distribution of ∑As concentration induced by leaked CO2 at different leakage intensities (QL) after 20 years. (d) Case 5, (e) Case 6, and (f) Case 7.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 QL (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 QL (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 C-, which was the maximum average of arsenic concentration monitored at the municipal well (see (7)). The size of QL varied with C- 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).Figure 12
(a)∑As profiles with different leakage intensities (QL) of 0.01, 0.02, and 0.03 kg/s (Cases 5–7) and (b) probabilistic health risk for Cases 5–7 plotted as cumulative density functions (CDF). The red dotted line represents the risk level of concern (10−4).
(a) (b)


### 7.2.3. Effect of Horizontal Permeability (kh) (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 (kh), 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 kh 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)).Figure 13
Saturation of leaked CO2 at different horizontal permeabilities (kh) after 1 year (365 days). (a) Case 8, (b) Case 9, and Case 10. Black circles with labels indicate the locations of mass center for CO2 plume. Distribution of ∑As Concentration induced by leaked CO2 at different horizontal permeabilities (kh) after 20 years: Cases 8, 9, and 10.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 (Nb=kvLΔρg/Hvμ) reflects the change in CO2 plume shape, where g is gravitational acceleration, L and H are the length and height of the model domain, respectively, v is the ambient groundwater flow rate, and μ is dynamic viscosity of CO2 (1.47 × 10−5Pa·s). Finally, Δρ is the density difference between groundwater (1,000 kg/m3) and CO2 (1.842 kg/m3) [81]. The calculated Nb 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 kh decreased.Even if the velocity of the CO2 plume was the lowest in Case 8 (or the largest Nb), 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 kh 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 kh. 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.Figure 14
(a)∑As profiles with different horizontal permeabilities (kh) of 0.2 × 10−13, 1.0× 10−13, and 5.0× 10–13 m2 (Cases 8–10) and (b) probabilistic health risk for Cases 8–10 plotted as cumulative density functions (CDF). The red dotted line represents the risk level of concern (10−4).
(a) (b)The calculated health risk for each case is plotted in Figure14(b). As expected, Case 10, which showed the smallest breakthrough of the ∑As concentration (e.g., the smallest C-), 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), C- calculated from the ∑As concentration profile is typically dependent on the duration of ED [31, 82, 83]. For example, C- 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, C- 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, C- for these two cases were similar to each other, and subsequently the predicted risk levels only showed a slight difference.


## 7.2.1. Effect of Sorption Intensity (Kd) (Cases 2–4)

In this sensitivity study, the intensity of the sorption effect (Kd 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 (Sg) were influenced by multiphase parameters such as capillary pressure and relative permeability (Table 1); the role of Kd 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 Kd, 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 Kd 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 Kd, 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 Kd, 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).Figure 10
(a)∑As profiles with different Kd of 25, 50, and 100 L/kg−1 (Cases 2–4) and (b) probabilistic health risk plotted as cumulative density functions. The red dotted line represents the risk level of concern (10−4).
(a) (b)In these simulations, with increasingKd, 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 (Kd = 25 L/kg), ∑As concentration revealed the earliest recovery (69.3 years) at the municipal well (Figure 10(a)). As Kd increased to 50 and 100 L/kg, the recovery time was delayed to 72.1 and 73.9 years, respectively.Figure10(b) represents the predicted probabilistic health risk for Cases 2–4. While calculating the health risk, C-, 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 C- 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 Kd. This result implies that variation in Kd was a less influential parameter for assessing health risk for arsenic species. This was because variation in Kd 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 (QL) (Cases 5–7)

In Cases 5–7, the effect of CO2 leakage rate (QL) 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 QL 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 QL 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 QL governed the size of the CO2 plume where the active dissolution of arsenopyrite, the source of ∑As, primarily occurred. Therefore, as QL 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.Figure 11
Saturation of leaked CO2 at different leakage rates after 1 year (365 days). (a) Case 5, (b) Case 6, and (c) Case 7. Black circles with labels indicate the locations of mass center for CO2 plume. Distribution of ∑As concentration induced by leaked CO2 at different leakage intensities (QL) after 20 years. (d) Case 5, (e) Case 6, and (f) Case 7.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 QL (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 QL (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 C-, which was the maximum average of arsenic concentration monitored at the municipal well (see (7)). The size of QL varied with C- 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).Figure 12
(a)∑As profiles with different leakage intensities (QL) of 0.01, 0.02, and 0.03 kg/s (Cases 5–7) and (b) probabilistic health risk for Cases 5–7 plotted as cumulative density functions (CDF). The red dotted line represents the risk level of concern (10−4).
(a) (b)


## 7.2.3. Effect of Horizontal Permeability (kh) (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 (kh), 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 kh 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)).Figure 13
Saturation of leaked CO2 at different horizontal permeabilities (kh) after 1 year (365 days). (a) Case 8, (b) Case 9, and Case 10. Black circles with labels indicate the locations of mass center for CO2 plume. Distribution of ∑As Concentration induced by leaked CO2 at different horizontal permeabilities (kh) after 20 years: Cases 8, 9, and 10.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 (Nb=kvLΔρg/Hvμ) reflects the change in CO2 plume shape, where g is gravitational acceleration, L and H are the length and height of the model domain, respectively, v is the ambient groundwater flow rate, and μ is dynamic viscosity of CO2 (1.47 × 10−5Pa·s). Finally, Δρ is the density difference between groundwater (1,000 kg/m3) and CO2 (1.842 kg/m3) [81]. The calculated Nb 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 kh decreased.Even if the velocity of the CO2 plume was the lowest in Case 8 (or the largest Nb), 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 kh 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 kh. 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.Figure 14
(a)∑As profiles with different horizontal permeabilities (kh) of 0.2 × 10−13, 1.0× 10−13, and 5.0× 10–13 m2 (Cases 8–10) and (b) probabilistic health risk for Cases 8–10 plotted as cumulative density functions (CDF). The red dotted line represents the risk level of concern (10−4).
(a) (b)The calculated health risk for each case is plotted in Figure14(b). As expected, Case 10, which showed the smallest breakthrough of the ∑As concentration (e.g., the smallest C-), 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), C- calculated from the ∑As concentration profile is typically dependent on the duration of ED [31, 82, 83]. For example, C- 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, C- 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, C- 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, 86–88]. 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.



---

*Source: 4834601-2018-04-05.xml*