Deepening roots can enhance carbonate weathering by amplifying CO₂-rich recharge

Carbonate weathering is essential in regulating atmospheric CO₂ and carbon cycle at the century timescale. Plant roots accelerate weathering by elevating soil CO₂ via respiration. It however remains poorly understood how and how much rooting characteristics (e.g., depth and density distribution) modify flow paths and weathering. We address this knowledge gap using field data from and reactive transport numerical experiments at the Konza Prairie Biological Station (Konza), Kansas (USA), a site where woody encroachment into grasslands is surmised to deepen roots.

Results indicate that deepening roots can enhance weathering in two ways. First, deepening roots can control thermodynamic limits of carbonate dissolution by regulating how much CO₂ transports vertical downward to the deeper carbonate-rich zone. The base-case data and model from Konza reveal that concentrations of Ca and dissolved inorganic carbon (DIC) are regulated by soil pCO₂ driven by the seasonal soil respiration. This relationship can be encapsulated in equations derived in this work describing the dependence of Ca and DIC on temperature and soil CO₂. The relationship can explain spring water Ca and DIC concentrations from multiple carbonate-dominated catchments. Second, numerical experiments show that roots control weathering rates by regulating recharge (or vertical water fluxes) into the deeper carbonate zone and export reaction products at dissolution equilibrium. The numerical experiments explored the potential effects of partitioning 40 % of infiltrated water to depth in woodlands compared to 5 % in grasslands. Soil CO₂ data suggest relatively similar soil CO₂ distribution over depth, which in woodlands and grasslands leads only to 1 % to ∼ 12 % difference in weathering rates if flow partitioning was kept the same between the two land covers. In contrast, deepening roots can enhance weathering by ∼ 17 % to 200 % as infiltration rates increased from 3.7 × 10⁻² to 3.7 m/a. Weathering rates in these cases however are more than an order of magnitude higher than a case without roots at all, underscoring the essential role of roots in general. Numerical experiments also indicate that weathering fronts in woodlands propagated > 2 times deeper compared to grasslands after 300 years at an infiltration rate of 0.37 m/a. These differences in weathering fronts are ultimately caused by the differences in the contact times of CO₂-charged water with carbonate in the deep subsurface. Within the limitation of modeling exercises, these data and numerical experiments prompt the hypothesis that (1) deepening roots in woodlands can enhance carbonate weathering by promoting recharge and CO₂–carbonate contact in the deep subsurface and (2) the hydrological impacts of rooting characteristics can be more influential than those of soil CO₂ distribution in modulating weathering rates. We call for colocated characterizations of roots, subsurface structure, and soil CO₂ levels, as well as their linkage to water and water chemistry. These measurements will be essential to illuminate feedback mechanisms of land cover changes, chemical weathering, global carbon cycle, and climate.

1 Introduction

Carbonate weathering has long been considered negligible as a long-term control of atmospheric CO₂ (0.5 to 1 Ma; Berner and Berner, 2012; Winnick and Maher, 2018). Recent studies, however, have underscored its significance in controlling the global carbon cycle at the century timescale that is relevant to modern climate change, owing to its rapid dissolution, its fast response to perturbations, and the order-of-magnitude-higher carbon store in carbonate reservoirs compared to the atmosphere (Gaillardet et al., 2019; Sullivan et al., 2019b). Carbonate weathering is influenced by many factors, including temperature (Romero-Mujalli et al., 2019b), hydrological regimes (Romero-Mujalli et al., 2019a; Wen and Li, 2018), and soil CO₂ concentrations (Covington et al., 2015) arising from different vegetation types (Calmels et al., 2014). Rapid alteration to any of these factors, either human or climate induced, may change global carbonate weathering fluxes and lead to a departure from the current global atmospheric CO₂ level. This is particularly important given that about 7 %–12 % of the Earth's continental area is carbonate based and about 25 % of the global population completely or partially depend on waters from karst aquifers (Hartmann et al., 2014).

Plant roots have long been recognized as a dominant biotic driver of chemical weathering and the global carbon cycle (Berner, 1992; Beerling et al., 1998; Brantley et al., 2017a). The growth of forests has been documented to elevate soil pCO₂ and amplify dissolved inorganic carbon (DIC) fluxes (Berner, 1997; Andrews and Schlesinger, 2001). Rooting structure can influence weathering in two ways (Fig. 1). First, rooting systems (e.g., grasslands, shrublands, and woodlands) may affect the distribution of soil carbon (both organic and inorganic), microbe biomass, and soil respiration (Drever, 1994; Jackson et al., 1996; Billings et al., 2018). The relatively deep root distributions of shrublands compared to grasslands may lead to deeper soil carbon profiles (Jackson et al., 1996; Jobbagy and Jackson, 2000), which may help elevate the deep CO₂ and acidity that determine carbonate solubility and weathering rates.

Second, plant roots may affect soil structure and hydrological processes. Root trenching and etching can develop porosity (Mottershead et al., 2003; Hasenmueller et al., 2017). Root death and decay can promote the generation of macropores or, more specifically, biopores with connected networks (Angers and Caron, 1998; Zhang et al., 2015). Root channels have been estimated to account for about 70 % of the total described macropores (Noguchi et al., 1997; Beven and Germann, 2013) and for over 70 % of water fluxes through soils (Watson and Luxmoore, 1986). In grasslands, the lateral, dense spread of roots in upper soil layers promotes the formation of horizontally oriented macropores that support near-surface lateral flow (Cheng et al., 2011). Highly dense fine roots also increase the abundance of organic matter and promote granular or sandy texture soil aggregates that facilitate shallow, near-surface water flow (Oades, 1993; Nippert et al., 2012). In contrast, in shrublands and forests, generally deeper and thicker roots tend to promote a high abundance of macropores and high connectivity to the deep subsurface (Canadell et al., 1996; Nardini et al., 2016), enhancing the drainage of water to the depth (Pawlik et al., 2016).

It is generally known that rooting characteristics vary among plant species and are critical in regulating water budgets, flow paths, and storage (Sadras, 2003; Nepstad et al., 1994; Jackson et al., 1996; Cheng et al., 2011; Brunner et al., 2015; Fan et al., 2017). Existing studies however have primarily focused on the role of soil CO₂ and organic acids (Drever, 1994; Lawrence et al., 2014; Gaillardet et al., 2019; Hauser et al., 2020). Systematic studies on coupled effects of hydrological flow paths and soil CO₂ distribution are missing, owing to the limitation in data that detail rooting effects on flow partitioning and complex hydrological–biogeochemical interactions (Li et al., 2020). Here we ask the following questions. How and to what degree do rooting characteristics influence carbonate weathering when considering both flow partitioning and soil CO₂ distribution? Which factor (flow partitioning or soil CO₂ distribution) predominantly controls weathering? We hypothesized that deepening roots in woodlands enhance carbonate weathering by promoting CO₂-enriched recharge into the deep, carbonate-abundant subsurface (Fig. 1).

We tested the hypotheses by a series of numerical experiments of reactive transport processes based on water chemistry data from an upland watershed in the Konza Prairie Biological Station, a tallgrass prairie and one of the Long-Term Ecological Research (LTER) sites in the US (Fig. S1, Macpherson and Sullivan, 2019b; Vero et al., 2018). We used the calibrated model to carry out numerical experiments for two end-members of vegetation covers, grasslands and woodlands, under flow-partitioning conditions that are characteristic of their rooting structure. These experiments differentiated the impacts of biogeochemical and hydrological drivers and bracketed the range of their potential impacts on weathering, thus providing insights on the missing quantitative link between rooting structure and chemical weathering. We recognize that rooting characteristics can have multiple influences on water flow paths and the water budget, for example, via water uptake and transpiration (Sadras, 2003; Fan et al., 2017; Pierret et al., 2016). This study focuses primarily on their potential influence via the alteration of hydrological flow paths.

Figure 1 A conceptual diagram of hydro-biogeochemical interactions in the grassland (a) and woodland (b). The shallow and dense fine roots in the grassland promote lateral macropore development and lateral water flow. In contrast, the woodlands induce vertical macropore development that supports vertical flow (recharge) into the deep, calcite-abundant subsurface compared to the grassland. The gray color gradient reflects the calcite abundance with more calcite in depth. M_L and M_G with a unit of moles per year (= flow rate Q (m/a) × species concentration C (mol/m³) × cross-section area (m²)) represent mass fluxes from lateral flow (soil water) and from vertical flow into groundwater, respectively. Soil water and groundwater were assumed to eventually flow into stream, adding up to the total discharge (infiltration) rate $Q_{\text{T}}=Q_{\text{L}}+Q_{\text{G}}$.

2 Research site and data sources

2.1 Site description

Details of the Konza site are in the Supplement and references therein. Here we provide brief information relevant to this work. Konza is a mesic native grassland where experimentally manipulated, long-term burning regimes have led to woody encroachment in up to 70 % of the catchment area in some catchments. The mean annual temperature and precipitation are 13 ^∘C and ∼ 835 mm, respectively (Tsypin and Macpherson, 2012). The bedrock contains repeating Permian couplets of limestone (1–2 m thick) and mudstone (2–4 m) (Macpherson et al., 2008). The limestone is primarily calcite with traces of dolomite, while the mudstone is dominated by illite, chlorite, and mixed layers of chlorite–illite and chlorite–vermiculite, varying in abundance from major to trace amounts. With an average thickness of 1–2 m in the lowlands, soils mostly have carbonate less than 25 % (Macpherson et al., 2008). Data suggest that the Konza landscape is undergoing a hydrogeochemical transition that coincides with and may be driven by woody encroachment. Parallel to these changes is a detectable decline in streamflow and an increase in weathering rates (Macpherson and Sullivan, 2019b) and groundwater pCO₂. The concentration–discharge relationships have exhibited chemodynamic patterns (i.e., solute concentrations are sensitive to changes in discharge) for geogenic species (e.g., Mg and Na) in woody-encroached sites compared to grass sites. Sullivan et al. (2019a) hypothesized that concentration–discharge relationships may be affected by woody species with deeper roots, which altered flow paths and mineral–water interactions (Fig. 1). We focus on the upland watershed N04d in Konza (Fig. S1) that has experienced a 4-year burning interval since 1990 and has seen considerable woody encroachment and changes in hydrologic fluxes (Sullivan et al., 2019a).

Table 1 Key reactions and kinetic and thermodynamic parameters.

^a Values of K_eq were interpolated using the EQ3/6 database (Wolery et al., 1990), except Reactions (1) and (4) (i.e., K₁ and K₄). ^b CO₂(aq) concentrations were calculated through $C_{{\mathrm{CO}}_{\mathrm{2}}}$(aq) = K₁pCO₂. The prescribed $C_{{\mathrm{CO}}_{\mathrm{2}}}$(aq) values were used as equilibrium constants in CrunchTope to describe how much soil CO₂ was available for weathering. ^c Calcite in the upper soil (above horizon B) is mixed with other minerals and has small particle size (Macpherson and Sullivan, 2019b), is relatively impure, and therefore has a lower K_eq value than those at depth with relatively pure calcite. The K_eq of the impure calcite was calibrated by fitting field data of Ca and alkalinity. ^d The kinetic rate parameters and specific surface areas were from Palandri and Kharaka (2004), except Reaction (0). The kinetic rate constant of the source CO₂(g^∗) dissolution (i.e., soil respiration rate constant, Reaction 0) was from Bengtson and Bengtsson (2007), Ahrens et al. (2015), and Carey et al. (2016); the specific surface area was referred to that of soil organic carbon (Pennell et al., 1995). ^e These reactions were only used in the base-case model as they occurred in upper soils in Konza. In the later numerical experiments, these reactions were not included so as to focus on carbonate weathering.

2.2 Data sources

Daily total meteoric precipitation and evapotranspiration were from the Konza data website (http://lter.konza.ksu.edu/data, last access: 22 November 2020). Wet chemistry deposition data were from the National Atmospheric Deposition Program (NADP, http://nadp.slh.wisc.edu/, last access: 22 November 2020). Data used in this work include monthly data of soil gases (at depths of 16, 84, and 152 cm from land surface), soil water (17 and 152 cm from land surface), and groundwater (366 cm from land surface) in 2009 and 2010 (Tsypin and Macpherson, 2012) (Fig. 2a). The sampling points were about 30 m away from the stream. More information on field and laboratory methods was included in the Supplement.

3 Reactive transport modeling

3.1 Base case: 1-D reactive transport model for the Konza grassland

A 1-D reactive transport model was developed using the code CrunchTope (Steefel et al., 2015). The code solves mass conservation equations integrating advective and diffusive/dispersive transport and geochemical reactions. It has been extensively used in understanding mineral dissolution, chemical weathering, and biogeochemical reactions (e.g., Lawrence et al., 2014; Wen et al., 2016; Deng et al., 2017). In this study, the base case had a porosity of 0.48 and a depth of 366.0 cm at a resolution 1.0 cm. Soil temperature was assumed to decrease linearly from 17 ^∘C at the land surface to 8 ^∘C at 366.0 cm, within typical ranges of field measurements (Tsypin and Macpherson, 2012). Detailed setup of domain, soil mineralogy, initial condition, and precipitation chemistry are in the Supplement.

3.1.1 Representation of soil CO₂

The model does not explicitly simulate soil respiration (microbial activities and root respiration) that produces soil CO₂. Instead, it approximates these processes by having a CO₂ source phase CO₂(g^∗) that continuously releases CO₂(g) and reproduced the observed soil pCO₂ levels at a kinetic rate constant of 10⁻⁹ mol/m²/s (Reactions 0–1 in Table 1). This value is at the low end of the reported soil respiration rates (10⁻⁹–10⁻⁵ mol/m²/s) (Bengtson and Bengtsson, 2007; Ahrens et al., 2015; Carey et al., 2016). The dissolution of CO₂(g) into CO₂(aq) follows Henry's law $C_{{\text{CO}}_{\mathrm{2}}}\left(\text{aq}\right)=K_{\mathrm{1}}p{\text{CO}}_{\mathrm{2}}$. Here K₁ is the equilibrium constant of Reaction (1), which equals to Henry's law constant. The extent of CO₂(g) dissolution was constrained by $C_{{\text{CO}}_{\mathrm{2}}\left(\text{aq}\right)}$, which was estimated using temperature-dependent K₁ (following the van 't Hoff equation in Eq. S1) and measured soil CO₂ data at different horizons (Table 2). These $C_{{\text{CO}}_{\mathrm{2}}\left(\text{aq}\right)}$ values were then linearly interpolated for individual grid blocks in the model. Finally, these prescribed $C_{{\text{CO}}_{\mathrm{2}}\left(\text{aq}\right)}$ values were used as the equilibrium constants of the coupled Reactions (0–1) in the form of CO₂(g${}^{\ast })\leftrightarrow {\text{CO}}_{\mathrm{2}}(\text{aq})$, such that the soil CO₂ values at different depth were represented (Sect. 4.1). In the base case, the soil profile of $C_{{\text{CO}}_{\mathrm{2}}}\left(\text{aq}\right)$ was updated monthly based on the monthly soil CO₂ data (Table 2). More details about the implementation in CrunchTope are included in the Supplement.

Table 2 Measured CO₂(g) at different depths and corresponding estimated $C_{{\text{CO}}_{\mathrm{2}}}$(aq).

^a Monthly measured soil CO₂ data for the Konza grassland (base case) were from Tsypin and Macpherson (2012); the soil CO₂ in the grassland experiments was averaged from monthly measurements. More information on measurements at the Konza grassland is detailed in the Supplement. The annual-average soil CO₂ in the woodland experiments was from the forested Calhoun site in South Carolina (Billings et al., 2018). ^b CO₂(aq) values were estimated using Henry's law: $C_{{\mathrm{CO}}_{\mathrm{2}}}$(aq) = K₁pCO₂; the temperature-dependent K₁ was calculated following Eq. (S1). The $C_{{\mathrm{CO}}_{\mathrm{2}}}$(aq) values at different soil depths were used to prescribe the available soil CO₂ for chemical weathering. ^c Soil temperature was estimated from the soil water and shallow groundwater temperature (Tsypin and Macpherson, 2012; Billings et al., 2018).

3.1.2 Reactions

In the model, the upper soil layers have more anorthite (CaAl₂Si₂O₈(s)) and K-feldspar (KAlSi₃O₈(s)), and the deeper subsurface contains more calcite (Table S1). The calcite volume increases from 0 % in the upper soil layer to 10 % in the deep subsurface. Table 1 summarizes reactions and thermodynamic and kinetic parameters. Soil CO₂ increases pore water acidity (Reactions 0–2) and accelerates mineral dissolution (Reactions 3–6). Silicate dissolution leads to the precipitation of clay (represented by kaolinite in Reaction 7). These reactions were included in the base case to reproduce field data. The kinetics follows the transition state theory (TST) rate law (Plummer et al., 1978) $r=\text{kA}\left(\mathrm{1}-\frac{\text{IAP}}{K_{\text{eq}}}\right)$, where k is the kinetic rate constant (mol/m²/s), A is the mineral surface area per unit volume (m²/m³), IAP is the ion activity product, and K_eq is the equilibrium constant. The term IAP ∕ K_eq quantifies the extent of disequilibrium: values close to 0 suggest far from equilibrium, whereas values close to 1.0 indicate close to equilibrium. At Konza, calcite in the upper soil (above horizon B) is mixed with other minerals, has small particle size, and is considered impure (Macpherson and Sullivan, 2019a) and therefore has a lower K_eq value than those at depth with relatively pure calcite. The K_eq of impure calcite (K₄) was calibrated by fitting field data of Ca and alkalinity. To reproduce the observed soil CO₂ profile, the soil CO₂ production rate (mol/m²/a) at the domain scale (calculated by the mass change in the solid phase CO₂(g^∗) over time) was assumed to increase with infiltration rates (Fig. S2). This is consistent with field observations that soil CO₂ production rate and efflux may increase with rainfall in grassland and forest ecosystems (Harper et al., 2005; Patrick et al., 2007; Wu et al., 2011; Vargas et al., 2012; Jiang et al., 2013). For example, Zhou et al. (2009) documented soil CO₂ production rates increasing from 3.2 to 63.0 mol/m²/a when the annual precipitation increased from 400 to 1200 mm. Wu et al. (2011) showed that increasing precipitation from 5 to 2148 mm enhanced soil respiration by 40 % and that a global increase of 2 mm precipitation per decade may lead to an increase of 3.8 mol/m²/a for soil CO₂ production. The simulated soil CO₂ production rates across different infiltration rates here (∼ 0.1–10 mol/m²/a shown in Fig. S2) were close to the reported belowground net primary production (belowground NPP) of typical ecosystems: ∼ 0.8–100 mol/m²/a for grasslands (Gill et al., 2002) and ∼ 0.4 –40 mol/m²/a for woodlands (Aragão et al., 2009).

Table 3 Physical and geochemical characteristics in numerical experiments.

3.1.3 Flow partitioning

Rainwater enters soil columns at the annual infiltration rate of 0.37 m/a, estimated based on the difference between measured precipitation (0.88 m/a) and evapotranspiration (0.51 m/a). At 50 cm, a lateral flow Q_L (soil water) exited the soil column to the stream at 0.35 m/a. The rest recharged the deeper domain beyond 50 cm (to the groundwater system) at 0.02 m/a (∼ 2 % of precipitation) and became groundwater (Fig. 1a), a conservative value compared to 2 %–15 % reported in another study (Steward et al., 2011). The groundwater flow Q_G eventually came out at 366.0 cm and was assumed to enter the stream as part of discharge ($=Q_{\text{L}}+Q_{\text{G}}$). The flow field was implemented in CrunchTope using the “PUMP” option.

3.1.4 Calibration

We used monthly alkalinity and Ca concentration data (Sect. 2.2) for model calibration. The monthly Nash–Sutcliffe efficiency (NSE) that quantified the residual variance of modeling output compared to measurements was used for model performance evaluation (Moriasi et al., 2007). NSE values higher than 0.5 are considered acceptable.

3.2 Numerical experiments

Numerical experiments were set up for grasslands and woodlands. The base case from Konza was used to represent grasslands: data from the Calhoun site (the Calhoun Critical Zone Observatory in South Carolina, USA) were used as representative for woody sites (Billings et al., 2018). Grasslands are typically characterized by a high proportion of horizontal macropores induced by dense, lateral-spread of roots mostly at depths less than 0.8 m (Jackson et al., 1996; Frank et al., 2010). These characteristics promote lateral flow (Q_L) at the shallow subsurface (Fig. 1a). At Konza, over 90 % of grass roots were at the top 0.5 m, leading to high hydrologic conductivity in top soils (Nippert et al., 2012). In the woodland, a greater proportion of deep roots enhances vertical macropore development (Canadell et al., 1996; Nardini et al., 2016), reduces permeability contrasts at different depths (Vergani and Graf, 2016), and is thought to facilitate more vertical water flow to the depth (Q_G). At the Calhoun site, over 50 % of roots are in the top 0.5 m in woodlands, with the rest penetrating deeper (Jackson et al., 1996; Eberbach, 2003; Billings et al., 2018).

The experiments aimed to compare the general, averaged behaviors rather than event-scale dynamics so the annual-average soil CO₂ data and corresponding prescribed CO₂(aq) concentrations were used (Table 2). The experiments focused on calcite weathering (Reactions 0–4) and excluded silicate weathering reactions (Reactions 5–7). The mineral-dissolution parameters from the base case were used for all experiments. We compared the relative significance of the hydrological (i.e., lateral versus vertical flow partitioning) vs. respiratory (i.e., CO₂ generation) influences of deepening roots. Though other potential differences might be induced by deepening roots (e.g., water uptake, water table, and transpiration) (Li, 2019) and influence weathering, we assumed they remain constant across the grassland and woodland simulations to examine the relative influences of flow path vs. CO₂ generation on carbonate weathering.

3.2.1 Hydrological and biogeochemical differences in grasslands and woodlands

Flow partitioning between lateral shallow flow and vertical recharge flow is challenging to quantify and is subject to large uncertainties under diverse climate, lithology, and land cover conditions. The ratios of lateral flow in upper soils versus the total flow inferred from a tracer study in a grassland vary from ∼ 70 % to ∼ 95 % (Weiler and Naef, 2003). Harman and Cosans (2019) found that the lateral flow rate at upper soils over the overall infiltration can vary between 50 % and 95 %. Deeper roots in woodlands can increase deep soil permeability by over 1 order of magnitude (Vergani and Graf, 2016). Assuming that the vertical, recharged flow water ultimately leaves the watershed as baseflow, the ratio of the lateral versus vertical flow has been reported with a wide range. In forests such as Shale Hills in Pennsylvania and Coal Creek in Colorado, ∼ 7 %–20 % of stream discharge is from groundwater, presumably recharged by vertical flow (Li et al., 2017; Zhi et al., 2019). Values of Q_G ∕ Q_T estimated through base flow separation vary from 20 % to 90 % in forest/wood-dominated watersheds (Price, 2011), often negatively correlating with the proportion of grasslands (Mazvimavi et al., 2004).

Soil respiration rates can vary between 10⁻⁹ and 10⁻⁵ mol/m²/s in both grasslands and woodlands (Bengtson and Bengtsson, 2007; Ahrens et al., 2015; Carey et al., 2016). Soil CO₂ levels may vary by 2–3 orders of magnitude depending on vegetation type and climate conditions (Neff and Hooper, 2002; Breecker et al., 2010). Soil CO₂ levels are further complicated by their dependence on topographic position, soil depth, and soil moisture, all of which determine the magnitude of microbial and root activities and CO₂ diffusion (Hasenmueller et al., 2015; Billings et al., 2018). There is however no consistent evidence suggesting which land cover exhibits higher soil respiration rate or soil CO₂ level. Below we describe details of the numerical experiments (Table 3) exploring the influence of hydrological versus biogeochemical impacts of roots.

Figure 2 (a) Schematic representation of sampling depths in the Konza grassland; corresponding monthly dynamics of (b) CO₂(aq) and CO₂(g), (c) alkalinity, (d) DIC, and (e) Ca concentrations. Lines represent modeling outputs at the corresponding sampling depth of monthly field measurements (dots), including soil water at horizons A (16 cm) and B (152 cm), as well as groundwater (366 cm). The lines of CO₂(aq) and CO₂(g) in panel (b) overlapped. Note that there are no DIC data so no dots in panel (d). The temporal trends of alkalinity, DIC, and Ca mirrored those of soil CO₂, indicating its predominant control on weathering.

3.3 Three numerical scenarios

Each scenario in Table 3 includes a grassland and woodland case, with their respective profiles of calcite distribution and soil pCO₂ kept the same (columns 3 and 4 in Table 3). The only difference in different scenarios is the flow partitioning (column 5). Soil pCO₂ in the grassland (red line in column 4) was set to reflect the annual average of the Konza site. Soil pCO₂ (blue line) in the woodland was set to reflect the annual average from the forest-dominant Calhoun site (Billings et al., 2018). In all scenarios, we assumed that grasslands and woodlands had the same total CO₂(g^∗) producing CO₂ gas and CO₂(aq) but differed in depth distributions. The CO₂(g^∗) depth distribution was constrained by CO₂(g) field data (Table 2). The distribution of CO₂(g^∗) in grasslands is slightly steeper than that in the woodland, with higher density of roots and more abundant CO₂(g^∗) in the upper soils. Below the rooting depth, CO₂(g^∗) was assumed to be smaller (by 10 times) to represent the potential soil CO₂ sources from microbial activities (Billings et al., 2018).

Figure 3 Simulated depth profiles of (A) calcite (vol %) and (B) soil CO₂ production rate; (C) CO₂(g) (%), (D) DIC, and (E) Ca at 0.37 m/a; and (F) effluent Ca concentrations at different infiltration rates. From top to bottom rows are Scenario 1 (Grass_PF and Wood_PF, with flow partitioning, first row), Scenario 2 (Grass_VF and Wood_VF, with 100 % vertical flow, second row), and Scenario 3 (Grass_HF and Wood_HF, with 100 % lateral flow, last row). Red and blue colors present grassland and woodland, respectively. Lines and empty circles represent modeling outputs, while filled circles with error bar in C are the annual-average CO₂(g) data (in Table 2). Arrows in A indicate flow conditions. In F, soil water and groundwater refer to concentrations at the lateral (50 cm) and vertical (366 cm) outlets, respectively. Higher fractions of vertical flow in Wood_PF led to higher stream Ca concentration compared to Grass_PF. In Scenario 2 and 3 without flow partitioning, stream Ca concentrations were similar. The concentrations of DIC versus discharge are very similar to the Ca concentrations in F.

Scenario 1 considered flow partitioning (Table 3). With the large permeability contrast of soil and bedrock (over 4 orders of magnitude) in Konza (Macpherson, 1996), the Grass_PF case (with Q_L=95 % × Q_T and Q_G=5 %  × Q_T) represents an end-member case for the grassland. The woodland (Wood_PF) case was set to have 60 % lateral flow and 40 % vertical flow into the deeper subsurface. This groundwater percentage is at the high end of flow partitioning and serves as an end-member case for woodlands (Vergani and Graf, 2016). Scenarios 2 and 3 had no flow partitioning. Scenario 2 had two cases with 100 % vertical flow via the bottom outlet (VF; Wood_VF and Grass_VF); Scenario 3 had two cases with 100 % horizontal flow (HF; Wood_HF and Grass_HF) via the shallow outlet at 50 cm. These cases represent the end-member flow cases with 100 % lateral flow or 100 % vertical flow. In addition, because the two cases have the same flow scheme, they enable the differentiation of effects of soil CO₂ distribution versus hydrology differences. All scenarios were run under infiltration rates from 3.7 × 10⁻² to 3.7 m/a (10⁻⁴–10⁻² m/d), the observed daily variation range at Konza. This was to explore the role of flow regimes and identify conditions where the most and least significant differences occur.

Each case was run until steady state, when concentrations at the domain outlet became constant (within ± 5 %) over time. The time to reach steady state varied from 0.1 to 30 a, depending on infiltration rates. The lateral flux (soil water, M_L=Q_L × C_L) and vertical fluxes (groundwater, M_G=Q_G × C_G) were calculated at 50 and 366 cm, in addition to total fluxes (weathering rates, M_T). These fluxes multiplied with unit cross-section area (m²) convert into rates in units of moles per year (mol/a).

3.4 Carbonate weathering over century timescales: soil property evolution

To compare the propagation of weathering fronts over longer timescales, we carried out two 300-year simulations for Scenario 1 with flow partitioning under the base-case infiltration rate of 0.37 m/a (i.e., Grass_PF and Wood_PF in Table 3). During this long-term simulation, we updated calcite volume, porosity, and permeability. The calcite volume changes were updated in each time step based on corresponding mass changes, which were used to update porosity. Permeability changes were updated based on changes in local porosity following the Kozeny–Carman equation: $\frac{{\mathit{\kappa }}_i}{{\mathit{\kappa }}_{i,\mathrm{0}}}={\left(\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\varnothing }}_{i,\mathrm{0}}}\right)}^{\mathrm{3}}\times {\left(\frac{\mathrm{1}-{\mathrm{\varnothing }}_{i,\mathrm{0}}}{\mathrm{1}-{\mathrm{\varnothing }}_i}\right)}^{\mathrm{2}}$ (Kozeny, 1927; Costa, 2006), where κ_i and ∅_i are permeability and porosity in grid i at time t, and κ_i,0 and ∅_i,0 are the initial permeability and porosity, respectively.

3.5 Quantification of weathering rates and their dependence on CO₂–carbonate contact

To quantify the overall weathering rates and CO₂–carbonate contact in each scenario, we used the framework from a previously developed upscaled rate law for dissolution of spatially heterogeneously distributed minerals (Wen and Li, 2018). The rate law says that three characteristic times are important. The equilibrium time τ_eq represents the characteristic timescale of mineral dissolution to reach equilibrium in a well-mixed system. The residence time τ_a, i.e., the timescale of advection, quantifies the overall water contact time with the whole domain. It was calculated by the product of domain length (L) and porosity divided by the overall infiltration rate (Q_T): ${\mathit{\tau }}_{\text{a}}=\frac{L\mathit{\varphi }}{Q_{\text{T}}}$. The reactive transport time τ_ad,r quantifies the water contact time with calcite as influenced by both advection and diffusion/dispersion. The upscaled rate law is as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& R_{\text{calcite}}=k_{\text{calcite}}{A}_{\text{T}}\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\tau }}_{\text{eq}}}{{\mathit{\tau }}_{\text{a}}}}\right)\right]\\ & \times {\left\{\mathrm{1}-\exp \left[-L\left({\displaystyle \frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}}\right)\right]\right\}}^{\mathit{\alpha }}\end{array}.\end{array}$$

Here k is the intrinsic rate constant measured for a mineral in a well-mixed reactor, A_T is the total surface area, L is domain length, α is geostatistical characteristics of spatial heterogeneity, and $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ is the reactive time ratio quantifying the relative magnitude of the water contact time with the whole domain versus the contact time with the reacting mineral. This rate law consists of two parts: the effective dissolution rates in homogeneous media represented by $k{A}_{\text{T}}\left[\mathrm{1}-\exp \left(-\frac{{\mathit{\tau }}_{\text{eq}}}{{\mathit{\tau }}_{\text{a}}}\right)\right]$ and the heterogeneity factor ${\left\{\mathrm{1}-\exp \left[-L\left(\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}\right)\right]\right\}}^{\mathit{\alpha }}$ that quantifies effects of preferential flow paths arising from heterogeneous distribution of minerals. When $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}>\mathrm{1}$, the water contact time with calcite zone is small, meaning the water is replenished quickly compared to the whole domain, leading to higher CO₂–carbonate interactions. In contrast, a small $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ ratio (<1) reflects that water is replenished slowly in the reactive calcite zones, leading to less CO₂–carbonate contact. These different timescales were calculated for Scenario 1–3 based on the flow characteristics and dissolution thermodynamics and kinetics, as detailed in the Supplement. Values of τ_eq, τ_a, and τ_ad,r for all experiments are listed in Table S2. Note that numerical experiments in Scenario 1–3 focused on the short-term scale, with negligible changes in the solid phase.

4 Results

4.1 The thermodynamics of carbonate dissolution: grassland at Konza as the base-case scenario

The calibrated model reproduced the temporal dynamics with a Nash–Sutcliffe efficiency (NSE) value > 0.6 and was considered satisfactory (Fig. 2). Note that the y axis is inverted to display upper soils at the top and deep soils at the bottom to be consistent with their subsurface position shown in Fig. 2a. The measured CO₂(g) varied between 0.24 % and 7.30 % (Fig. 2b right axis), 1 to 2 orders of magnitude higher than the atmospheric level of 0.04 %. The estimated CO₂(aq) (Fig. 2b left axis) generally increased with depth except in July and August when horizon B was at peak concentration. The timing of the peaks and valleys varied in different horizons. The CO₂(aq) reached maxima in summer in soil horizons A and B and decreased to less than 0.5 mM in winter and spring. The groundwater CO₂(aq) exhibited a delayed peak in September and October and dampened seasonal variation compared to the soil horizons. The temporal trends of alkalinity, DIC, and Ca in groundwater mirrored those of soil CO₂ at their corresponding depths, indicating the predominate control of soil CO₂ on carbonate weathering. The groundwater concentrations of these species were also higher than soil concentrations. The simulated groundwater DIC (approximately summation of CO₂(aq) and alkalinity) was > 6 times higher than that in upper soil (∼ 1.0 mM). The dissolved mineral volume was negligible for the simulation period (<0.5 % v∕v). Sensitivity analysis revealed that changes in flow velocities influenced concentrations in horizon A where anorthite is the dominating dissolving mineral (0–1.8 m); their effects are negligible in horizon B and groundwater where fast-dissolving calcite rapidly approaches equilibrium.

Several measurements/parameters were important in reproducing data. These include soil CO₂, which determined the CO₂(aq) level and its spatial variation, and equilibrium constant (K_eq) of calcite dissolution. Imposition of monthly variations and depth distributions of soil CO₂ were essential to capture the variation of alkalinity and Ca data at different horizons. The imposition of calcite K_eq was also critical for reproducing Ca concentrations. Impurities were suggested to affect K_eq of natural calcite by a factor of ∼ 2.0 (Macpherson and Sullivan, 2019a). Calcite K_eq in upper soils had to be reduced by a factor of 3.8 in the model to reproduce concentrations in horizon A. The alkalinity and Ca concentrations were not sensitive to kinetic parameters nor precipitation, because carbonate dissolution rapidly approaches equilibrium.

4.2 Numerical experiments: the significance of hydrological flow partitioning

Scenario 1 for hydro-biogeochemical effects with flow partitioning (Wood_PF and Grass_PF). Figures 3A1–E1 show depth profiles of calcite and soil CO₂ production rates and steady-state concentrations of reaction products. The soil CO₂ production rate was highest in the upper soil at around 10^−6.5 mol/s and decreased to ∼ 10^−9.5 mol/s at 366 cm (Fig. 3B1), consistent with the decline with soil depth observed in natural systems. The CO₂(g) level (released from CO₂(g^∗)) increased with depth due to lower efflux into the atmosphere in deeper zone and downward fluxes of CO₂–charge water from the upper soil (Fig. 3C1). Concentrations of reaction products (Ca, DIC) were lower in the top 40 cm, reflecting the lower carbonate-mineral background level, and higher at depths over 60 cm. The transition occurred between 35 and 60 cm in the vicinity of the calcite–no-calcite interface at 55 cm, where concentrations of Ca and DIC increased abruptly until reaching equilibrium. This thin transition was driven by fast calcite dissolution and rapid approach to equilibrium, resulting in a short equilibrium distance. The equilibrated DIC and Ca concentrations below 60 cm followed the similar increasing trend of CO₂(g) with depth in the deep zone (Fig. 3C1–E1).

Figure 3F1 shows that Ca concentrations in soil water (light color) were lower than groundwater (dark color) and varied with infiltration rates. The difference between soil water and groundwater Ca concentrations was relatively small at low infiltration rates because both reached equilibrium but diverged at high infiltration rates. Higher infiltration rates diluted soil water but not as much for groundwater. As expected, concentrations in the stream water, a mixture of soil water and groundwater (solid line), were in between these values but closely resembled soil water in the grassland. In both cases, stream concentration decreased as infiltration increased, indicating a dilution concentration–discharge relationship.

Scenario 2–3 for biogeochemical effects (without flow partitioning). Scenarios 2 and 3 were end-member cases that bracketed the range of rooting effects. Calcite and soil CO₂ were distributed the same way as their corresponding flow partitioning (PF) cases (Table 3). The Wood_VF and Grass_VF cases had 100 % flow going downward via the deeper calcite zone maximizing the CO₂–water–calcite contact. In the Wood_HF and Grass_HF cases (Table 3), all water exited at 50 cm, bypassing the deeper calcite-abundant zone and minimizing the CO₂–water–calcite contact. Figure 3A2–F2 and 3A3–F3 show the VF (Wood_VF and Grass_VF) and HF (Wood_HF and Grass_HF) cases, respectively. Similar to the flow partitioning cases, the concentrations of reaction products were low in the shallow zone and increased over 10 times within a short distance ∼ 5 cm at the depth of ∼ 50 cm. The woodland cases increased slightly more than the grassland cases because of the slightly steeper soil CO₂ distribution (Fig. 3C2 and 3C3). Figure 3F2 indicates that the effluent Ca concentrations were slightly higher in Wood_VF due to the high soil CO₂ level at the bottom outlet. The Grass_HF and Wood_HF case almost had the same effluent Ca concentrations at the upper soil (Fig. 3F3).

Concentration–discharge relationship and weathering rates in all cases. The VF cases had the highest effluent concentrations and weathering rates, whereas the HF cases had lowest concentrations and weathering rates, and the Grass_PF and Wood_PF cases fell in between (Fig. 4a, b). This is because the VF cases maximized the CO₂–calcite contact with 100 % flow through the calcite zone, whereas the HF cases had minimum CO₂–calcite contact with most water bypassing the calcite zone (column 3–4 in Table 3). The Grass_PF and Wood_PF cases allowed different extent of contact prescribed by the amount of flow via the calcite zone. A case run without soil respiration (i.e., no CO₂(g^∗), not shown) indicated that Ca and DIC concentrations were more than an order of magnitude lower than cases with soil respiration. In addition, the PF and HF cases generally showed dilution patterns with concentration decreasing with infiltration rates, as compared to the VF cases where a chemostatic pattern emerged with almost no changes as infiltration rates vary. This is because in the VF cases the concentrations mostly reached equilibrium concentration. In the PF and HF cases, a large proportion of water flows through soils with negligible calcite where the water moves away from equilibrium as infiltration rates increase.

Figure 4 (a–b) Stream water Ca and DIC concentrations, and (c–d) stream fluxes from all scenarios. Stream water was the overall effluent from soil water and groundwater. The differences caused by hydrological differences (VF, HF, and PF) were much larger than the differences within each pair with the same flow partitioning, indicating significant hydrological impacts on weathering.

Although weathering rates generally increased with infiltration rates, the woodland increased more (4.6 × 10⁻² to 2.4 mol/a in Wood_PF). The Ca fluxes in the VF cases (100 % vertical flow) were higher than flow partitioning cases because they enabled maximum CO₂-charged water content with unweathered calcite at depth. Within each pair with the same flow pattern, the difference was mainly due to the distribution of soil CO₂. The difference in weathering fluxes was 1 %–12 % and 1 %–5 % between HF cases and VF cases, respectively, much smaller than differences between the PF cases. Comparing the HF and VF cases, differences in weathering fluxes were 73 % at 3.7 × 10⁻² m/a and 721 % at 3.7 m/a, which is about 1–2 orders of magnitude higher than differences induced by soil CO₂ distribution. Between the Grass_PF and Wood_PF cases, the differences were in the range of 17 %–207 % at the flow range of 3.7 × 10⁻²–3.7 m/a. The DIC fluxes showed similar trends.

Development of reaction fronts at the century scale. To explore the longer-term effects, we ran the PF cases at 0.37 m/a for 300 years. More water flowing vertically through abundant calcite zones in Wood_PF resulted in faster weathering and a deeper reaction front at a depth of ∼ 210 cm compared to ∼ 95 cm in Grass_PF (Fig. 5a). The depletion of calcite led to an increase in porosity (Fig. 5b), which was over 1 order of magnitude higher at the domain scale of the woodland than that in the grassland (Fig. 5c). Permeability evolution had a similar trend to porosity (not shown here). This indicates that if deeper roots promoted more water into deeper soils, they would push reaction fronts deeper and control the position where chemically unweathered bedrock was transformed into weathered bedrock. At timescales longer than century scale, calcite may become depleted, which ultimately reduces weathering rates (White and Brantley, 2003) and lead to similar weathering fronts in grasslands and woodlands.

Figure 5 Predicted soil profiles of (a) calcite volume change (calcite = initial calcite volume – current calcite volume) and (b) porosity after 300 years; (c) predicted temporal evolution of domain-scale porosity in the grassland and woodland. The infiltration rate is 0.37 m/a.

4.3 The regulation of weathering rates by CO₂–carbonate contact time

Natural systems are characterized by preferential flow paths such that flow distribution is not uniform in space. Weathering in such systems with preferential flow in zones of differing reactivities has been shown to hinge on the contact time between water and the reacting minerals instead of all minerals that are present (Wen and Li, 2018) (Eq. 1). Here we contextualize the weathering rates in different scenarios (symbols) with predictions from an upscaled rate law developed by Wen and Li (2018) that incorporates the effects of heterogeneities in flow paths (Fig. 6). The time ratio $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ compares the domain water contact time (or residence time) with the contact time with dissolving calcite. Note that τ_a is total domain pore volume V_T∕total water flow Q_T, and τ_ad,r is total reactive pore volume V_r∕total water fluxes passing through reactive zone Q_r. Also note that water passing through the reactive zone stays in the subsurface longer, such that it is older water in general. The ratio $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ is therefore akin to the fraction of older water compared to the total water fraction. The older water fraction F_ow is the counterpart of the young water fraction F_yw discussed in the literature (Kirchner, 2016, 2019). In Grass_VF and Wood_VF (open circles) where all CO₂-charged water flew through the deeper calcite zones, CO₂–calcite interactions reached maximum such that values of $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ approached 1.0, meaning all water interacted with calcite. Under this condition, weathering rates were the highest among all cases. In contrast, in Wood_HF and Grass_HF (open diamonds) where all CO₂-charged water bypassed the deeper calcite zone, $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ can be 1–3 orders of magnitude lower, and weathering rates were at their minima. The $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ value in the water partitioning cases (Grass_PF and Wood_PF, filled circles) fell in between. At the same τ_a, the Wood_PF case with deepening roots promoted CO₂–water–calcite contact (i.e., larger $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$) and dissolved calcite at higher rates.

Figure 6 Calcite weathering rate as a function of the reactive time ratio $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$, a proxy of the fraction of the older, reactive water compared to the total water fluxes. Symbols are rates from numerical experiments. Gray lines are predictions from the rate law equation at α=0.8 (Wen and Li, 2018). Large to small dots and thick to thin lines are for infiltration rates from 10^0.6 to 10^−1.4 m/a. The red filled stars represent the monthly rates in September (highest soil CO₂) and March (lowest soil CO₂) in the base case at 10^−0.4 (0.37) m/a; the black to gray dashed lines represent predictions with increasing soil CO₂ level (i.e., larger τ_eq) from Eq. (1). The gray filled star is for the case without soil CO₂. At any specific infiltration rate or τ_a, the VF and HF cases bracket the two ends, whereas the PF cases fall in between. The Wood_PF case with deeper roots enhanced the CO₂–water–calcite contact (i.e., larger $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ ratios), leading to higher calcite weathering rates compared to grasslands. The differences of weathering rates induced by different soil CO₂ level were relatively small compared to those hydrological changes induced by rooting depth.

The magnitude of the rate difference also depends on the overall flow rates (or domain contact time τ_a). At fast flow with small τ_a (large symbols and thick lines in Fig. 6), flow partitioning has a larger influence. At τ_a of 0.47 years, the weathering rates in the VF cases were more than an order of magnitude higher than those in the PF cases. The weathering rate in the woodland was over 7 times that of the grassland. In contrast, at τ_a of 47 years, the rate differences between grass and woody cases were less than 25 %, largely because the dissolution has reached equilibrium. In the VF and HF cases where flow conditions were the same, the soil CO₂ distribution in grassland and woodland differed only slightly, leading to similar values of $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ and weathering rates and indicating minimal impacts of soil CO₂ distribution. These rates from numerical experiments closely follow the prediction from the upscaled reaction rate law (gray lines, Eq. 1). The rate law predicted that weathering rates increased from HF cases with small $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ values to VF cases where $\frac{{\mathit{\tau }}_{\text{a}}}{{\mathit{\tau }}_{\text{ad,r}}}$ approached 1. It also showed that weathering rates reached their maxima when all water flows through the reactive zone, i.e., when the fraction of older, reactive water is essentially 1.0.

Figure 7 (a) Ca and (b) DIC data from the literature (gray dots) and prediction lines of Eqs. (2)–(3) at 10 ^∘C (solid line) and 17 ^∘C (dashed line). Measured spring water concentrations are from carbonate-dominant catchments in the literature: Abongwa and Atekwana (2015), Lopez-Chicano et al. (2001), Moral et al. (2008), Huang et al. (2015), Ozkul et al. (2010), Dandurand et al. (1982), Calmels et al. (2014), Kanduc et al. (2012), and Tsypin and Macpherson (2012).

5 Discussion

Because carbonate dissolution is thermodynamically controlled and transport limited, the overall weathering rates depend on how much CO₂-charged water flushes through the carbonate zone. This work indicates that deepening roots potentially enhance weathering rates in two ways. First, roots can control thermodynamic limits of carbonate dissolution by regulating how much CO₂ is transported downward and enters the carbonate-rich zone. In fact, the base-case grassland data and model reveal that the concentrations of Ca and DIC are regulated by seasonal fluctuation of pCO₂ and soil respiration. Second, roots can control how much water fluxes through the carbonate zone and export reaction products at equilibrium such that more dissolution can occur. The numerical experiments indicate that carbonate weathering at depth hinges on the recharge rate of delivery of CO₂-enriched water. Deepening roots in woodlands that channel more water into unweathered carbonate at depths can elevate weathering rates by more than an order of magnitude compared to grasslands. Below we elaborate and discuss these two messages.

5.1 The thermodynamics of carbonate weathering: control of temperature and pCO₂

The base-case data and simulation showed that calcite dissolution reaches equilibrium rapidly and is thermodynamically controlled (Sect. 4.1), which is generally well known and echoes observations from literature (Tsypin and Macpherson, 2012; Gaillardet et al., 2019). The extent of dissolution, or solubility indicated in Ca and DIC concentrations, is determined by soil CO₂ levels that in turn hinge on ecosystem functioning and climate. In hot, dry summer, soil respiration rates reach maxima in upper soil horizons and pCO₂ peaks (Fig. 2). In wet and cold winter, soil pCO₂ plummets, leading to much lower Ca and DIC concentrations. Based on Reactions (0–4) (Table 1) and temperature dependence of equilibrium constants, the following equations can be derived (detailed derivation in the Supplement) for Ca and DIC concentrations in carbonate-dominated landscapes:

$$\begin{array}{}\text{(2)}& {\displaystyle }{C}_{\text{Ca}}=[K_{\text{t,25}}\times e^{f\left(T\right)}\times p{\text{CO}}_{\mathrm{2}}{]}^{\frac{\mathrm{1}}{\mathrm{3}}}\text{(3)}& {\displaystyle }\begin{array}{rl}& C_{\text{DIC}}=\mathrm{2}[K_{\text{t,25}}\times e^{f\left(T\right)}\times p{\text{CO}}_{\mathrm{2}}{]}^{\frac{\mathrm{1}}{\mathrm{3}}}+K_{\mathrm{1},\mathrm{25}}\times e^{f\left(T\right)}\times p{\text{CO}}_{\mathrm{2}}\\ & \text{where}f\left(T\right)=-{\displaystyle \frac{\mathrm{\Delta }{H}_{\mathrm{1}}^{{}^{\circ }}}{\mathrm{4}R}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{298.15}}}\right)\end{array}\end{array}$$

Here C_Ca and C_DIC are concentrations of Ca and DIC (mol/L); pCO₂ is in units of atm; K_t,25 is the total equilibrium constant $K_t(=\frac{a_{{\text{Ca}}^{\mathrm{2}+}}{a}_{{\text{HCO}}_{\mathrm{3}}^{-}}^{\mathrm{2}}}{p{\text{CO}}_{\mathrm{2}}}=K_{\mathrm{1}}{K}_{\mathrm{4}}$) of the combined Reactions (1) and (4) at 25 ^∘C; $\mathrm{\Delta }{H}_{\text{t}}^{{}^{\circ }}$ is the corresponding standard enthalpy (−35.83 kJ/mol); K_1,25 and $\mathrm{\Delta }{H}_{\mathrm{1}}^{{}^{\circ }}$ (−19.98 kJ/mol) are the equilibrium constant and standard enthalpy of Reaction (1) (in Table 1); R is the gas constant (= 8.314 × 10⁻³ kJ/K/mol). Equations (2) and (3) imply that DIC and Ca in upper soil water (0.2 m) are lower compared to groundwater (3.6 m) in the base case at Konza, due to lower dissolved CO₂(aq) at higher temperature and higher diffusion rates in upper soil.

Equations (2) and (3) were tested with soil pCO₂ and spring water chemistry data from eight carbonate-dominated catchments (Dandurand et al., 1982; Lopez-Chicano et al., 2001; Moral et al., 2008; Ozkul et al., 2010; Kanduc et al., 2012; Calmels et al., 2014; Abongwa and Atekwana, 2015; Huang et al., 2015) (also see the Supplement for details). As shown in Fig. 7, spring water (representing groundwater) DIC and Ca concentrations increase with pCO₂. Equations (2)–(3) can describe these data and DIC and Ca levels from numerical simulations from this work (empty dots in Fig. 7). The lines of Eqs. (2)–(3) describe the relationship at 10 ^∘C, the mean annual temperature, confirming the thermodynamics control of carbonate dissolution by soil CO₂ levels. The equation lines closely predicted Ca and DIC under high-pCO₂ and lower-pH conditions (<8.0), because these conditions ensure the validity of the assumption of negligible CO${}_{\mathrm{3}}^{\mathrm{2}-}$ in the derivation of the equations. The presence of cations and anions other than Ca and DIC can complicate the solution and can bring significant variations of DIC and Ca concentrations under the same pCO₂ conditions.

High temperature (T=17 ^∘C) leads to lower DIC and Ca concentrations by about 10 % due to the lower calcite and CO₂ solubility at higher T. Higher T also elevates pCO₂ by enhancing soil respiration. Various equations have been developed in literature to predict soil pCO₂ based on climate and ecosystem functioning indicators such as net primary production (NPP) (Cerling, 1984; Goddéris et al., 2010; Romero-Mujalli et al., 2019a). These equations can be used together with Eqs. (2–3) for the estimation of Ca and DIC concentrations in carbonate-derived waters. Gaillardet et al. (2019) showed that pCO₂ can increase by 2 times with T increasing from 9 to 17 ^∘C, which can elevate DIC and Ca concentrations over 50 %. Macpherson et al. (2008) observed a 20 % increase in groundwater pCO₂ in Konza over a 15-year period and suggested that increased soil respiration under warmer climate may have elevated soil and groundwater pCO₂. Zhi et al. (2020) showed that the stream concentrations of dissolved organic carbon, a product of soil respiration, tripled in warmer years as minimum average temperature approaches 0 ^∘C in a high elevation mountain. Hasenmueller et al. (2015) demonstrated topographic controls on soil CO₂. Soil CO₂ and Ca and DIC concentrations are an integrated outcome of climate, soil respiration, subsurface structures, and hydrological conditions.

5.2 Hydrological controls of root-enhanced carbonate weathering

Evidence from field data. Data in Konza show that although soil pCO₂ peaks in summer, these peaks do not occur right away in deeper groundwater until about 2 months later. Macpherson et al. (2008) and Tsypin and Macpherson (2012) contributed this lag time to the water travel time from the soil to the groundwater aquifer. The calculation of travel time based on depth difference in soil and groundwater sampling location (214 cm) and average velocity (0.37 m/a) indicates that it will take 7–8 months on average for water to reach deeper groundwater. The 2-month delay, much shorter than the estimated 7–8 months, suggests that CO₂-charged water may arrive at deeper zones via preferential flow facilitated by roots or other types of macropores such as large conduits that are often observed in carbonate formations (Hartmann et al., 2014; Husic et al., 2019).

The numerical experiments suggest that the root-relevant hydrology can play an essential role in enhancing chemical weathering (Figs. 4 and 6). In general, the hydrological enhancement of weathering rates also alludes to the key connection between weathering and the transit (travel) time distribution (McGuire and McDonnell, 2006; Sprenger et al., 2019). In particular, water that routes through lateral paths is younger; water that penetrates deeper is older. Figure 6 says that weathering rates increase with increasing reactive water fraction, until reaching their maxima when all water is in contact with reactive minerals. This aligns with observations at Konza that woody-encroached watersheds exhibit higher Ca fluxes in streams and supports the hypothesis that deeper roots can enhance mineral–water interaction via deeper flow paths (Sullivan et al., 2019a). Deepening roots can also enhance connectivity between shallow and deeper zones, therefore reducing concentration contrasts between soil water and groundwater. This can lead to more chemostatic C–Q relationships as shown in Wood_PF ($b=-\mathrm{0.14}$ in C=aQ^b) compared to Grass_PF ($b=-\mathrm{0.33}$) (Fig. 4). These findings echo conclusions from Zhi et al. (2019) and Zhi and Li (2020) that chemistry differences in shallow versus deeper waters regulate C–Q patterns. The b values from the C–Q relationships coming out of the 1-D modeling (Fig. 4) suggest that if flow partitioning is the only difference between the grassland and woody watersheds, a C–Q relationship exhibiting dilution with negative b values is expected in the grassland. The Konza stream data in Sullivan et al. (2019a) however showed that C–Q slopes in grasslands ($b=-\mathrm{0.003}$) and in woody-encroached lands (b=0.013) are both close to zero (Figs. 7 and 8 in Sullivan et al., 2019a). The root influence on the C–Q relationship therefore remains equivocal.

Limitations of the model. This discrepancy may suggest that catchment features that are not represented in the simple 1-D model can influence C–Q relationships. The model does not explicitly simulate how and to what degree root distribution at depth alters flow pathways. Instead we focus on the first-order principles of the hydrological ramification of roots. The numerical experiments took general observations of rooting characteristics in grasslands and woodlands and assumed that deepening roots in woodlands promote higher flow partitioning into the deep subsurface (Canadell et al., 1996; Nardini et al., 2016; Pawlik et al., 2016). In natural systems, however, other factors can also influence flow partitioning. For example, contrasts in flow-conducting properties (i.e., porosity and permeability) in shallow and deep zones, physical and chemical heterogeneity in carbonate distribution (Wen and Li, 2018; Salehikhoo and Li, 2015), connectivity between different areas of the catchment in dry and wet times (Wen et al., 2020), and water table and corresponding lateral flow depth associated with the rainfall frequency and intensity (Li et al., 2017; Harman and Cosans, 2019). All these factors may affect the water–calcite contact, leading to changes in weathering rates (Fig. 6). In addition to the alterations of hydrological flow paths, deepening roots may also affect the proportion of plant water uptake or soil water loss through transpiration (Pierret et al., 2016; Zhu et al., 2018), further modifying recharge into the deep subsurface layers. For example, studies show that woody species in semiarid places predominantly use deep subsurface water while grasses predominantly use soil water from the upper soil layer (Ward et al., 2013). Furthermore, shallower root distributions in grasslands may be more efficient in using water from small rainfall events than forests with their deeper root distributions (Mazzacavallo and Kulmatiski, 2015). Other studies have also documented significant competition for water uptake in upper soil layers among both woody and grass species (Scholes and Archer, 1997). It remains inconclusive how these water uptake characteristics are best represented in numerical experiments. These processes are therefore not included at this point.

In addition, distributions of microbes and organic acids associated with rooting structures vary with sites and seasons and root channels. Dry conditions may trigger calcite reprecipitation. Microbial activities surrounding living and dead roots can also lead to calcite precipitation and infilling of fractures and other macropores, as well as alter flow pathways (Lambers et al., 2009). Organic acids can decrease soil water pH, increase mineral solubilities through organic–metal complexations, and accelerate chemical weathering (Pittman and Lewan, 2012; Lawrence et al., 2014). Their impacts on carbonate weathering kinetics might be smaller (due to the fast kinetics) compared to their alteration of calcite solubility in natural systems. Such processes further complicate how to represent biogeochemical processes in models. Although we do not explicitly simulate these competing processes, the model was constrained by the field data in grasslands and woodland that have already integrated these effects in natural systems. Indeed, the weathering rates can be understood as the net weathering rates that are the net difference between dissolution and reprecipitation. This is reflected in lower carbonate weathering rates under low infiltration (low flow) conditions.

The need for root characterization. The data–model discrepancy highlights the limitation of a simple model but also points to the need for measurements. In fact, because carbonate weathering is transport limited and depends largely on water flow via carbonate-rich zones, colocated measurements of rooting characteristics, flow, and water chemistry at depth are essential. Root measurements however rarely go deeper than 30 cm (Richter and Billings, 2015). Existing work exploring root influence on weathering focused primarily on effects of soil CO₂ and root exudates (Drever, 1994; Lawrence et al., 2014; Gaillardet et al., 2019). Although it is well known that roots play a paramount role in modulating macropores and subsurface flow (Fan et al., 2017), the interactions between root characteristics, flow partitioning, and chemical weathering have remained poorly understood. Rooting characteristics depend on climate, plant species, topography, soil properties, and geology (Canadell et al., 1996; Mazvimavi et al., 2004; Price, 2011; Nardini et al., 2016). Further studies are needed to characterize root distribution beyond 30 cm, how they vary with intrinsic plant species and external conditions, and how and to what extent they alter subsurface flow. The first step could link rooting characteristics, including density and depth, to soil properties and borrow insights from existing relationships between soil properties and subsurface structure. For example, images of roots and pore structures can be used to characterize the spatiotemporal heterogeneity of fluxes (Renard and Allard, 2013). Geostatistical indices such as permeability variance and correlation length can be used to quantify rooting structure and relate to flow partitioning and mineral weathering via numerical experiments (Wen and Li, 2017). The combination of numerical reactive transport experiments built on realistic rooting structure can help develop models for estimating the influence of rooting dynamics on water and carbon cycles at the catchment scale.

5.3 Deepening roots enhance carbon fluxes into the deep subsurface: a potential carbon sink?

Mounting evidence has shown that the terrestrial system has become a stronger carbon sink in recent decades (Heimann and Reichstein, 2008), potentially accounting for the missing carbon sink as large as ∼ 1 Pg C/a in the global carbon budget (Houghton, 2007; Cole et al., 2007). Although still much debated, recent studies have proposed the downward transport of soil-respired CO₂ and DIC into groundwater aquifers in deserts as a possible carbon sink in the global carbon cycle (Ma et al., 2014; Li et al., 2015). Considering the longer residence time of DIC in groundwater (10²–10⁴ years) than in the atmosphere (∼ 10¹–10² years; Archer and Brovkin, 2008), groundwater may act as a CO₂ storage sink. This work indicated that deepening roots can potentially reroute DIC fluxes to deeper groundwater storage. In particular, vegetation in dry places like deserts often has deep roots (Gupta et al., 2020). In fact, plants are known for growing deeper roots to tap groundwater during droughts (Brunner et al., 2015). Deepening roots can enhance downward water drainage (i.e., high vertical connectivity) to the depth and potentially facilitate the transport of DIC fluxes into the deep subsurface. As the pace of climate change accelerates, summer droughts are expected to intensify, which can potentially channel more CO₂ into the deeper subsurface via deepening roots.

With constraints from soil CO₂ data, the simulated CO₂ production rates in Grass_PF and Wood_PF are similarly at ∼ 0.6 mol C/m²/a under the infiltration rate of 0.37 m/a (Fig. S2). This is at the low end of the belowground net production (NPP) estimations at the Konza site (∼ 2.0–30.0 mol C/m²/a), assuming that belowground NPP accounts for 50 % of the total NPP (Lett et al., 2004; Knapp and Ojima, 2014). The simulations showed that in woodlands the DIC downward fluxes can be > 2.0 times higher than those in grasslands (Fig. 4). In other words, more soil-respired CO₂ can transport to groundwater and become stored there for centuries to millennia before entering a stream. At the short timescale, this enhanced downward transport will reduce CO₂ escape back into the atmosphere. This may explain the observations at Konza that woody encroachment increased NPP; however, soil CO₂ flux was significantly reduced compared to the open grassland (Lett et al., 2004). Changing land cover (e.g., woody encroachment and boreal forest creep) however is not the only mechanism for a carbon sink (Stevens et al., 2017; Wang et al., 2020). Older aged forests also tend to have deepening roots and may act as the carbon sink (Luyssaert et al., 2008), although this is not the case in Konza.

6 Conclusions

This work aims to understand thermodynamic and hydrological control of carbonate weathering driven by deepening roots. Field data and reactive transport simulation for a grassland in the Konza Prairie LTER site suggest that carbonate dissolution is thermodynamically controlled, and seasonal changes in temperature and pCO₂ drive variations of Ca and DIC concentrations in soil water and groundwater. We derived equations based on reaction thermodynamics (Eqs. 2–3) to estimate Ca and DIC as a function of pCO₂ and temperature, which have been shown to be applicable in other carbonate-dominated systems. The numerical experiments probed the potential effects of deepening roots on weathering by channeling a higher proportion of vertically downward flow (40 % of the total) into the deep subsurface with abundant calcite. The results show that deeper penetration of roots and higher vertical flow (recharge) enhanced CO₂–carbonate contact. At an infiltration rate of 3.7 m/a, calcite weathering flux in woodlands was about 200 % higher than that in grasslands. At 3.7 × 10⁻² m/a, the weathering flux in woodland was 17 % higher. The hydrological impacts on carbonate weathering were much higher than the biogeochemical impacts via elevated soil CO₂ alone, underscoring the importance of rooting depth in weathering. The modeling demonstrates that weathering rates depend on flow partitioning (the older water fraction that penetrates deeper) and relative magnitude of water contact time in the deep carbonate zone. At the century scale, with the higher proportion of vertical flow, the deeper roots pushed the weathering fronts 2 times deeper and resulted in a 10-times greater increase in porosity and permeability. Broadly, this deeper propagation of reaction fronts may accelerate rates of channel incision and hillslope erosion (Lebedeva and Brantley, 2013; Brantley et al., 2017b) and therefore speed up landscape evolution (Phillips, 2005). It alludes to the importance of considering changes in subsurface hydrological flows associated with shifts in vegetation types and rooting characteristics. This is particularly relevant as we assess the effects of climate change, land cover, and elevated atmosphere CO₂ concentrations on chemical weathering and carbon cycling.