Deepening roots can enhance carbonate weathering by amplifying CO2-rich recharge

Carbonate weathering is essential in regulating atmospheric CO2 and carbon cycle at the century timescale. Plant roots accelerate weathering by elevating soil CO2 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 CO2 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 pCO2 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 CO2. 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 CO2 data suggest relatively similar soil CO2 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−2 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 CO2-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 CO2–carbonate contact in the deep subsurface and (2) the hydrological impacts of rooting characteristics can be more influential than those of soil CO2 distribution in modulating weathering rates. We call for colocated characterizations of roots, subsurface structure, and soil CO2 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. Published by Copernicus Publications on behalf of the European Geosciences Union. 56 H. Wen et al.: Deepening roots can enhance carbonate weathering

Abstract. Carbonate weathering is essential in regulating atmospheric CO 2 and carbon cycle at the century timescale. Plant roots accelerate weathering by elevating soil CO 2 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 2 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 2 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 2 . 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 2 data suggest relatively similar soil CO 2 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 −2 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 2 -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 2 -carbonate contact in the deep subsurface and (2) the hydrological impacts of rooting characteristics can be more influential than those of soil CO 2 distribution in modulating weathering rates. We call for colocated characterizations of roots, subsurface structure, and soil CO 2 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.

Introduction
Carbonate weathering has long been considered negligible as a long-term control of atmospheric CO 2 (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-ofmagnitude-higher carbon store in carbonate reservoirs compared to the atmosphere (Gaillardet et al., 2019;. 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 2 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 2 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 2 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 Jobbagy and Jackson, 2000), which may help elevate the deep CO 2 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 tex-ture 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 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 2 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 2 distribution are missing, owing to the limitation in data that detail rooting effects on flow partitioning and complex hydrologicalbiogeochemical interactions . 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 2 distribution? Which factor (flow partitioning or soil CO 2 distribution) predominantly controls weathering? We hypothesized that deepening roots in woodlands enhance carbonate weathering by promoting CO 2 -enriched recharge into the deep, carbonateabundant 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. . 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 3 ) × cross-section area (m 2 )) 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 T = Q L + Q G .
2 Research site and data sources

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 chloriteillite 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  and groundwater pCO 2 . The concentrationdischarge 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 .

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. Soil CO 2 production through CO 2 (g * ) and dissolving into CO 2 (aq) (0) CO 2 (g * ) ↔ CO 2 (g) --−9.00 1.0 (1) CO 2 (g)↔ CO 2 (aq) a Values of K eq were interpolated using the EQ3/6 database (Wolery et al., 1990), except Reactions (1) and (4) (i.e., K 1 and K 4 ). b CO 2 (aq) concentrations were calculated through C CO 2 (aq) = K 1 pCO 2 . The prescribed C CO 2 (aq) values were used as equilibrium constants in CrunchTope to describe how much soil CO 2 was available for weathering. c Calcite in the upper soil (above horizon B) is mixed with other minerals and has small particle size , 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 2 (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.
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.

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

Reactions
In the model, the upper soil layers have more anorthite (CaAl 2 Si 2 O 8 (s)) and K-feldspar (KAlSi 3 O 8 (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 2 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 repro-duce field data. The kinetics follows the transition state theory (TST) rate law (Plummer et al., 1978) r = kA 1 − IAP K eq , where k is the kinetic rate constant (mol/m 2 /s), A is the mineral surface area per unit volume (m 2 /m 3 ), 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 4 ) was calibrated by fitting field data of Ca and alkalinity. To reproduce the observed soil CO 2 profile, the soil CO 2 production rate (mol/m 2 /a) at the domain scale (calculated by the mass change in the solid phase CO 2 (g * ) over time) was assumed to increase with infiltration rates (Fig. S2). This is consistent with field observations that soil CO 2 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 2 production rates increasing from 3.2 to 63.0 mol/m 2 /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 2 /a for soil CO 2 production. The simulated soil CO 2 production rates across different infiltration rates here (∼ 0.1-10 mol/m 2 /a shown in Fig. S2) were close to the reported belowground net primary production (belowground NPP) of typical ecosystems: ∼ 0.8-100 mol/m 2 /a for grasslands (Gill et al., 2002) and ∼ 0.4 -40 mol/m 2 /a for woodlands (Aragão et al., 2009).

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 L + Q G ). The flow field was implemented in Crunch-Tope using the "PUMP" option.

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.

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 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 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 Eberbach, 2003;Billings et al., 2018).
The experiments aimed to compare the general, averaged behaviors rather than event-scale dynamics so the annualaverage soil CO 2 data and corresponding prescribed CO 2 (aq) concentrations were used ( Table 2). The experiments focused on calcite weathering (Reactions 0-4) and excluded silicate weathering reactions (Reactions 5-7). The mineraldissolution 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 2 generation) influences of deepening roots. Though other potential differences might be induced by deepening roots (e.g., water uptake, water table, and transpiration)  and influence weathering, we assumed they remain constant across the grassland and woodland simulations to examine the relative influences of flow path vs. CO 2 generation on carbonate weathering.

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 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 −9 and 10 −5 mol/m 2 /s in both grasslands and woodlands (Bengtson and Bengtsson, 2007;Ahrens et al., 2015;Carey et al., 2016). Soil CO 2 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 2 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 2 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 2 level. Below we describe details of the numerical experiments (Table 3) exploring the influence of hydrological versus biogeochemical impacts of roots.

Three numerical scenarios
Each scenario in Table 3 includes a grassland and woodland case, with their respective profiles of calcite distribution and soil pCO 2 kept the same (columns 3 and 4 in Table 3). The only difference in different scenarios is the flow partitioning (column 5). Soil pCO 2 in the grassland (red line in column 4) was set to reflect the annual average of the Konza site. Soil pCO 2 (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 2 (g * ) producing CO 2 gas and CO 2 (aq) but differed in depth distributions. The CO 2 (g * ) depth distribution was constrained by CO 2 (g) field data ( Table 2). The distribution of CO 2 (g * ) in grasslands is slightly steeper than that in the woodland, with higher density of roots and more abundant CO 2 (g * ) in the upper soils. Below the rooting depth, CO 2 (g * ) was assumed to be smaller (by 10 times) to represent the potential soil CO 2 sources from microbial activities (Billings et al., 2018).
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 endmember 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 2 distribution versus hydrology differences. All scenarios were run under infiltration rates from 3.7 × 10 −2 to 3.7 m/a (10 −4 -10 −2 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 2 ) convert into rates in units of moles per year (mol/a).

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 (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.

Quantification of weathering rates and their dependence on CO 2 -carbonate contact
To quantify the overall weathering rates and CO 2 -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 charac-teristic 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 ): τ a = Lφ Q 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: 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 τ a τ 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 kA T 1 − exp − τ eq τ a and the heterogeneity factor 1 − exp −L τ a τ ad,r α that quantifies effects of preferential flow paths arising from heterogeneous distribution of minerals. When τ a τ ad,r > 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 2 - Figure 3. Simulated depth profiles of (A) calcite (vol %) and (B) soil CO 2 production rate; (C) CO 2 (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 2 (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. carbonate interactions. In contrast, a small τ a τ ad,r ratio (< 1) reflects that water is replenished slowly in the reactive calcite zones, leading to less CO 2 -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 shortterm scale, with negligible changes in the solid phase.

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 2 (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 2 (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 2 (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 2 (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 2 at their corresponding depths, indicating the predominate control of soil CO 2 on carbonate weathering. The groundwater concentrations of these species were also higher than soil concentrations. The simulated groundwater DIC (approximately summation of CO 2 (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 2 , which determined the CO 2 (aq) level and its spatial variation, and equilibrium constant (K eq ) of calcite dissolution. Imposition of monthly variations and depth distributions of soil CO 2 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 . 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.

Numerical experiments: the significance of hydrological flow partitioning
Scenario 1 for hydro-biogeochemical effects with flow partitioning (Wood P F and Grass P F ). Figures 3A1-E1 show depth profiles of calcite and soil CO 2 production rates and steady-state concentrations of reaction products. The soil CO 2 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 2 (g) level (released from CO 2 (g * )) increased with depth due to lower efflux into the atmosphere in deeper zone and downward fluxes of CO 2charge 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 2 (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 2 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 2 -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 2 -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 2 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 2 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 2 -calcite contact with 100 % flow through the calcite zone, whereas the HF cases had minimum CO 2 -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 with-  Ca and DIC concentrations, and (cd) 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. out soil respiration (i.e., no CO 2 (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.
Although weathering rates generally increased with infiltration rates, the woodland increased more (4.6 × 10 −2 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 2 -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 2 . 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 −2 m/a and 721 % at 3.7 m/a, which is about 1-2 orders of magnitude higher than differences induced by soil CO 2 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 −2 -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 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 domainscale porosity in the grassland and woodland. The infiltration rate is 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.

The regulation of weathering rates by CO 2 -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 τ a τ 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 τ a τ 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 wa- Figure 6. Calcite weathering rate as a function of the reactive time ratio τ a τ 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 2 ) and March (lowest soil CO 2 ) in the base case at 10 −0.4 (0.37) m/a; the black to gray dashed lines represent predictions with increasing soil CO 2 level (i.e., larger τ eq ) from Eq. (1). The gray filled star is for the case without soil CO 2 . 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 2 -water-calcite contact (i.e., larger τ a τ ad,r ratios), leading to higher calcite weathering rates compared to grasslands. The differences of weathering rates induced by different soil CO 2 level were relatively small compared to those hydrological changes induced by rooting depth.
ter fraction F yw discussed in the literature (Kirchner, 2016(Kirchner, , 2019. In Grass VF and Wood VF (open circles) where all CO 2charged water flew through the deeper calcite zones, CO 2calcite interactions reached maximum such that values of τ a τ 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 2 -charged water bypassed the deeper calcite zone, τ a τ ad,r can be 1-3 orders of magnitude lower, and weathering rates were at their minima. The τ a τ 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 2 -water-calcite contact (i.e., larger τ a τ ad,r ) and dissolved calcite at higher rates.
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 2 distribution in grassland and woodland differed only slightly, leading to similar values of τ a τ ad,r and weathering rates and indicating minimal impacts of soil CO 2 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 τ a τ ad,r values to VF cases where τ a τ 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.

Discussion
Because carbonate dissolution is thermodynamically controlled and transport limited, the overall weathering rates depend on how much CO 2 -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 2 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 2 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 2enriched 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.

The thermodynamics of carbonate weathering:
control of temperature and pCO 2 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 2 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 2 peaks ( Fig. 2). In wet and cold winter, soil pCO 2 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: (2) Here  Table 1); R is the gas constant (= 8.314 × 10 −3 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 2 (aq) at higher temperature and higher diffusion rates in upper soil.
Equations (2) and (3) were tested with soil pCO 2 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 2 . 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 2 levels. The equation lines closely predicted Ca and DIC under high-pCO 2 and lower-pH conditions (< 8.0), because these conditions ensure the validity of the assumption of negligible CO 2− 3 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 2 conditions. High temperature (T = 17 • C) leads to lower DIC and Ca concentrations by about 10 % due to the lower calcite and CO 2 solubility at higher T . Higher T also elevates pCO 2 by enhancing soil respiration. Various equations have been developed in literature to predict soil pCO 2 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 2 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 2 in Konza over a 15-year period and suggested that increased soil respiration under warmer climate may have elevated soil and groundwater pCO 2 .  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 2 . Soil CO 2 and Ca and DIC concentrations are an integrated outcome of climate, soil respiration, subsurface structures, and hydrological conditions. Evidence from field data. Data in Konza show that although soil pCO 2 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 2 -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 . 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 = −0.14 in C = aQ b ) compared to Grass PF (b = −0.33) (Fig. 4). These findings echo conclusions from Zhi et al. (2019) and  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  however showed that C-Q slopes in grasslands (b = −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 dis-tribution 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 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 , and water table and corresponding lateral flow depth associated with the rainfall frequency and intensity 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 organicmetal 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 2 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 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 2 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 2 -10 4 years) than in the atmosphere (∼ 10 1 -10 2 years; Archer and Brovkin, 2008), groundwater may act as a CO 2 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 2 into the deeper subsurface via deepening roots. With constraints from soil CO 2 data, the simulated CO 2 production rates in Grass PF and Wood PF are similarly at ∼ 0.6 mol C/m 2 /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 2 /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 2 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 2 escape back into the atmosphere. This may explain the observations at Konza that woody encroachment increased NPP; however, soil CO 2 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.

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 2 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 2 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 2 -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 −2 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 2 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 2 concentrations on chemical weathering and carbon cycling.
Data availability. All data used for model parameterization can be acquired from http://lter.konza.ksu.edu/data (last access: 22 November 2020, Macpherson, 2019). The input files necessary to reproduce the results are available from the authors upon request by emailing lili@engr.psu.edu.
Author contributions. HW, PS, and LL initiated the idea and designed the numerical experiments. GLM and SB provided the field data. HW ran the simulations, analyzed simulation results, and wrote the first draft of the manuscript. All coauthors participated in editing the manuscript.