Deepening roots can enhance carbonate weathering

Carbonate weathering is essential in regulating atmospheric CO2 and carbon cycle at the century time scale. Plant roots have been known to accelerate weathering by elevating soil CO2 via respiration. It however remains poorly understood how and 10 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 potentially enhance weathering in two ways. First, deepening roots can control thermodynamic limits of carbonate dissolution by regulating how much CO2 transports downward to the deeper carbonate-rich 15 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 fluctuation of 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, which has been shown to apply in multiple carbonate-dominated catchments. Second, numerical experiments show that roots control weathering rates by regulating the amount of water fluxes that flush through the carbonate zone and export reaction products at dissolution equilibrium. Numerical 20 experiments explored the potential effects of partitioning 40% of infiltrated water to depth in woodlands compared to 5% in grasslands. Soil CO2 data from woodand grasslands suggest relatively similar soil CO2 distribution over depth, and only led 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 207% as infiltration rates increased from 3.7×10-2 to 3.7 m/yr. Numerical experiments also indicated that weathering fronts in woodlands propagated > 2 times deeper compared to grasslands after 300 years at the 25 infiltration rate of 0.37 m/yr. These differences in weathering fronts are ultimately caused by the contact time of CO2-charged water with carbonate rocks. We recognize that modeling results are subject to limitations in representing processes and parameters, but we propose that the data and numerical experiments allude to the hypothesis that 1) deepening roots can enhance carbonate weathering; 2) the hydrological impacts of rooting characteristics can be more influential than those of soil CO2 distribution in modulating weathering rates. We call for co-located characterizations of roots, subsurface structure, soil CO2 levels, and their 30 linkage to water and water chemistry. These measurements will be essential to improve models and illuminate feedback mechanisms of land cover changes, chemical weathering, global carbon cycle, and climate.

however have underscored its significance in controlling global carbon cycle at the century time scale that is relevant to modern climate change, owing to its rapid dissolution, fast response to perturbations, and the order of magnitude 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 40 soil CO2 concentrations (Covington et al., 2015) arising from vegetation type (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 CO2 level.
Plant roots have long been recognized as a dominant biotic-driver of chemical weathering and global carbon cycle (Berner, 1992;Beerling et al., 1998;Brantley et al., 2017a). The growth of forests has been documented to elevate soil pCO2 and amplify 45 Dissolved Inorganic Carbon (DIC) fluxes (Berner, 1997;Andrews and Schlesinger, 2001). Rooting structure can influence weathering in two ways. First, rooting systems (e.g., grass-, shrub-and woodlands) may affect the distribution of soil carbon (both organic and inorganic), microbe biomass, and soil respiration at depth (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 CO2 and acidity that determine carbonate solubility and weathering 50 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 >70% of 55 water fluxes through soils (Watson and Luxmoore, 1986). In grasslands, the lateral, dense spread of roots in top soils promote horizontal 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, shrublands and forests with deep and thick roots tend to have a high proportion of macropores and high connectivity to the deep subsurface Nardini et al., 2016), enhancing the drainage of 60 water to the depth (Pawlik et al., 2016).
It is generally known that rooting characteristics vary among plant species and are critical in regulating flow paths and storage (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 CO2 and organic acids (Drever, 1994;Lawrence et al., 2014;Gaillardet et al., 2019).
Systematic studies on coupled effects of flow partitioning and soil CO2 distribution are missing, owing to the limitation in data 65 that detail rooting effects on flow partitioning and complex hydrological-biogeochemical interactions. Here we ask the questions: how and to what degree do rooting characteristics influence carbonate weathering when considering both flow partitioning and soil CO2 distribution? Which factor (flow partitioning or soil CO2 distribution) predominantly controls weathering? We hypothesized that deepening roots in woodlands enhance carbonate weathering by promoting deeper recharge and CO2-carbonate contact at depth (Figure 1). 70 We tested the hypotheses by a series of numerical experiments integrating reactive transport modelling and 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 ( Figure S1). 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 impact of biogeochemical and hydrological drivers and bracketed the range 75 of their potential impacts on weathering, thus providing insights on the missing quantitative link between rooting structure, flow paths, and chemical weathering.

Research site and data sources
Site description. Details of the Konza site are in the Supporting Information (SI) and references therein. Here we provide brief information most relevant to this work. Konza is a mesic native grassland where experimentally manipulated, long-term burning 80 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 85 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. In parallel to these changes is a detectable decline in stream flow and an increase in weathering rates  and groundwater pCO2 (Macpherson et al., 2008). Stream water concentration-discharge (C-Q) patterns 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. (2019) 90 hypothesized that concentration-discharge relationships may be affected by woody species with deeper roots, which altered flow paths and mineral-water interactions ( Figure 1). We focus on the upland watershed N04d in Konza ( Figure 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 95 (lter.konza.ksu.edu/data). Wet chemistry deposition data was from the National Atmospheric Deposition Program NADP (nadp.sws.uiuc.edu). Monthly data of soil gases (16, 84 and 152 cm), soil water (17 and 152 cm), and groundwater (366 cm) from 2009-2010 were used for this work (Tsypin and Macpherson, 2012) (Figure 2A). The sampling points were about 30 m away from the stream. More information on field and laboratory methods were included in the SI. 100 https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. 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 promotes lateral macropore development and lateral water flow; the deeper roots in the woodland induces vertical macropore development that supports vertical flow into groundwater. The gray color gradient reflects the calcite abundance with more calcite in depth. ML and MG with a unit of mol/yr (= 105 flow rate Q (m/yr) × species concentration C (mol/m 3 ) × unit 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. Infiltration rate QT = QL + QG.
3 Reactive transport modeling 3.1 Base case: 1D reactive transport model for the Konza grassland (calibrated with field data) 110 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 115 measurements (Tsypin and Macpherson, 2012). Detailed set up of domain, soil mineralogy, initial condition, and precipitation chemistry, are in the SI.
prescribed !" ! ($%) values were used as the equilibrium constants of the coupled Reaction 0-1 in the form of ( ( * ) ↔ ( ( ), such that the soil CO2 at different depth were represented (Section 4.1). In the base case, the soil profile of !" ! ($%) values were updated monthly based on the monthly soil CO2 data (  b. By combining the temporal-and spatial-dependent soil pCO2 data, CO2(aq) concentrations were calculated through !" ! ($%) = ' ( . The prescribed !" ! ($%) values were used as equilibrium constants in CrunchTope to describe how much soil CO2 was available for weathering. 135 c. Calcite in the shallow soil (above Horizon B) is mixed with other minerals and has small particle size , relatively "impure" and therefore has a lower Keq value than those at depth with relatively "pure" calcite. The Keq 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 solid phase CO2(g*) dissolution (i.e., soil CO2 production rate constant) was from (Bengtson and Bengtsson, 140 2007;Ahrens et al., 2015;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 shallow soils in Konza. In the later numerical experiments, these reactions were not included so as to focus on carbonate weathering. 145 Monthly measured soil CO2 data for the Konza grassland (base case) were from (Tsypin and Macpherson, 2012); The annualaverage soil CO2 used in the grassland experiments was averaged from the corresponding monthly measurements. More information on measurements at the Konza grassland were detailed in the SI. The annual-average soil CO2 in the woodland experiments was from the Calhoun site in South Carolina where trees (forest) are dominant (Billings et al., 2018). 150 b. CO2(aq) were estimated using Henry's law !" ! ($%) = ' ( ; The temperature-dependent K1 was calculated following Eq. S1. The !" ! ($%) values at different soil depths were used to prescribe the available soil CO2 for chemical weathering. c. Soil temperature was estimated from the soil water and shallow groundwater temperature (Tsypin and Macpherson, 2012;Billings et al., 2018). 155 Reactions. The shallow soils have more anorthite ( ( ( . ( )) and K-feldspar ( , . ( )) and deeper subsurface contains more calcite (Table S1). Table 1 summarizes reactions and thermodynamic and kinetic parameters. Soil CO2 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 160 the transition state theory (TST) rate law (Plummer et al., 1978) = ;1 − Flow partitioning. Rain water entered soil columns at the annual infiltration rate of 0.37 m/yr, estimated based on the difference 165 between measured precipitation (0.88 m/yr) and evapotranspiration (0.51 m/yr). At 50 cm, a lateral flow QL (soil water) exited the soil column to the stream at 0.35 m/yr. The rest recharged the deeper domain beyond 50 cm (to the groundwater system) at 0.02 m/yr (~ 2% of precipitation) and became groundwater ( Figure 1A), a conservative value compared to 2-15% reported in another study (Steward et al., 2011). The groundwater flow QG eventually came out at 366.0 cm and was assumed to enter the stream as part of discharge (= QL + QG). The flow field was implemented in CrunchTope using the "PUMP" option. 170 Calibration. We used monthly measured alkalinity and Ca concentration data in different horizons (Horizon A, B and groundwater in Figure 2A) from 2009-2010 at the Konza grassland for model calibration. The monthly Nash-Sutcliffe efficiency (NSE), which 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. At Konza, calcite in the shallow soil (above Horizon B) is mixed with other minerals, has small particle size, and is considered "impure"  and therefore 175 has a lower Keq value than those at depth with relatively "pure" calcite. The Keq of "impure" calcite (K4) was calibrated by fitting field data of Ca and alkalinity.
Numerical experiments were set up to understand how and to what extent roots regulate weathering rates and solute transport in grasslands and in woodlands. The base case from Konza was used to represent grasslands, and the Calhoun site (the 180 Calhoun Critical Zone Observatory in South Carolina, USA) was used as a representative woody site (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 (QL) at the shallow subsurface ( Figure 1A). At Konza, > 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 185 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 (QG). At the Calhoun site, > 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 CO2 data and corresponding prescribed CO2(aq) concentrations were used ( Table 2). The experiments focused on 190 calcite weathering (Reaction 0-4) and excluded silicate weathering reactions (Reaction 5-7). The mineral-dissolution-related parameters from the base case were used for all experiments. We compared the relative significance of hydrological and biogeochemical effects (soil CO2 level and distribution) of rooting depths.
Hydrological and bigeochemical 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 195 land cover conditions. The ratios of lateral flow in shallow 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 shallow soils over the overall infiltration can vary between 50% and 95%. Deeper roots in woodlands can increase deep soil permeability by > one 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 200 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 QG/QT 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 ~ 10 -5 mol/m 2 /s for both grasslands and woodlands (Bengtson and Bengtsson, 205 2007;Ahrens et al., 2015;Carey et al., 2016). Soil CO2 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 CO2 levels are further complicated also by their dependence on topographic position, soil depth, and soil moisture, all of which determine the magnitude of microbial and root activities and CO2 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 CO2 level. Below we describe details of the numerical experiments 210 (Table 3) exploring the influence of hydrological versus biogeochemical impacts of roots. b. Calcite distribution (black line) in all cases increases with depth; c. The flow field was implemented in CrunchTope using the "PUMP" option; d. Annual-average soil CO2 used for constraining the grassland (red line) was from the Konza grassland (Tsypin and Macpherson, 2012). The annual-average soil CO2 in the woodland (blue line) was from the Calhoun forest (Billings et al., 2018). 225

215
Three scenarios. Each scenario in Table 3 includes a grassland and woodland case, with their respective profiles of calcite distribution and soil pCO2 kept the same (column 3 and 4 in Table 3). The only difference in different scenarios is the flow partitioning (column 5 in Table 3). Soil pCO2 in the grassland (red line in Table 3) were set to reflect the annual average of the Konza site. Soil pCO2 (blue line) in the woodland is 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 solid phase "CO2(g*)" producing CO2 gas 230 and CO2(aq) but differed in depth distributions. The CO2(g*) served as the source of soil CO2 and was constrained by CO2(g) field data ( Table 2). In grasslands, the modeled distribution of CO2(g*) was steep, with higher density of roots and more abundant CO2(g*) in the shallow soils; while the distribution was less steep in woodlands . Below the rooting depth, CO2(g*) was assumed to be smaller (by 10 times) to represent the potential soil CO2 sources from microbial activities (Billings et al., 2018). 235 Scenario 1 considered flow partitioning (Table 3). With the large permeability contrast of soil and bedrock (over four orders of magnitude) in Konza (Macpherson, 1996), the GrassPF case (with QL = 95%×QT and QG = 5%×QT) represents an endmember case for the grassland. The woodland (WoodPF) case was set to have 60% lateral flow and 40% vertical flow into 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 240 bottom outlet (VF; WoodVF and GrassVF); Scenarios 3 had two cases with 100% horizontal flow (HF; WoodHF and GrassHF) 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 CO2 distribution versus hydrology differences. All scenarios were run under infiltration rates from 3.7×10 -2 to 3.7 m/yr (10 -4 -10 -2 m/day), the observed daily variation range at Konza. This was to explore the role of flow regimes and identify conditions where most and least significant 245 differences occur. To reproduce the observed soil CO2 profile, the soil CO2 production rate (mol/m 2 /yr) at the domain scale (calculated by the mass change of the solid phase CO2(g*) over time) was assumed to increase with infiltration rates ( Figure S2). https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. This is consistent with field observations that soil CO2 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 CO2 production rates increasing from 3.2 to 63.0 mol/m 2 /yr when the annual precipitation increased 250 from 400 to 1,200 mm. In addition, Wu et al. (2011) showed that increasing precipitation from 5 to 2,148 mm enhanced soil respiration by 40% and that a global increase in 2 mm precipitation per decade may lead to an increase of 3.8 mol/m 2 /yr for soil CO2 production. The simulated soil CO2 production rates across different infiltration rates here (~0.1-10 mol/m 2 /yr shown in Figure   S2) were close to the reported belowground net primary production (belowground NPP) of typical ecosystems: ~0.8-100 mol/m 2 /yr for grasslands (Gill et al., 2002) and ~ 0.4 -40 mol/m 2 /yr for woodlands (Aragao et al., 2009). 255 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 yrs, depending on infiltration rates. The lateral flux (soil water, ML=QL×CL) and vertical fluxes (groundwater, MG=QG×CG) were calculated at 50 and 366 cm, in addition to total fluxes (weathering rates, MT). These fluxes multiplied with unit cross-section area (m 2 ) convert into rates in units of mol/yr.
Carbonate weathering over century time scales: soil property evolution. To compare the propagation of weathering fronts 260 over longer time scales, we carried out two 300-yr simulations for Scenario 1 with flow partitioning under the base-case infiltration rate of 0.37 m/yr (i.e., GrassPF and WoodPF 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 7 and ∅ 7 is permeability and porosity in grid i at time t, and 7,9 and 265 ∅ 7,9 are the initial permeability and porosity, respectively.

Quantification of weathering rates and their dependence on CO2-carbonate contact
To quantify the overall weathering rates and CO2-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 characteristics times are important. The equilibrium time teq represents the characteristic timescale of . The reactive transport time $=,> quantified 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, AT is the total surface area, L is domain length, is geostatistical characteristics of spatial heterogeneity, and  Table S2. 285

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 (Figure 2). Note that y axis is inverted to display shallow soils at the top and deep soils at the bottom to be consistent with their subsurface position shown in Figure 2A. The measured CO2(g) varied between 0.24 -7.30% ( Figure 2B right 290 axis), one to two orders of magnitude higher than the atmospheric level of 0.04%. The estimated CO2(aq) ( Figure 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 CO2(aq) reached maxima in summer in soil horizons A and B and decreased to < 0.5 mM in winter and spring. The groundwater CO2(aq) exhibited a delayed peak in September and October, and dampened seasonal variation compared to soil horizons. The temporal trends of alkalinity, DIC, and Ca in groundwater mirrored those of soil CO2 at 295 their corresponding depths, indicating the predominate control of soil CO2 on carbonate weathering. The groundwater concentrations of these species were also higher than soil concentrations. The simulated groundwater DIC (approximately summation of CO2(aq) and alkalinity) was > 6 times higher than that in shallow 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 with anorthite as the dominating dissolving mineral (0 -1.8 m); their effects are negligible in Horizon 300 B and groundwater where fast-dissolving calcite has reached equilibrium. corresponding sampling depth of monthly field measurements (dots), including soil water at Horizons A (16 cm) and B 305 (152 cm), and groundwater (366 cm). The lines of CO2(aq) and CO2(g) in Figure B overlapped. Note that there are no DIC data so no dots in Figure D. The temporal trends of alkalinity, DIC, and Ca mirrored those of soil CO2, indicating its predominant control on weathering.
Several measurements/parameters were important in reproducing data, especially soil CO2 that determined the CO2(aq) level and its spatial variation, and equilibrium constant (Keq) of calcite dissolution. Imposition of monthly variations and depth 310 distributions of soil CO2 were essential to capture the variation of alkalinity and Ca data at different horizons. The imposition of calcite Keq was also critical for reproducing Ca concentrations. Impurities were suggested to affect Keq of natural calcite by a factor of ~2.0 . Calcite Keq in shallow 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 was fast and thermodynamically controlled. 315

Scenario 1 for hydro-biogeochemical effects with flow partitioning (WoodPF and GrassPF). Figures 3A1-E1 show
depth profiles of calcite and soil CO2 production rates and steady-state concentrations of reaction products. Although the model did not explicitly simulate soil respiration, the soil CO2 production rate was highest in the shallow soil at around 10 -6.5 mol/s and decreased to ~ 10 -9.5 mol/s at 366 cm ( Figure 3B1), consistent with the decline with soil depth observed in natural systems. The 320 CO2(g) level (released from the solid phase CO2(g*)) increased with depth due to autotrophic and heterotrophic respiration and downward fluxes of CO2-charge water from the shallow soil ( Figure 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 > 60 cm. The transition occurred between 35 cm -60 cm at the vicinity of 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 325 equilibrium, resulting in a short equilibrium distance. The equilibrated DIC and Ca concentrations below 60 cm followed the similar increasing trend of CO2(g) with depth in the deep zone ( Figure 3C1-E1).
Figures 3F1 show 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 330 much for groundwater. As expected, stream concentrations, 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.
https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License.  Figure  C are the annual average CO2(g) data (in Table 2). Arrows in Figure A indicate flow conditions. In F, soil water and 340 groundwater refer to concentrations at the lateral (50 cm) and vertical (366 cm) outlets, respectively. Higher fractions of vertical flow in WoodPF led to higher stream Ca concentration compared to GrassPF. In Scenarios 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.

345
Scenarios 2-3 for biogeochemical effects (without flow partitioning). Scenarios 2 and 3 were end-member cases that bracket the range of rooting effects. Calcite and soil CO2 were distributed the same way as their corresponding PF cases (Table 3). The WoodVF and GrassVF cases had 100% flow going downward via the deeper calcite zone maximizing the CO2-water-calcite contact. In the WoodHF and GrassHF cases (Table 3), all water exited at 50 cm, bypassing the deeper calcite-abundant zone, minimizing the CO2-water-calcite contact. Figures 3A2-F2 and 3A3-F3 show the VF (WoodVF and GrassVF) and HF (WoodHF and GrassHF) cases, 350 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 CO2 distribution ( Figure 3C2 and 3C3). Figure 3F2 indicates that the effluent Ca concentrations were slightly higher in WoodVF due to the high soil CO2 level at the bottom outlet. The GrassHF and WoodHF case almost had the same effluent Ca concentrations at the shallow soil ( Figure 3F3). 355 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 GrassPF and WoodPF cases fell in between ( Figure 4A-B). This is because the VF cases maximized the CO2-calcite contact with 100% flow through the calcite zone, whereas the HF cases had minimum CO2-calcite contact with water bypassing the calcite zone (column 3-4 in Table 3). The GrassPF and WoodPF cases allowed different extent of contact prescribed by the amount of flow via the calcite zone. A case run 360 without soil respiration (i.e., no CO2(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 emerge with almost no changes with infiltration rates. This is because in the VF cases, the concentrations mostly reached equilibrium concentration. In PF and HF cases, a large proportion of water flow through soils with negligible calcite where the water becomes further away from 365 equilibrium as infiltration rates increase. 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 370 hydrological impacts on weathering. The Ca fluxes in the VF cases (100% vertical flow) were higher than flow partitioning cases because they enabled maximum CO2-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 CO2. The difference of weathering fluxes was 1-12% and 1-5% between HF cases and VF cases, respectively, much smaller than differences between 375 the PF cases. Comparing the HF and VF cases, differences in weathering fluxes was 73% at 3.7 ×10 -2 m/yr and 721% at 3.7 m/yr, about 1-2 orders of magnitude higher than differences induced by soil CO2 distribution. Between the GrassPF and WoodPF cases, the differences were in the range of 17%-207% at the flow range of 3.7 ×10 -2 -3.7 m/yr. The DIC fluxes showed similar trends ( Figure 4D).
Although weathering rates generally increased with infiltration rates, the woodland increased more (4.6×10 -2 to 2.4 380 mol/yr in WoodPF). Estimates of the overall soil CO2 production rate in GrassPF and WoodPF suggest that the values under the same infiltration rate were all within a difference of 10% ( Figure S2), indicating that the difference of weathering rates are mainly driven by the flow partitioning. As a reference, the weathering rate in cases with soil respiration was up to 10 times higher compared to that without soil respiration.

Development of reaction fronts at the century scale.
To explore the longer-term effects, we ran the PF cases at 0.37 m/yr for 385 300 years. More water flowing vertically through abundant calcite zones in WoodPF resulted in faster weathering and deeper reaction front at a depth of ~210 cm compared to ~95 cm in GrassPF ( Figure 5A). The depletion of calcite led to an increase of porosity ( Figure 5B), which was over one order of magnitude higher at the domain scale of the woodland than that in the grassland ( Figure 5C). Permeability evolution had the similar trend with porosity (not shown here). This indicates that if deeper roots promoted more water into deeper soils, they would push reaction fronts deeper and controlled the position where chemically 390 unweathered bedrock was transformed into weathered bedrock. At time scales 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 CO2-carbonate contact
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 have been shown to hinge on the contact time 400 https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. 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 incorporate the effects of heterogeneities in flow paths ( Figure 6). The time ratio The magnitude of rate difference also depends on the overall flow rates (or domain contact time ta). At fast flow with small ta (large symbols and thick lines in Figure 6), flow partitioning has a larger influence. At ta of 0.47 year, 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 ta 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, 415 the soil CO2 distribution in grassland and woodland differed only slightly, leading to similar values of  (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/yr. The red filled stars represent the monthly rates in September (highest soil CO2) and March (lowest soil CO2) in the base case at 10 -0.4 (0.37) m/yr; the black to grey dash 425 lines represent predictions with increasing soil CO2 level (i.e., larger B% ) from Eq. (1). Grey filled star is for the case τ a τ ad,r 10 -3 10 -2 10 -1 10 0 10 1 10 -2 10 -1 https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. without soil CO2. At any specific infiltration rate or $ , the VF and HF cases bracket the two ends whereas the PF cases fall in between. The WoodPF case with deeper roots enhanced the CO2-water-calcite contact (i.e., larger D ( D (),* ratios), leading to higher calcite weathering rates compared to grasslands. The differences of weathering rates induced by different soil CO2 level were relatively small compared to those hydrological changes induced by rooting depth. 430

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

The thermodynamics of carbonate weathering: control of temperature and pCO2
The base case data and simulation showed that calcite dissolution reaches equilibrium rapidly and is thermodynamically controlled (Section 4.1), which echoes observations from other weathering studies (Tsypin and Macpherson, 2012;Gaillardet et al., 2019). The extent of dissolution, or solubility indicated in Ca and DIC concentrations, is determined by soil CO2 levels that are a function of ecosystem functioning and climate. In hot, dry summer, soil respiration reaches its maximum rates in shallow soil 445 horizons and pCO2 peaks (Figure 2). In wet and cold winter, soil pCO2 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 SI) for Ca and DIC concentrations in carbonate-dominated landscapes:  Table 1); R is the gas constant (=8.314×10 -3 kJ/K/mol). Eq. 2 and 3 imply that DIC and Ca in shallow soil water (0.2 m) are lower compared to groundwater (3.6 m) in the base case at Konza, due to lower dissolved CO2(aq) arising from higher temperature and higher diffusion rates in shallow soil. 455 Eq. 2 and 3 were tested with compiled soil pCO2 and spring water chemistry data from 8 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 SI for details). Plotting spring-water chemistry against pCO2 indicates that spring water (as representing groundwater) DIC and Ca concentrations increase with pCO2 (Figure 7), and Eq. 2-3 https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License.
can describe these data and also DIC and Ca levels from numerical simulations from this work (empty dots in Figure 7). High temperature (T = 17°C) leads to lower DIC and Ca concentrations by about 10% due to the lower calcite and CO2 solubility at higher T. Higher T also elevates pCO2 due to higher soil respiration. Various equations exist for predicting soil pCO2 based on climate and ecosystem functioning indicators such as Net Primary Production (NPP) (Cerling, 1984;Goddéris et al., 475 2010;Romero-Mujalli et al., 2019a). These equations can be used together with Eq. 2-3 for the estimation of Ca and DIC concentrations in carbonate-derived waters. Gaillardet et al. (2019) shows that pCO2 can increase by 2 times with T increase from 9°C to 17°C, which can elevate DIC and Ca concentrations over 50%. Macpherson et al. (2008) observed a 20% increase in groundwater pCO2 in Konza over a 15-year period and suggested that increased soil respiration due to climate warming may have elevated soil and groundwater pCO2. Hasenmueller et al. (2015) demonstrated topographic controls on soil CO2. Variations of soil 480 CO2 and Ca and DIC concentrations therefore are an integrated outcome of climate, soil respiration, subsurface structures, and hydrological conditions.

Hydrological controls of root-enhanced carbonate weathering
Evidence from field data. Data in Konza show that although soil pCO2 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 the delay 485 to the water travel time from soil to 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/year) indicate that it will takes 7-8 months for water to reach deeper groundwater. The 2 months delay, much shorter than the estimated 7-8 months, suggests that CO2-charged water may arrive deeper zones via preferential flow facilitated by roots or via large storm events with much faster flow. https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License.
The numerical experiments suggest that the root-relevant hydrology can play an essential role in enhancing chemical 490 weathering (Figure 4 and 6). The enhanced weathering rates aligns with observations at Konza that woody-encroached watersheds exhibit higher Ca fluxes in stream 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 WoodPF (b = -0.14 in C=aQ b ) compared to GrassPF (b = -0.33) (Figure 4) experiments took general observations of rooting characteristics in grasslands and woodlands, and assumed that deepening roots 505 in woodlands promote higher flow partitioning into deep subsurface Nardini et al., 2016;Pawlik et al., 2016).
In natural systems, however, other factors can also influence flow partitioning. This includes, for example, contrasts in flowconducting property (i.e., porosity and permeability) in shallow and deep zones, physical and chemical heterogeneity in carbonate distribution (Wen and Li, 2018), and connectivity between different areas of the catchment in dry and wet times (Wen et al., 2020).
Representing these dynamics require a large number of processes and parameters that we do not have data to constrain at this point. 510 The need for root-relevant measurements. 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, co-located 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 CO2 and root exudates (Drever, 1994;Lawrence et al., 2014;Gaillardet et al., 2019). 515 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 has 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 520 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 them 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 525 quantitative relationship between roots and flow partitioning, and could support models for estimating the influence of rooting dynamics on water and carbon cycles at the catchment scale.
Mounting evidence has shown that the terrestrial system has become a stronger carbon sink in recent decades (Heimann and Reichstein, 2008), potentially counting for the missing carbon sink as large as ~1 Pg C yr -1 in the global carbon budget 530 (Houghton, 2007;Cole et al., 2007). Although still much debated, recent studies have contributed the downward transport of soilrespired CO2 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 long 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 CO2 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 have deep roots (Gupta 535 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 deep subsurface. As the pace of climate change accelerates, summer droughts are expected to intensify, which can potentially channel more CO2 into the deeper subsurface via deepening roots.
With constraints from soil CO2 data, the simulated CO2 production rates in GrassPF and WoodPF are similarly at ~ 0.6 mol 540 C/m 2 /yr under the infiltration rate of 0.37 m/yr ( Figure 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 /yr) 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 (Figure 4). In other words, more soil-respired CO2 can transport to groundwater and become stored there for centuries to millennium before entering a stream. At the short time scale, this enhanced downward transport will reduce CO2 545 escape back into the atmosphere. This may explain the observations at Konza that woody encroachment increased NPP however soil CO2 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 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. 550

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 pCO2 drives variations of Ca and DIC concentrations in soil water and groundwater. We derived equations based on reaction thermodynamics (Eq. 2-3) to estimate Ca and DIC as a 555 function of pCO2 and temperature, which has 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 deep subsurface with abundant calcite. The results show that deeper penetration of roots and higher recharge flow enhanced CO2-carbonate contact. At an infiltration rate of 3.7 m/yr, calcite weathering flux in woodlands was 207% higher than that in grasslands. At 3.7×10 -2 m/yr, the weathering flux in woodland was 17% higher. The hydrological 560 impacts on carbonate weathering were more than 10 times higher than the biogeochemical impacts via elevated soil CO2 alone.
The modeling demonstrates that the weathering rates depend on flow partitioning and relative magnitude of water contact time in the deep carbonate zone. At the century scale, with the assumed 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 565 et al., 2017b), and therefore speed up landscape evolution (Phillips, 2005). It alludes to the importance of considering changes in https://doi.org/10.5194/bg-2020-180 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. subsurface hydrological flows associated with shifts in vegetation types and measuring rooting characteristics. This is particularly relevant as we assess the effects of climate change, land cover and elevated atmosphere CO2 concentrations on chemical weathering and carbon cycles.
Data availability. All data used for model parameterization can be acquired from http://lter.konza.ksu.edu/data. The input files 570 necessary to reproduce the results are available from the authors upon request (lili@engr.psu.edu). Competing interests. The authors report no conflicts of interest. 575