Quantification of Carbon Cycling in a Large Aquifer Using Reactive Transport Modelling

Miotliński, Konrad; Peeters, Luk J. M.

doi:10.3389/frwa.2021.714075

ORIGINAL RESEARCH article

Front. Water, 10 September 2021

Sec. Water and Hydrocomplexity

Volume 3 - 2021 | https://doi.org/10.3389/frwa.2021.714075

This article is part of the Research Topic Modeling-based Approaches for Water Resources Problems View all 7 articles

Quantification of Carbon Cycling in a Large Aquifer Using Reactive Transport Modelling

$\nKonrad Miotli&#x;ski$ Konrad Miotliński¹^*

Luk J. M. Peeters²

¹Centre for People, Place and Planet, Edith Cowan University, Joondalup, WA, Australia
²CSIRO Land and Water, Glen Osmond, SA, Australia

Continental scale aquifers can store significant amounts of carbon as a result of immense water volumes, substantial concentrations of dissolved inorganic carbon (DIC) and its reactions with a matrix, thus contributing the global carbon storage and cycle. However, concentration of dissolved solutes may vary significantly over distances, which causes interpretative challenges and difficulties in process quantification. This occurs in the Guarani Aquifer System in South America, which is a subject of extensive research due to a significant strategic role in water supply. Dissolved CO₂ is expected to dissociate and undergo reactions with aluminosilicate minerals, but it is unknown how much DIC may get immobilised in the aquifer. To quantify the processes, we performed reactive transport modelling which combines hydrological and geochemical information followed by global sensitivity analysis. We show that more than a half of the infiltrated CO₂ may be consistently precipitated as CaCO₃. The DIC concentrations across the aquifer depend primarily on the input carbon concentrations and the plagioclase hydrolysis rate, while other parameters including hydraulic conductivity, recharge rate and mineral stability are of the minor importance. We present how advanced modelling techniques may be used to interpret and quantify processes controlling water quality in continental scale groundwater systems.

1. Introduction

Groundwater systems are important inorganic carbon reservoirs, and they are gaining attention as to which extent they contribute to the global carbon cycle. Estimates of carbon fluxes in the subsurface are hugely uncertain due to difficulties in (1) calculating groundwater and chemical fluxes (Fontaine et al., 2007; McCallum et al., 2012; Li et al., 2015), (2) subsurface conceptualisation (Vilhelmsen et al., 2018; Enemark et al., 2019), and (3) the dynamic role of soils in CO₂ storage and deep recharge (Kessler and Harvey, 2001; Fontaine et al., 2007; Sánchez-Cañete et al., 2018).

Reactive transport models (RTMs) are useful in addressing these uncertainties because they (1) combine complexity of water flow with a set of biogeochemical reactions (Maher and Navarre-Sitchler, 2019), (2) enable to test a variety of conceptual scenarios when they are run in a predictive mode (Jessen et al., 2017; Jakobsen et al., 2018), (3) may be evaluated using global sensitivity analysis (Dai et al., 2014). Consequently, RTMs can play an important role at every stage of a research project: initially, when little data is available to test a general conceptual model; at the stage of data collection, when combination with global sensitivity analysis offers an efficient way to allocate financial resources or when new conceptual models are available; and finally, to integrate a range of datasets to draw quantitative conclusions and to make predictions.

The data integration remains a major challenge in water chemistry and quality evaluation, which may be independent of the scale of evaluation. Nevertheless, as the modelling techniques improve, new hypotheses may be drawn, particularly when combined with new observations. The DIC reaction with aluminosilicate minerals, which is referred as chemical rock weathering may be associated with precipitation of secondary minerals resulting in a long-term carbon storage (Berner, 1998; Banks and Frengstad, 2006; Steffen et al., 2007; Maher and Chamberlain, 2014). Nevertheless, the reaction pathways are complex and heterogeneous and may depend on a variety of factors including: (1) groundwater flow velocity (Maher and Chamberlain, 2014), (2) thermodynamic equilibria and dissolved ion availability (Zhang and Planavsky, 2020) or (3) rates of dissolution and precipitation (Pogge von Strandmann et al., 2019). These pathways can be explored with RTMs that rely on both hydrological and biogeochemical principles (Steefel et al., 2005).

In a small scale of soil thickness, in situ measurements of a CO₂ and O₂ couple indicate that CO₂ is produced in a shallow subsurface which may not, at least immediately, return to the atmosphere due to (1) dissolution in groundwater, (2) transport in gaseous or dissolved forms, (3) mineral weathering (Sánchez-Cañete et al., 2018). The dynamic interactions at a scale of 10⁰–10¹ m influence what happens along the way of water flow. Consequently, the CO₂ dissolution coupled with advective transport may lead to significant amounts of dissolved inorganic carbon (DIC) stored in groundwater bodies (Li et al., 2015; Zhang and Planavsky, 2020) that may be transported over long distances, undergo biogeochemical transformations and, finally, be discharged to rivers (Gaillardet et al., 1999) or sea (Zhang and Planavsky, 2020).

In aquifers, the process dynamics are subdued. However, water quality patterns observed along flow paths may give an insight of reaction controls. The spatial patterns have been quantified at a scale from 10¹–10² m Jakobsen et al. (2018) to 10⁵–10⁶ m (Appelo, 1994; Bethke et al., 2002). The large-scale example includes one of the world's most voluminous Guarani Aquifer System (GAS). The spatial patterns are the result of reactions within the aquifer (Sracek and Hirata, 2002; Gastmans et al., 2010) and mixing with other quality groundwater (Aggarwal et al., 2015; Fernandes et al., 2016; Teramoto et al., 2020).

To estimate the carbon fluxes in GAS, we developed a coupled RTM in which we integrate physical properties with the thermodynamic reactions assuming DIC reacts with the abundant aluminosilicate minerals. We verified our model by comparing its results with the groundwater quality data observed along the flow paths in the northeastern part of the aquifer. Independent studies imply dissolution of primary silicates and precipitation of carbonates. Finally, we performed global sensitivity analysis to identify the importance of individual model parameters.

2. Site Description and Methodology

Hydrogeology and geochemical evolution of the north-eastern part of GAS are described in detail elsewhere. Here only a brief description is provided. This is followed by the conceptual and numerical modelling.

2.1. Hydrogeology

The Guarani Aquifer System (GAS) is one of the most voluminous groundwater bodies on the Earth due to its great geographical extend (area of 1,195,500 km²) and a large thickness (Foster et al., 2009). In the NE part, the aquifer is a part of the Paraná Basin which consists of stratified sandstone series of Triassic (Piramboia unit) and Cretaceous (Botucatu unit) ages (Araujo et al., 1999; Gastmans et al., 2010; Hirata et al., 2011). The formation is up to 832 m thick, with the average thickness of 400 m (Araujo et al., 1999). Reported hydraulic conductivity and porosity values are in the range 0.0864–0.864 m/d (De Paula E Silva and Cavaguti, 1994; Rebouças, 1994) and 0.16–0.24 (Rebouças, 1976; Hirata et al., 2011), respectively. About 10% of the formation outcrops, whereas the remaining part is confined by a sequence of basalts of the Serra Geral formation (Fernandes et al., 2016). The GAS is underlain by siltstones containing highly mineralised groundwater of the Passa Dois Group (Silva, 1983; Meng and Maynard, 2001).

Natural aquifer recharge of GAS dominates as direct percolation of precipitation through the outcrop areas. The deep groundwater recharge estimated with the water level fluctuation method accounts locally for up to 50 mm/year, which is 3.5% of the mean annual precipitation (Wendland et al., 2007). No estimates have been performed for the confined part of the aquifer, but since the basalt formation is intersected with dikes (Fernandes et al., 2016) and groundwater samples contain ¹⁴C (Aggarwal et al., 2015), infiltration of recent rainwater through the overlying basalt unit seems to be taking place. Downstream, in the deeply confined part of the aquifer, elevated Cl⁻ concentrations along with crust ⁴He gas dissolved in groundwater suggest recharge from the underlying formations (Aggarwal et al., 2015). The groundwater flow is radially diverging and the discharge takes place in the Paraná River valley, where artesian conditions are apparent (Rebouças, 1994).

Groundwater residence times evaluated using the ⁸¹Kr isotope vary in the broad range of 80–820 ky (Aggarwal et al., 2015). The groundwater dating implies that the average linear groundwater velocities are from 0.15 to 1.25 m/year and, consequently, that GAS has been flushed at least 180 times since its deposition which was followed by a number of tectonic events with the most recent one between 88 and 65 million years ago (Araujo et al., 1999).

2.2. Spatial Water Quality Patterns

The chemical composition of groundwater exhibits a zonation from fresh Ca-HCO₃ and Ca-Mg-HCO₃ water types in the unconfined and marginal part of the basin towards brackish Na-HCO₃-Cl and Na-Cl in the confined part (Gastmans et al., 2010). The evolution has been suggested to be a result of either aluminosilicate weathering (Silva, 1983; Gastmans et al., 2010) or cation exchange coupled to carbonate dissolution (Sracek and Hirata, 2002). There is more evidence for the former mechanism, which includes: (1) multiple flushing since the formation of the aquifer (Araujo et al., 1999); (2) the presence of primary (plagioclase, K-feldspar) and secondary aluminosilicates (kaolinite, smectites) (Gesicki, 2007); (3) the presence of chalcedony near the outcrop and authigenic calcite cement in the deep confining zone (Hirata et al., 2011). A similar trend was observed and modelled elsewhere (Banks and Frengstad, 2006; Pogge von Strandmann et al., 2019).

2.3. Conceptual and Numerical Models

We conceptualised a groundwater system as a two-dimensional (2D) vertical cross-section parallel to the direction of groundwater flow (Freeze and Witherspoon, 1968). The system is a rectangular with the dimensions 300 km and 400 m in length and thickness, respectively (Figure 1). We imposed groundwater flow field by using the flux boundary condition along the upper and lower boundaries of the system.

FIGURE 1

Figure 1. Conceptual 2D model of the water flow in the GAS with shaded polygons depicting water fluxes through the top and bottom boundaries.

The USGS PHAST code (Parkhurst et al., 2010) and PHREEQC database (Parkhurst and Appelo, 1999) were used to simulate reactive transport in the 2D domain. The methodology was used in smaller scale studies (Jessen et al., 2017; Jakobsen et al., 2018). Here, the domain was divided into 31 and 401 nodes in a vertical and horizontal direction, respectively, which correspond to the depth and distance. In addition to flux boundary conditions at the upper boundary, we assigned a single node at 300 km with a constant boundary condition to assure a numerical stability of the solution (Figure 1).

Chemical fluxes through the system are products of water flux and concentrations. The recharge water at the outcrop contains dissolved CO₂, which is meant to react with the aquifer material. The basement water is highly-mineralised Na-Cl groundwater of the Passa Dois Group which discharges Na⁺ and Cl⁻ to GAS (Meng and Maynard, 2001; Hirata et al., 2011; Soler i Gil and Bonotto, 2014; Aggarwal et al., 2015).

The aquifer matrix consists of plagioclase (Ca_0.5Na_0.5Al_1.5Si_2.5O₈) which is allowed to develop: (1) chalcedony (SiO₂), (2) kaolinite (Al₂Si₂O₅(OH)₄) and calcite (CaCO₃) precipitation if mineral saturation is reached. The plagioclase hydrolysis follows the kinetic Transition State Theory (Appelo and Postma, 2005; Maher and Navarre-Sitchler, 2019):

\begin{array}{l} R = r (1 - \frac{I A P}{K}) \frac{A_{0}}{V} {(\frac{m}{m_{0}})}^{n} & (1) \end{array}

where R is the overall reaction rate (M/L/s), r is the rate constant (M/m²/s), IAP is the ion activity product (-), K is the solubility product (-), A₀ is the initial surface area of the solid (m²), V is the volume of solution (m³), m₀ initial moles of the solid, m moles of solid at a given time, n rate exponent. We assumed that m₀ and m are equal 10 M, A₀ = 2.26 m²/kg and V = 0.162 L/kg (Appelo and Postma, 2005). This is a conservative approach if concentrations achieve steady state. The m₀, A₀, and the IAP change over the course of a simulation making the overall reaction rate variable in time and space.

The IAP/K ratio controls the mineral saturation. The mineral saturation index is:

\begin{array}{l} S I_{m i n e r a l} = l o g (\frac{I A P}{K}) & (2) \end{array}

To ensure that the geochemical system reaches steady state in terms of water quality, each model run was 1 million years. All model parameters are listed in Table 1). A few parameters of the RTM were fixed in the model: (1) the pH value at of recharge water (4.4) as mean value observed in the outcrop zone (Gastmans et al., 2010), (2) the maximum recharge rate at the bottom (0.00015 m/yr) and (3) the concentrations at the bottom recharge (5 mM Na⁺ and Cl⁻, respectively).

TABLE 1

Table 1. Parametrization of the reactive transport models with ranges explored in global sensitivity analysis.

An ensemble of 1,000 parameter combinations was generated with Latin Hypercube Sampling and evaluated with the RTM to carry out a global sensitivity analysis. Table 1 includes the mean (X), standard deviation (gamma) and additional comments for each model parameter. We used the delta moment-independent measure technique developed by Borgonovo (2007) and Plischke et al. (2013), as implemented in the SALib python library (Herman and Usher, 2017). The technique calculates a total sensitivity index value between a parameter of the model and a model output. Higher values indicate that the model outcome is more sensitive to parameter values.

2.4. Water Quality Data

To compare the modelling results with the field data (Gastmans et al., 2010) we selected a number of monitoring wells that are located along interpreted flow paths and we estimated distances from the basin's margin (Table 2 in Appendix).

3. Results and Discussion

3.1. Hydrogeochemical Evolution in a 2D Siliciclastic System

The flux of water with CO₂ solution together with thermodynamic feasibility for secondary reactions create (1) distinct peaks of Ca²⁺ and total alkalinity ( $H C O_{3}^{-} + C O_{3}^{2 -} + H_{3} S i O_{4}^{-} + H_{2} S i O_{4}^{2 -}$ ) and (2) a gradual increase in pH (Figure 2). The pattern on a smaller spatial scale was identified in CO₂ driven silicate weathering systems. That included plagioclase-rich anorthosites in Norway (Banks and Frengstad, 2006) naturally recharged with CO₂-rich water, and in engineered CO₂ injection into basalts in Iceland (Pogge von Strandmann et al., 2019).

FIGURE 2

Figure 2. Profile average (measured and modelled) and 2D distribution of concentrations (mM) and pH due to plagioclase hydrolysis and associated secondary reactions. Input parameters for the model are: Rech_0 = 0.0025, Rech_50 = 0.0002, Kx = 0.33, Kz = 0.033, pCO2 = –1.5, log r_Plag = –17.22, log K_Cal = –8.48, n = 0.22, CEC = 2.5. Different symbols for measured data indicate separate flow paths.

Plagioclase hydrolysis drives the stoichiometric increases in Ca²⁺ and Na⁺ and the intensity of the process depends on the chemical rate (R in Equation 1). At low concentration levels, the process consumes dissolved CO₂ and hydrogen ions producing silicic acid and carbonate alkalinity (Equation 3):

\begin{array}{l} 2 C a_{0.5} N a_{0.5} A l_{1.5} S i_{2.5} O_{8} + C O_{2} + 6.5 H_{2} O + 2 H^{+} \to \\ 1.5 A l_{2} S i_{2} O_{5} {(O H)}_{4} + 2 H_{4} S i O_{4}^{0} + N a^{+} + C a^{2 +} + H C O_{3}^{-} & (3) \end{array}

Once groundwater becomes saturated with respect to chalcedony, its precipitation takes place which prevents an ongoing increase of SiO₂ in the solution (Equation 4):

\begin{array}{l} 2 C a_{0.5} N a_{0.5} A l_{1.5} S i_{2.5} O_{8} + C O_{2} + 2.5 H_{2} O + H^{+} \to \\ 1.5 A l_{2} S i_{2} O_{5} {(O H)}_{4} + 2 S i O_{2} + N a^{+} + C a^{2 +} + H C O_{3}^{-} & (4) \end{array}

As the evolution continues, at a pH ≈ 8 (x=140 km) groundwater becomes saturated with respect to calcite. At this pH value, the $HC O_{3}^{-}$ is a dominant dissolve carbonate form which is being consumed, consequently, leading to calcite precipitation (Equation 5):

\begin{array}{l} 2 C a_{0.5} N a_{0.5} A l_{1.5} S i_{2.5} O_{8} + H C O_{3}^{-} + 2.5 H_{2} O + H^{+} \to \\ 1.5 A l_{2} S i_{2} {(O H)}_{4} + 2 S i O_{2} + N a^{+} + C a C O_{3} & (5) \end{array}

The concentrations of carbonate species and calcite saturation index (Equation 2) and precipitation are depicted in Figure 3. Concentrations of Ca², carbonate alkalinity and pH observed in the field enable the comparison with the model. Banks and Frengstad (2006) found the same.

FIGURE 3

Figure 3. Modelled concentrations of dissolved carbonate species (upper), saturation index (SI) for calcite (line, middle), and rate of calcite accumulation in a function of distance (lower). SI calculated with PHREEQC (Parkhurst and Appelo, 1999) from the observations is also shown (points, middle). Input parameters for the model are: Rech_0 = 0.0025, Rech_50 = 0.0002, Kx = 0.33, Kz = 0.033, pCO2 = –1.5, log r_Plag = –17.22, log K_Cal = -8.48, n = 0.22, CEC = 2.5. Different symbols for measured data indicate separate flow paths.

Plagioclase continues on being dissolved at higher pH values producing the dissociated $H_{3} Si O_{4}^{-}$ increasing total alkalinity (Equation 6):

\begin{array}{l} 2 C a_{0.5} N a_{0.5} A l_{1.5} S i_{2.5} O_{8} + H C O_{3}^{-} + 5.5 H_{2} O + H^{+} \to \\ 1.5 A l_{2} S i_{2} O_{5} {(O H)}_{4} + 2 H_{3} S i O_{4}^{-} + N a^{+} + C a C O_{3} & (6) \end{array}

Production of $H_{3} Si O_{4}^{-}$ (Equation 6) causes an increase in the total Si concentration (Figure 2). The increases in Na⁺ and Cl⁻ concentration from 150 to 300 km in Figure 2 are chiefly a result of the inflow from the underlying strata as the boundary concentrations are 5 mM. Since the bottom recharge Na:Cl ratio is 1, any number larger than 1 is a combination of plagioclase hydrolysis and ion exchange.

3.2. Sensitivity of the Model Output to Selected Parameters

Figure 4 depicts how sensitive the model outputs are in relation to the model inputs at different distances. Overall, the carbon concentrations depend primarily on partial pressure CO₂ and the plagioclase hydrolysis rate (Figure 4 and Appendix). These two input parameters affect a number of model outputs including Alkalinity, Ca²⁺, Na²⁺, pH, and SiO₂. Although high groundwater recharge rates indicate elevated carbon loading, the global influences of Rech_0 and Rech_50 on dissolved carbon concentration are minor.

FIGURE 4

Figure 4. Sensitivity indices and confidence intervals for different distances from the model edge. Inputs are on the horizontal axes, whereas outputs are on the vertical axis. Large SI for the DIC output indicates the most influential input parameters, chiefly pCO2 and r_Plag. The exact SI and Conf values are given in Appendix.

The solubility constants of secondary minerals including calcite do not have major implications on the DIC concentrations. This indicates that even a relatively wide range of calcite solubility values (SI_Cal variations in a range of ± 1.0) does not affect the carbon cycling to this extent as recharge CO₂ and the plagioclase dissolution rate do. The cation exchange capacity has a minor significance too, although it is often regarded as a major parameter affecting solute concentrations in large groundwater basins (Appelo, 1994; Sracek and Hirata, 2002; Gastmans et al., 2010). Indeed it influences concentrations of other species, including Na⁺ and total alkalinity.

The comparison of sensitivities at 50, 150, and 250 km show the increasing importance of weathering and cation exchange over the distance. Nevertheless, relatively similar sensitivities at these distances indicate that it does not really matter which part of the aquifer groundwater is sampled to identify the most influential global parameters.

The input parameters were evaluated in a range of values. The ranges were chosen either from the field observations (Rebouças, 1994; Gastmans et al., 2010) or literature (Appelo and Postma, 2005; Gudbrandsson et al., 2014). It is widely accepted that a saturation index (Equation 2) of a precipitating mineral phase is positive in aquifers due to kinetic inhibition of the process (Appelo and Postma, 2005). This is the reason; we used a wide range of log solubilities.

Since the field mineral dissolution rates are consistently smaller in the large scale than in the lab or small scale experiments (Gudbrandsson et al., 2014; Maher and Navarre-Sitchler, 2019), we assumed a range which enabled us to develop the concentration patterns observed in the field, chiefly the oversaturation for calcite in the intermediate range (120–160 km) of distances (Figure 3 middle).

3.3. Distribution of Carbon Flux Through the Boundaries and Transformation in the Aquifer

The DIC flux, which is a product of water recharge and DIC concentration, takes place primarily in the outcrop boundary due to relatively high recharge rates (Figure 5). The fluxes through the outcrop (Recharge from 0 to 50 km) account for approximately 75% of total carbon fluxes to the system. This is consistent with the current understanding of the GAS system (Teramoto et al., 2020). Only a minor part of the total flux is discharged.

FIGURE 5

Figure 5. Fluxes through the boundaries and the difference between input and output attributed to the transformation for 1,000 model realisations.

Thus, the transformation of DIC in the aquifer is significantly larger than the DIC discharge rate, indicating that the aquifer serves as an effective carbon sink in the wide range of parameters (Figure 5).

Extrapolating the results to the extend of the aquifer indicates that if the GAS area is 1,195,500 km², more than 4 × 10¹¹ M/yr gets immobilised within the aquifer assuming homogeneous reactions. Another assumption is that there is no carbon consumption during the leakage through the overlying basalt.

3.4. Limitations and Additional Considerations

The presented approach of estimated carbon cycling is based on the assumption that CO₂ removal is driven by plagioclase hydrolysis as it was found in highly reactive anothrosites (Banks and Frengstad, 2006) or basalts (Pogge von Strandmann et al., 2019). Although the model can be used in a variety of settings and at different spatio-temporal scales, a critical evaluation of the conceptual model is warranted.

The evaluation of the conceptual model requires water quality and mineralogical or petrographic measurements. The process of data collection is an expensive task which may not unequivocally show whether the conceptual model is correct (Enemark et al., 2019). In GAS, an alternative conceptual scenario of the carbonate dissolution process coupled with ion exchange was considered (Sracek and Hirata, 2002; Gastmans et al., 2010) as mineralogical studies showed that carbonates were being dissolved in certain locations of the aquifer. In different locations, however, there is evidence of carbonate accumulation (Gastmans et al., 2010; Hirata et al., 2011). Geochemical processes that may lead to different conceptualisations were studied and modelled at small (pore) scales (Jakobsen, 2007; Molins et al., 2012). RTMs should be useful to calculate net precipitation or dissolution.

A detailed process quantification requires not only models, but extensive multicomponent analysis. This can be accomplished in focused engineered studies, in which controlled injection and monitoring is possible. At the CarbFix2 site in Iceland a number of conservative and reactive tracers including calcium isotopes were used to study the extend of the carbonate precipitation process (Pogge von Strandmann et al., 2019). Fortunately, the geological systems can be studied by analogy across scales.

The major limitations of our model include:

• steady-state infiltration and constant carbon input. Although relatively small cyclical variations of climate in Brazil over the last few 100, 1,000 years (Cruz et al., 2005; Rodríguez-Zorro et al., 2020) were encountered, the uncertainty of our estimates remains very high;

• the plagioclase rate dissolution is pH-independent. Some studies point to the significant dependence on pH values (Gudbrandsson et al., 2014);

• the geometry of the model domain is very simple, although it captures the concentration patterns observed in NE part of the aquifer;

• only Ca-Na plagioclase is reacting with DIC, whereas in reality there is a competition among the mineral assemblage for a reactant (Oelkers et al., 2018) leading to geochemical heterogeneity and uneven evolution of porous media (Seigneur et al., 2019).

The interpreted rates of plagioclase hydrolysis rates from our model are 2–3 orders of magnitude smaller than in the small scale and laboratory studies, which is consistent with findings (Molins et al., 2012; Maher and Navarre-Sitchler, 2019). However, the weathering mineral phase should be interpreted here as a mineral assemblage which would have different reactivity than a pure mineral.

The composition of the mineral assemblage, concentrations throughout the aquifer and along the vertical profiles as well as an estimation of distribution of recharge and discharge using tracers, could be helpful in reducing uncertainties of the carbon flux. The model could further be extended to 3D representation of the aquifer.

Our numerical model with its all simplifications underpins the importance of depth specific rather than depth integrated groundwater sampling for process identification (Appelo and Postma, 2005). Even in a relatively simple chemical configuration of a homogeneous aquifer composed of plagioclase receiving infiltration CO₂, the concentration patterns reflect complexity both laterally and vertically. The effect is exacerbated by low flow rates in the aquifer.

4. Conclusions

The role of aquifers in the global carbon cycle poses significant interpretative and quantitative challenges. The interpretation is related with the conceptual model, which needs to be developed on the basis of the available datasets and the geological history.

The process quantification may be performed using a reactive transport model which encompasses geochemical and hydrogeological data. Such a model offers partial and integrated water and geochemical fluxes that can be further critically evaluated using observations. Here, we used the distribution of water quality, water residence times and qualitative information on matrix composition from the Guarani Aquifer System as an example of the continental scale aquifer. This enabled us to calculate the water and geochemical balances assuming that plagioclase weathering is thermodynamically coupled with carbonate precipitation.

A simple conceptualisation with a plagioclase mineral only suggests that 50% of carbon is immobilised as CaCO₃. By using global sensitivity analysis, we found that the DIC fluxes and weathering rates of aluminosilicate minerals have a major significance in controlling solute concentrations and, consequently, DIC outflux at the discharge zone. Other parameters including horizontal hydraulic conductivity, porosity and stability of secondary minerals play a minor or negligible role in controlling solute concentrations.

In order to better understand the role of large aquifers in carbon cycling, the effort should be placed on:

• the rate of the mineral assemblage dissolution;

• the refined water sampling along the flow paths and in the vertical profile;

• the evaluation of the DIC concentration and its temporal variations in recharge water;

• the estimation of recharge rates;

• the incorporation of isotopes (e.g., Ca²⁺ and Mg²⁺) in water quality monitoring to constrain chemical rates and hydrogeological pathways of DIC cycling;

• the distribution of cation exchange capacity in the aquifer and the sorption site composition along the flow paths.

The refinement of the models with aforementioned data would further reduce predictive uncertainty of the CO₂ immobilisation in aquifers.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

Author Contributions

KM designed the manuscript, developed the reactive model, and wrote the manuscript. LP developed the global sensitivity methodology and revised the manuscript. Both authors contributed to the article and approved the submitted version.

Supplemental Data

The numerical model repository is located on https://www.hydroshare.org/resource/b0c0815fbc4949ad9d0c21102df93a39/. Measured concentrations from GAS are from Gastmans et al. (2010).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

A number of colleagues are acknowledged for friendly and useful remarks on an earlier version of the manuscript. Two Frontiers reviewers greatly contributed to its final version.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2021.714075/full#supplementary-material

References

Aggarwal, P. K., Matsumoto, T., Sturchio, N. C., Chang, H. K., Gastmans, D., Araguas-Araguas, L. J., et al. (2015). Continental degassing of 4He by surficial discharge of deep groundwater. Nat. Geosci. 8, 35–39. doi: 10.1038/ngeo2302