Abstract
It has long been recognized that the advent of vascular plants in the Paleozoic must have changed silicate weathering and fundamentally altered the long-term carbon cycle. Carbon cycle models are frequently employed to quantify the effect of this state change in the Earth system. These models have suggested that plants likely played a key role in modulating atmospheric CO2 in the past, with the largest plant-induced changes to the system in the late Paleozoic. These studies have, for the most part, focused on the effects of plants on weathering via their impacts on the soil environment. Yet, plants also modify the hydrological cycle, which may also have had implications for global weathering. Here, we evaluate the consequences of plant evolutionary innovation that have not been previously incorporated into carbon cycle models by coupling a one-dimensional vapor transport model to a reactive transport model of silicate weathering. Using this cascade of models, we investigate: 1) how evolutionary shifts in plant transpiration may have enhanced silicate weathering through increased downwind transport of water vapor to continental interiors; 2) the importance of deeply-rooted plants and their associated microbial communities in increasing soil CO2 and weathering zone length scales; and, 3) the coupled effect of these two processes on weathering rates.
The hydrologic balance for our modeling approach is framed by energy/supply constraints (encapsulated through Budyko relationships) calibrated for minimally vegetated-, vascular plant forested-, and angiosperm-worlds. Using constraints for atmospheric vapor-transport over continents and terrestrial weathering fluxes, we find that the emergence of widespread transpiration and associated inland vapor recycling associated with the advent of deep-rooted vascular plants over the later Devonian increases weathering solute concentrations by a factor of 1.64 to 2.39 for a fixed atmospheric CO2. The later domination of ecosystems by angiosperms in the late Cretaceous and Cenozoic, and the subsequent increase in transpiration fluxes, increased weathering solute concentrations by 7 to 55 percent leading to a cumulative total increase in weathering solute concentrations by a factor of 1.70 to 2.55. Partitioning of the vapor recycling effects and feedbacks on the weathering thermodynamics indicates that 55 percent of the increases in weathering concentrations are due to the thermodynamic effects of soil CO2 with the remaining 45 percent attributable to hydrologic effects. Our estimates of the relative changes in weathering solute concentrations caused by land plant evolution are of a similar magnitude to relative flux scaling relationships implemented in existing carbon cycle models such as GEOCARBSULF (Berner, 2006), COPSE (Bergman and others, 2004), and GEOCLIM (Le Hir and others, 2011). As Phanerozoic plant evolution resulted in a more efficient generation of weathering solutes, a weaker dependence of silicate weathering on physical forcing factors such as temperature, pH of rainwater, tectonics, and prevailing lithology, must have resulted. Consequently, we postulate that terrestrial plant evolution likely contributed to a more stable climate system over the Phanerozoic.
INTRODUCTION
Earth's climate has been remarkably stable over the last 4 billion year history. The relative stability of the Earth's climate has been attributed to a negative feedback on climate imposed by the weathering of silicate rocks (for example, Walker and others, 1981; Berner and others, 1983; Berner 1991, 2004, 2006; Berner and Caldeira, 1997; Maher and Chamberlain, 2014). Weathering reactions occur subaerially, and to a lesser degree in association with hydrothermal systems in the ocean (Coogan and Gillis, 2013; Mills and others, 2014). Weathering reactions convert carbon in the ocean-atmosphere system to carbonate alkalinity (bicarbonate, HCO3−, and carbonate, CO32−, ions), which then combine with cations (primarily Ca2+ and Mg2+) that are deposited as calcium and magnesium carbonate in marine sediments. The net result is the removal of carbon from the ocean-atmosphere to the solid Earth (Urey, 1952; Zeebe and Caldeira, 2008)—a process that is often referred to as the silicate weathering negative feedback where increased CO2 concentrations lead to more efficient removal of CO2 via silicate weathering. The overall weatherability of the Earth surface determines the efficiency with which silicate weathering converts CO2 to bicarbonate through reaction with subaerially exposed silicate rocks. This weatherability is determined by two primary factors: (1) the reactivity of rocks within the weathering zone and (2) the areal extent and distribution of continents subject to weathering processes. All else being equal, greater reactivity of rocks and larger continental area will result in enhanced weatherability, which results in a lower steady-state pCO2. Moreover, the stronger negative feedback, the coupling between silicate weathering fluxes and climate, should provide greater resistance to CO2 perturbations (Kump and Arthur, 1997; Caves and others, 2016a).
Yet, physical determinants of weatherability such as continental emergence, tectonics, and variations in the lithologic character of exposed rocks, are cyclical processes over geologic time lacking unidirectional state changes (Kump and Arthur, 1997; Gibbs and others, 1999; Kump and others, 2000; Berner, 2006, 2008; Mills and others, 2014; Ibarra and others, 2016; Lee and others, 2018). In contrast, there is unambiguous evidence for biological influences on the weathering zone, including the physical and chemical influence of land plants, which have dramatically and unidirectionally changed through geologic time (Berner, 1987; Schwartzman and Volk, 1989; Volk, 1989; Berner, 1992; Algeo and others, 1995; Algeo and Scheckler, 1998; Knoll, 2003; Pagani and others, 2009; Le Hir and others, 2011; Taylor and others, 2012; Boyce and Lee, 2017).
Land plants are thought to alter weatherability through two key mechanisms (Berner, 1987, 1992; Drever, 1994). First, plants and associated soil microbes enhance weathering directly (Barker and others, 1997). This biota respires CO2, resulting in pCO2 levels in soils that are often several orders of magnitude greater than atmospheric pCO2 (Brook and others, 1983; Kelly and others, 1998; Montañez, 2013) and providing a critical source of protons that contribute to the weathering of rocks (Moulton and Berner, 1998; Moulton and others, 2000; Andrews and Schlesinger, 2001; Berner and others, 2003). Breakdown of rocks is also enhanced by attack of mineral grains by plant fungal symbionts (Blum and others, 2002; Quirk and others, 2014), which increases the mineral surface area accessible for chemical weathering processes. Second, plants may enhance weathering indirectly by altering global hydrology through transpiration due to the evolution of stomata, vasculature, and deep roots that can access groundwater that would otherwise be shielded from evaporation (Shukla and Mintz, 1982; Lee and Boyce, 2010; Boyce and Lee, 2010, 2017). Transpiration increases precipitation globally (Shukla and Mintz, 1982; Kleidon and others, 2000; Lee and Boyce, 2010; Boyce and Lee, 2010, 2017; Alkama and others 2012), permitting both the transport of water vapor into continental interiors and extensive local recycling of water vapor. Because silicate weathering fluxes are fundamentally mediated by the hydrologic cycle (Bluth and Kump, 1994; White and Blum, 1995; Dessert and others, 2003; Donnadieu and others, 2004; Labat and others, 2004; Le Hir and others, 2011; Maher, 2011; Maher and Chamberlain, 2014; Ibarra and others, 2016), we hypothesize that increases in transpiration may have increased the continental area that is accessed for silicate weathering. In short, not only do plants alter the local weathering environment, but they may also expand the area over which silicate weathering can occur through downstream hydrologic feedbacks and teleconnections (for example, Boyce and Lee, 2010, 2011; Green and others, 2017).
Previous studies have attempted to quantify this role of plants on the silicate weathering feedback in a number of ways. The general approach of early studies was to use empirically (for example, Knoll and James, 1987; Moulton and Berner, 1998) or experimentally (for example, Blum and Lasaga, 1988) derived weathering rates to adjust flux scaling terms in 1-D carbon cycle models for the colonization of continents by vascular plant trees in the Devonian and Carboniferous (Schwartzman and Volk, 1989; Berner 1991, 1992, 2006) and the evolution and expansion of angiosperms in the Cretaceous and Early Tertiary (Volk, 1989). More recent work has utilized global climate models configured with paleo-boundary conditions to examine holistically the effect of land plant evolution and colonization on planetary albedo and hydrology (for example Le Hir and others, 2011). Our approach is distinct from these studies as we explicitly examine the effect of plants on silicate weathering through changing hydrology and compare this to the effect due to changing reaction rates. As first hypothesized by Berner (1992) and Drever (1994), based on early climate modeling results by Shukla and Mintz (1982), colonization of continents by plants must have increased rainfall in the interior of continents through the downstream recycling of water vapor by transpiration and thus have increased global weathering fluxes. For example, plants in a coastal region will locally alter the weathering environment through increased weathering rates. Additionally, their transpiration will increase precipitation in continental interiors. This downwind recycling of water will both supply water to the weathering zone and increase plant productivity, thus increasing soil pCO2 across wide swaths of continents. Consequently, the presence of land plants may increase the quantity of water in continental interiors available for weathering and increase primary productivity, further contributing to an increase in global weatherability.
To test this hypothesis, we examine the relative contributions of these two mechanisms to silicate weathering using a one-dimensional vapor transport model (Hendricks and others, 2000; Mix and others, 2013; Chamberlain and others, 2014; Winnick and others, 2014) that is linked to a reactive transport model of silicate weathering (Maher, 2011; Maher and Chamberlain, 2014; Ibarra and others, 2016, 2017; Wymore and others, 2017).
MODELING FRAMEWORK
In this paper, we are explicitly interested in modeling steady-state solute concentrations from silicate weathering on a continent before and after two major plant evolutionary innovations. Following Berner (1992), we assume a baseline condition (referred to as a minimally vegetated world) of precursor microbial communities that likely colonized continents before the advent of vascular plants in the Silurian (DesMarais, 1990). Specifically, this refers to a land surface with no transpiration or biotic amplification of soil CO2, but with a vegetated albedo similar to modern. First, we test how weathering changed after the proliferation of deeply rooted (greater than a few centimeters) vascular land plants (particularly gymnosperm seed plants and their likely progymnosperm ancestors) that occurred in the later Devonian and Carboniferous. These non-angiosperm vascular land plants would have allowed for increased total evapotranspiration rates and vapor recycling through the evolution of deep roots to access groundwater and vasculature to transport that water to the stomata to support photosynthetic gas exchange. Second, we examine how the spread and dominance of angiosperms in the Late Cretaceous and Cenozoic (Crane and Lidgard, 1989; Burnham and Johnson, 2004; Jaramillo and others, 2006) affected weathering. We do not attempt to model an un-vegetated ‘desert world’ scenario (for example, Fraedrich and others, 1999; Kleidon and others, 2000). The large changes in albedo involved in such a scenario generate large changes in temperature and potential evapotranspiration that would obscure our assessment of the importance of plant transpiration (for example, Kleidon and others, 2000; Le Hir and others, 2011; Boyce and Lee, 2017). Thus, we do not consider changes in albedo in our modeling and do not simulate ‘desert world’ or ‘bare world’ scenarios.
We focus on how plant evolution affects the concentration of weathering solutes, rather than addressing how plant evolution changes weathering fluxes, to isolate the effect of plants on weatherability (independent of changes in runoff). The weathering flux is related to the concentration of weathering solutes:
where Fsil is the weathering flux [moles/t], q is runoff [L/t], and C is the concentration of weathered solutes [moles/L]. On geologic timescales (> 1 Ma), Fsil must equal the volcanic and metamorphic flux of carbon (Berner and Caldeira, 1997; Broecker and Sanyal, 1998; Zeebe and Caldeira, 2008; Caves and others, 2016a). Runoff, in turn, is controlled by climate, including both the strength of the greenhouse effect and solar insolation, though we note that the colonization of continents by plants will change the radiative balance of the Earth through albedo changes (Fraedrich and others, 1999; Kleidon and others, 2000; Alkama and others, 2012), a factor that we do not consider in this study.
For each of the scenarios, we calculate the change in the concentration of weathered solutes through two linked models (fig. 1): The first model (Model 1) allows (evapo-) transpiration to be treated as a variable with integrated precipitation and vapor recycling parameters as outputs. Specifically, Model 1 consists of an advection-diffusion-reaction atmospheric vapor transport model (Hendricks and others, 2000; Winnick and others, 2014; Chamberlain and others, 2014; see Appendix 1 for derivation and details) with a mass-balance scheme to partition surface water derived from precipitation into runoff and evapotranspiration. The second model (Model 2) resolves the weathering concentrations based on changes to plant/bedrock interaction. Model 2 is based on a recently developed reactive transport-weathering model (see Maher, 2011; Maher and Chamberlain, 2014; Ibarra and others, 2016; Wymore and others, 2017 for details). The coupling of Model 1 and 2 is through a runoff term (q, fig. 1), common to both models (output of Model 1, input into Model 2). Model variables are described in table 1 with assumed values given in table 2. (The Code is available at http://earth.geology.yale.edu/%7eajs/SupplementaryData/2019/IbarraCode.)
Schematic of physical set-up of modeling approach for the primary set of boundary conditions (constant flux to constant gradient) considered for Model 1 and the framework for weathering used in Model 2 (see text for details). ET in the lower panels refers to evapotranspiration (ET = E+T, evaporation plus transpiration). Note that throughout the text we denote evaporation (with or without plants) as ET for simplicity.
Definition and units of model variables
Variable constraints for input parameters
Model 1: Vapor Transport Model
To capture the effect of plants on vapor recycling and downwind precipitation, evapotranspiration, and runoff, we use an advection-diffusion-reaction scheme coupled to a mass-balance formulation of evapotranspiration (ET) (sensu Budyko) to ensure mass conservation of water (fig. 1, left panels). Following Hendricks and others (2000), we track water vapor in the atmosphere with the following advection-diffusion-reaction equation:
We solve this in steady state, such that dw/dt = 0. Here, ∇ is the gradient operator, w [kg/m2] is the integrated column water vapor, D [m2/s] is the diffusion coefficient for atmospheric moisture, representing transport by atmospheric eddies, u [m/s] is the velocity, representing the average motion of the atmosphere, P [kg/m2/s] is precipitation, and ET [kg/m2/s] is evaporation plus transpiration. We simplify equation 2 by assuming transport along a one-dimensional storm track (left to right in fig. 1, left panel). We further simplify equation 2 by assuming that precipitation (P) is a function of the available water vapor (w), such that:
where kp [1/sec] is a precipitation efficiency term scaled to the global average residence time of modern atmospheric vapor (∼ 4.5 days) (Läderach and Sodemann, 2016).
To partition precipitation into ET and q and conserve water mass, we use the Budyko equation (Budyko, 1974; Choudhury, 1999; Zhang and others, 2004; Koster and others, 2006; Gerrits and others, 2009; Greve and others, 2015). This framework allows us to place quantitative bounds on ET and P and is based on water and energy balance equations (Budyko, 1974; Arora, 2002). The form of the one-parameter Budyko equation we use (Choudhury, 1999; Zhang and others, 2004; Greve and others, 2015) is:
Here, Ep [kg/m2/s] is potential evapotranspiration and ω is a dimensionless variable ranging from 1 to >5 that characterizes the efficiency with which P is converted to ET. Increasing ω represents increasing efficiency of converting P to ET (Zhang and others, 2004; Koster and others, 2006; Gerrits and others, 2009; Greves and others, 2015). Thus, a geologically implausible end-member with no evapotranspiration would be represented by ω = 1 and ET = 0, where all water is partitioned to runoff. In contrast, the estimate of ω for the modern world is 2.6 and represents the modern average end member (Greve and others, 2015). To mimic the evolution of plants, we increase ω to capture the effect of increasing transpirative capacity from a minimally-vegetated world (fig. 1, upper panels) to a vascular plant-world to a flora dominated by angiosperms (both fig. 1, lower panels).
Boundary Conditions for Model 1
We apply our hydrologic model over a domain of 10,000 km in length, as this is similar to the zonal length of the large continental landmasses in the late Paleozoic (Domeier and Torsvik, 2014). Our parameterization of the diffusional constant (D) encapsulates a number of complicated processes that are not explicitly resolved in our one-dimensional advection-diffusion-reaction equation and is, further, dependent upon the precise length-scale of interest. Here we calculate D based on modern water isotope observations (Winnick and others, 2014) which suggest that, over mid-latitude areas of continents, both advection and diffusion are roughly equally important in equation 1, corresponding to a Peclet number (Pe = u/D) of 1. Note that the use of ‘diffusion’ here is more traditionally referred to as ‘dispersion’ in hydrologic reactive transport modeling, although the use of ‘diffusion’ to describe this process is consistent with literature on transport of water vapor by atmospheric eddies (for example, Hendricks and others, 2000; Winnick and others, 2014).
To solve the advection-diffusion-reaction equation, we evaluated a series of boundary conditions (see Appendix 1) and present the most realistic set in the main text. We set the boundary condition of a constant water flux (constant flux boundary condition: win = constant) entering the continent on the western (left) coastal boundary, and for the second boundary condition, we impose no gradient in water vapor exiting the continent on the eastern (right) coastal boundary (no gradient boundary condition: dw/dx = 0). This simulates a situation where the water vapor content has reached equilibrium with respect to transpiration and precipitation and no further net reduction occurs. Some of the precipitated water is returned to the atmosphere via ET and is advected off the continent; consequently, total runoff decreases with increased transpiration as more water vapor is ‘lost’ at the eastern edge of the continent by advection. Such a decrease in globally-averaged runoff with increasing transpiration is also observed in climate models that evaluate the hydrologic impact of plant evolution (Kleidon and others, 2000; Alkama and others, 2012). In the Appendix, we also show simulations that do not cause total runoff to decrease by imposing constant flux boundary conditions on both edges of the domain, where the western (left) coastal boundary is a positive constant water flux and the eastern (right) boundary is a negative constant water flux leaving the domain.
Because P is a function of w, ET is therefore also a function of w; consequently, there is no analytical solution to the combined set of equations above (equations 2–4). Therefore, we use a numerical discretization scheme to solve the full advection-diffusion-reaction equation (Appendix 1). Values for the above equations are given in table 2 and are derived from the mean terrestrial values from well-constrained climate models and observations (Boyce and Lee 2010, 2017; Greve and others, 2015). We calculate runoff (q) as the difference between P and ET at each grid cell to calculate weathering rates per grid cell (Model 2).
Model 2: Chemical Weathering in Soils
To calculate weathering solute concentrations (C [μmol/L]), we use the solute production equation derived in Maher (2011) and Maher and Chamberlain (2014) (see also Ibarra and others, 2016, 2017; Wymore and others, 2017):
where Cmax [μmol/L] is the maximum solute concentration for a given aqueous species (for example, SiO2, HCO3−) dictated by the “thermodynamic limit” of the reaction between primary and secondary minerals, q is runoff [m/yr] from Model 1, Dw [m/yr] is the Damköhler coefficient, and τ is an exponential scaling factor derived from an assumption of the water transit time distribution (τ = e2) (see derivation in Maher and Chamberlain, 2014). The Damköhler coefficient is a measure of the reactivity of the weathering zone and is derived by factoring out runoff from the Damköhler number, which is the ratio of the fluid travel time to the time required to reach equilibrium (see further explanation in Johnson and DePaolo, 1994; Maher and Chamberlain, 2014; Li and others, 2014; Ibarra and others, 2016; Wymore and others, 2017). The Damköhler coefficient [m/yr] is given by:
where L is the reactive flowpath length [m], φ is the porosity, and Teq is the equilibrium time [years] of the net weathering reaction, which can be related to Cmax as Teq = Cmax/Rn, where Rn[μmol/L/years] is the net dissolution rate of the weathering reaction of interest (Maher, 2010, 2011; Maher and Chamberlain, 2014; Ibarra and others, 2016; Wymore and others, 2017). The reactive flowpath length and porosity describe the length of the weathering zone (soil depth in fig. 1).
Using runoff (q) from Model 1, the weathering solute concentration C (eq. 5) is then calculated for each grid cell over our continental landmass. The methods by which Cmax is calculated are described in the section “Effect c” below. Within this framework, we aim to model the consequences of plants through (a) an expansion of water transport across a continent, (b) greater local recycling of water vapor, and (c) a higher [thermodynamic potential] for chemical weathering. Effects (a) and (b) capture the effects of plants on hydrology, whereas effect (c) focuses on the effects of plants on weathering reactions. Below, we describe how we alter model variables to simulate silicate weathering for our three plant worlds (minimally vegetated world, deeply-rooted vascular plant world, and angiosperm world).
Effect a: Continent-scale Changes in Hydrology (ω) Through Plant Transpiration
In Model 1, we vary ω to capture the effects of increased areal extent of precipitation as transpired water advects inland over 1) a low evaporation, minimally vegetated land surface, 2) a vascular plant world with intermediate evapotranspiration, and 3) an angiosperm world with maximum transpiration capacity (Boyce and Lee, 2010, 2017; fig. 1).
To constrain ω for minimally vegetated-, vascular plant-, and angiosperm-worlds, we use output from climate model simulations that examined the effect of vascular plants and angiosperms on the global hydrological cycle (Boyce and Lee, 2010, 2011; Lee and others, 2012; Boyce and Lee, 2017). As shown in figure 2A, there appears to be a functional limit on ET for vascular plants compared to the angiosperm world. We fit equation 4 by solving for the dimensionless quantity (ω) that dictates the form of the equation for different hydrologic environments (figs. 2A, 2B). The best-fit curves (shown in fig. 2A) show that ω increases from minimally vegetated- to vascular plant- to angiosperm-worlds within generations of climate model experiments, although the specific values of ω differ greatly (fig. 2A; table 2). In addition to using climate model simulations, we also use a compilation of direct, hydrological measurements (Bosch and Hewlett, 1982; Brown and others, 2005) at paired watersheds that compare vegetated vs. non-vegetated (that is, clear-cut) treatments (irrespective of plant composition). These modern vegetated (both angiosperm and gymnosperm-dominated) watersheds have higher ET relative to their paired, clear-cut counterparts (fig. 2B). To estimate ω from these paired watersheds, we combine estimates of Ep (Trabucco and Zomer, 2009) with the observed change in runoff. The vegetated catchments have a ω value of 2.6, which is in concordance with estimates from larger datasets of modern catchments (Greve and others, 2015), whereas clear-cut, non-vegetated sites are characterized by a lower ω value of 1.9 (fig. 2B). Based on these constraints from our Budyko analysis, we select a range of values to represent end-member scenarios for minimally vegetated (ω = 1.5; fig. 1 upper panels), and angiosperm (ω = 2.5; fig. 1 lower panels) worlds, with an intermediate value representing the vascular plant world (ω = 2.2). As a sensitivity test, we also present a broader range of ω values in our analysis from no evapotranspiration (ω = 1.0) to highly optimized inland vapor recycling (ω > 3.5). As such, our modeling is intended to capture the direction and trends in weathering, despite uncertainties in constraining ω.
Determination of the effect of plants on water fluxes in Budyko space. (A) Climate model simulations from Boyce and Lee (2010) (see also Lee and Boyce, 2010; Lee and others, 2012) using NCAR CAM3.0/CLM 3.5. Each land grid cell is represented by a gray dot. Net irradiance (Rnet) is the sum of incoming and outgoing longwave and shortwave radiation and can be equated with potential evapotranspiration, Ep, (for explanation see Roderick and others, 2014), which was then converted into water units using the latent heat of vaporization for water. Black dots are the binned average ET/P for every increment of 0.25 Rnet/P, and the black line is the best-fit line to the binned averages (see reasoning in Koster and others, 2006). Note that in the no transpiration case there are less land points shown because they either violate mass balance (E>P) and were not included, or sit to the right (Rnet/P > 3) of the axes shown. (B) Compilation of field experiments of deforestation or clear-cutting where the necessary climatological parameters are available to assess the impact of decreased transpiration (see compilation in table A2). The modern observation of ω (2.643) agrees well with the median value from Greve and others (2015), which analyzed more than 400 modern catchments. Budyko curves were fit using the ‘nls2’ R package (Grothendieck, 2013).
Effect b: Local Water Vapor Recycling (f) Through the Depth of Plant Rooting
Two key, coupled changes for hydrologic cycling are induced by the evolution of deep-rooted plants. First, because plants increase ω, both ET and P increase, resulting in greater recycling of precipitation (Kleidon and others, 2000; Boyce and Lee, 2010, 2011; Alkama and others, 2012). Second, this recycled precipitation can more frequently transit the entirety of the weathering zone without entering runoff due to the presence of deep-rooted vascular plants that can draw upon deep water for transpiration. Consequently, a fundamentally new weathering-zone dynamic appears after the colonization of land by plants whereby water availability in the weathering zone may no longer be strictly delineated by runoff, but also by the frequency with which water is recycled through the weathering zone (Drever, 1992).
To account for deep-rooted vascular plants, we modify equation 5 in Model 2 and add a recycling factor (f):
Here, the recycling factor (f) is the ratio of precipitation to runoff (P/q). This form of the equation applies to a situation in which plant roots extend to or beyond the actively weathering zone (that is, near the soil-saprolite-bedrock interface)—such as large trees do today (Stone and Kalisz, 1991; Ehleringer and Dawson, 1992; Jongmans and others, 1997; Landeweert and others, 2001; Lucas, 2001). In this scenario, the transpired water moves through the weathering zone but does not enter runoff and instead would be recharged to the atmosphere (green ET arrow in lower right panel of fig. 1), effectively increasing the residence time of each water parcel in each grid cell as vapor is transported into the continental interior. Thus, equations 5 (sans recycling factor) and 7 represent two end-members that mimic a shallow rooting versus a deep rooting environment.
Since f is calculated for every grid cell of the domain, and thus evolves spatially along the continental transect, this formulation satisfies the observational and modeling evidence (primarily from precipitation stable isotope measurements and analyses of wind and humidity gradients) that downstream recycling is important, particularly in inland continental settings such as Eurasia and southern Africa today (Salati and others, 1979; Brubaker and others, 1993; van der Ent and others, 2010; Caves and others, 2015; Bershaw, 2018). For example, an analysis of atmospheric convergence by Brubaker and others (1993) found that the local contribution from recycled evapotranspiration to monthly precipitation in many regions accounts for up to 40 percent of precipitation; in Central Asia, this percentage may reach beyond 80 percent (van der Ent and others, 2010). Additionally, both atmospheric vapor content and hence continental vapor recycling increase in the simulations (Boyce and Lee, 2017; Lee and Boyce, 2010; Boyce and Lee, 2010, 2011) that are used to parameterize Model 1. Our specific formulation for incorporating the appearance of deep-rooted plants may not precisely capture the complex dynamics that occur between precipitation, rooting-depth, and soil moisture; nevertheless, we view equation 7 as a first-order attempt to capture a fundamentally novel switch between a shallow-rooted, low recycling world and a deep-rooted, high recycling world.
Effect c: Enhancement of Chemical Weathering Reactions Through Productivity
In addition to generating a more active terrestrial hydrological cycle, plants also increase the thermodynamic limits of chemical weathering reactions. This effect has been shown by previous studies that emphasize the effects of plants on precipitation, soil moisture, and primary productivity, which affect soil pCO2 and thus chemical weathering reactions (for example, Berner, 1992; Drever, 1994; Kelly and others, 1998). In our model framework, these effects are encapsulated in the Cmax term (Model 2; eqs. 5 and 7)—the thermodynamic maximum weathering solute concentration as determined by soil conditions. We carried out equilibrium reaction path simulations to calculate Cmax from weathering of granitic bedrock to secondary aluminosilicate minerals as a function of soil pCO2 with the “React” module of The Geochemist's Workbench (GWB; Bethke, 2007). The following standard states were adopted for the thermodynamic calculations: unit activity for pure minerals and H2O at any temperature and pressure; for aqueous species, unit activity of the species in a hypothetical 1 mol/kg solution referenced to infinite dilution at any temperature and pressure; for gases, unit fugacity of a hypothetical ideal gas at 1 bar of pressure. Soil pCO2 and fugacity are related by the fugacity coefficient, λ, by fCO2 = λ × pCO2. Similarly, Cmax is obtained from the aqueous species concentration (μmol/l) of HCO3− as the weathering solute representative of silicate weathering, which is related to activity through the activity coefficient, γ, by a = γ × C. Thermodynamic data is from the Thermoddem database (version llvl2; Blanc and others, 2012).
Using a fixed fugacity of CO2 representing an open system with an infinite reservoir of CO2, pH was determined by charge balance and allowed to evolve over the course of the reaction. An initial mineral phase of Na-rich plagioclase (An20Al80, 20 mol% anorthite and 80 mol% albite) was added to water for a fluid:rock ratio of 100. The secondary mineral phase allowed to precipitate was the hydrated aluminosilicate halloysite (Al2Si2O5(OH)4). This reaction is
We perform this calculation across a wide range of soil pCO2 (200–10,000 ppm). The CO2 (aq) is set to be in equilibrium with the specified pCO2 at the given temperature. To test the sensitivity of our calculations to different primary and secondary mineral phases and to temperature, we vary the composition of the initial feldspar in equation 8, the secondary mineral phase, and the reaction temperature (5–25 °C). The initial water composition is from Maher and others (2009) (see their Table 2) and is representative of dilute rainwater. Note that the simulations shown here using GWB are analogous to the ‘open-system’ soil column reactive transport simulations presented by Winnick and Maher (2018) (see their Figure 3).
Results of GWB titration experiments (Bethke, 2007) reacting rainwater (see Maher and others, 2009) with feldspars to form secondary minerals to determine Cmax values of HCO3−. Batch calculations were carried out to completion at a fixed fugacity of pCO2 (that is, an infinite reservoir) to explore the effect of temperature (left panel), secondary mineral end-member, and primary mineral starting material.
Several results follow from this analysis and confirm previous modeling (Maher, 2011; Maher and Chamberlain, 2014; Ibarra and others, 2016; Thomas and others, 2016; Winnick and Maher, 2018). First, these calculations show that in all cases, HCO3− and SiO2(aq) Cmax values will increase with increasing soil CO2 regardless of parent material, reaction temperature, or specific clay forming reaction (fig. 3). Second, the Cmax values are relatively insensitive to temperature over reasonable Earth surface values (fig. 3A) (as previously shown by Maher and Chamberlain, 2014; Winnick and Maher, 2018). Third, the Cmax values are relatively sensitive to the specific soil forming reaction and parent material (figs. 3B and 3C). For example, we use the reaction that forms halloysite (following Maher, 2011; Maher and Chamberlain, 2014) because it results in HCO3− and SiO2(aq) Cmax values in close agreement (factor of 2 to 10 scaling) with modern solute data (Drever and Zobrist, 1992; Arthur and Fahey, 1993; Berner and Rao, 1996; Bormann and others, 1998; Moulton and Berner, 1998; Gaillardet and others, 1999; Moulton and others, 2000; Maher and Chamberlain, 2014; Ibarra and others, 2016, 2017; Wymore and others, 2017). This reaction estimates the lowest values of Cmax for a given soil CO2 (fig. 3), whereas using kaolinite as the secondary mineral predicts much higher values of Cmax. Using Na-rich plagioclase feldspar An20/Al80 gives intermediate Cmax values between that of weathering of potassium feldspar (lowest Cmax) versus pure anorthite (highest Cmax). Critically, our sensitivity analysis shows that the slopes of the relationship between soil CO2 and Cmax for different primary and secondary minerals are very similar (with the exception of 100% anorthite as the primary mineral). Thus, although the absolute Cmax value for each scenario may depend on the assumed generalized reaction stoichiometry, the change in Cmax—and therefore the effect of plants on weathering—as soil CO2 increases is not sensitive to the reaction chosen.
In addition to using GWB simulations, we also compare to data from Moulton and Berner (1998) (as well as Moulton and others, 2000; Oskarsson and others, 2012) that compared weathering fluxes from a barren (minimally vegetated) and a vegetated basaltic catchment in Iceland. We use the distribution of concentration measurements and calculate the 95 percent maximum (percentile) of the data (fig. 4) to estimate Cmax values (following Dixon and von Blanckenburg, 2012; Maher and Chamberlain, 2014; Ibarra and others, 2016). We find that the Cmax values increase ∼1.8-fold between barren and vegetated catchments, in agreement with the range suggested by our GWB simulations, although Iceland soil pCO2 for the catchments studied by Moulton and Berner (1998) are unknown and represent relatively low Cmax concentrations relative to a synthesis of modern basaltic catchments (Ibarra and others, 2016).
Summary of stream concentration data by Moulton and Berner (1998) for barren (minimally vegetated) and vegetated modern landscapes. Following previous work (Dixon and von Blanckenburg, 2012; Maher and Chamberlain, 2014; Ibarra and others, 2016), Cmax is calculated assuming the 95% maximum (horizontal black lines) calculated from the mean concentrations (white dots) and observed distributions (Table 1 of Moulton and Berner, 1998). Average values are taken for Cmax values for plant and barren landscapes (dashed lines).
Because Cmax is sensitive to soil pCO2 and soil pCO2 is sensitive to climate (Cotton and Sheldon, 2012; Oskarsson and others, 2012; Montañez, 2013), we vary Cmax across our continental transect as a function of precipitation. That is, Cmax is allowed to evolve (decrease) laterally across the domain as precipitation rates decrease due to the loss of water as runoff. To provide Cmax estimates for each grid cell on our continental transect, we first calculate soil pCO2 as a function of precipitation using data compiled by Cotton and Sheldon (2012) and Cotton and others (2013), which show that soil pCO2 is positively correlated with precipitation for the mid-latitude sites they examined, primarily in western North America (see fig. A4). We note that, particularly in semi-arid midlatitude sites, there are large seasonal variations in soil CO2 (Breecker and others, 2009; Oerter and Amundson, 2016), however, given the steady-state nature of our calculations, we are unable to account for seasonal fluctuations. Further, weathering likely occurs during the period when precipitation occurs, which is correlated with increasing productivity and hence higher soil pCO2. Using the regression shown in figure A4, we calculate how Cmax of HCO3− varies with soil pCO2, assuming that Cmax varies between approximately 250 μmol/L in low soil pCO2 (200 ppm) locations to approximately 1200 μmol/L in high soil pCO2 environments using a regression between the reaction path Cmax of HCO3− results and the simulated soil pCO2 values (fig. 3, solid black line in each panel; fig. A4). As discussed above, this relationship is in qualitative agreement (in direction and magnitude) with the Moulton and Berner (1998) observations.
Modeling Setup and Assumptions
Our modeling is a simplified representation of how plants impact chemical weathering and the hydrologic cycle. As a result, we have made several key assumptions that may limit the potential for this work to be scalable to global weathering relationships. The limitations of this modeling setup include:
Vapor transport and land surface-atmosphere exchange is constrained to follow the Budyko curve formulation. While representative datasets spanning different climates (for example Koster and others, 2006; Greve and others, 2015) and climate model experiments (for example Roderick and others, 2014; see also fig. 2) follow Budyko landscape parameters (ω values) centered on ∼2.6 today, continental transects as simulated in our model do not necessarily follow this trend. In fact, some evidence from the mid-latitudes suggests lower ω values in mountainous continental interiors (see Greve and others, 2015).
Ocean evaporation and thus the flux into the domain is held constant as prescribed by the constant flux boundary condition. This flux may change in response to the biophysical effects of plant evolution, including changes in albedo and surface roughness, which we do not consider in this study. Nevertheless, this boundary condition and consequent model results are supported by a climate model (whose results are utilized here; from Boyce and Lee, 2010): despite predicting a slight decrease in global runoff, this model also predicts significant changes in precipitation patterns with changing plant transpiration.
The two models are calibrated to, and thus representative of, a mid-latitude continent-dominated world for two reasons. First, we have chosen a Peclet number of 1 (see above) for the vapor transport module, which influences the degree of rainout and advection of moisture inland. Second, in modern mid-latitude settings, many surfaces are young and not ‘supply limited’. This is important because flat tropical environments today produce some of the lowest area-normalized weathering rates despite high precipitation and runoff due to thick regolith (weathering zone) development in humid, tropical regions (for example, West and others, 2005; Godderis and others, 2006; Hewawasam and others, 2013).
Future work that extends our results may include coupling global climate model output (for example, Boyce and Lee, 2017) with weathering fluxes (for example von Blanckenburg and others, 2015), or directly modeling the influence of plants on transpiration and weathering reactions in intermediate complexity models (such as cGENIE; Colbourn and others, 2015).
MODEL RESULTS
We find that the impacts of plants agree well with previously estimated effects of plant evolution on weathering; however, our model framework allows us to calculate how this cumulative effect is partitioned between hydrologic vapor recycling effects and thermodynamic effects. Our results support the hypothesis that plants will increase the areal distribution over which precipitation – and hence chemical weathering – occurs as vapor is recycled downwind. Due to only this hydrological change, we estimate that the colonization of vascular plants over a minimally vegetated world will increase solute concentrations by approximately a factor of 1.6 for a fixed atmospheric CO2 concentration. Further, we attempt to quantify the increase in chemical weathering due to plant/rock/groundwater interactions in the weathering zone. Our model demonstrates that this feedback plays a significant role in increasing weathering solute concentrations, potentially by a factor of about 1.5. Finally, we calculate that the evolution of angiosperms, depending on the efficiency of hydrologic recycling effects and thermodynamic effects by gymnosperms, increased weathering concentrations by an additional 7 to 55 percent.
Continental Transects of Hydrology and Weathering
Our calculations of the role of vapor transport by plants (Model 1) show that plants significantly increase the areal extent over which precipitation occurs over continents. In figure 5, we show simulations using ω values representative of no transpiration (minimally vegetated world; ω = 1.5), non-angiosperm vascular plants (ω = 2.2), and modern angiosperms (ω = 2.5) (table 2). As shown in figure 5, water vapor content in the atmosphere, precipitation, and evapotranspiration increase at each transition from a minimally vegetated landscape to a vascular plant world and to an angiosperm world. The effect of water replenishment of the atmosphere by evapotranspiration (fig. 5C) leads to a flatter distribution of runoff along the moisture transport pathway in both the vascular plant and angiosperm worlds (fig. 5D). Total runoff decreases relative to a minimally vegetated world by 17 percent and 23 percent, respectively. Additionally, at any point along the continental transect, the land surface has a greater recycling ratio in plant covered landscapes, with a greater effect for angiosperms than for gymnosperms (fig. 5E). Together, these results show that the hydrologic cycles in minimally vegetated-, in vascular plant-, and in angiosperm-worlds are distinct (fig. 5).
Water fluxes simulated by Model 1 using constant flux (left-hand boundary) to constant gradient (right-hand boundary) boundary conditions (see fig. 1 and text for explanation). Results show ω values representing the simulations from Boyce and Lee (2010) and modern observations as in figure 2. (A) Total water vapor content, (B) precipitation, (C) evapotranspiration, (D) runoff, (E) recycling ratio (P/runoff) and (F) the evolution of the vapor transport profiles in Budyko space.
In figure 6, we show how weathering-related parameters change across the continental transect as a function of hydrologic changes alone (panel A; effects a and b) or as a function of both hydrologic changes and thermodynamic feedbacks (panels B, C and D; effects a, b and c). Dw values, a measure of the reactivity of the weathering zone, are held at 0.03 m/yr in figure 6 (global average determined by Maher and Chamberlain, 2014). In a minimally vegetated world, soil pCO2 remains low (1700 ppm) relative to highly productive soils (10,000 ppm; Cotton and Sheldon, 2012; Montañez, 2013; Cotton and others, 2013) and is unchanged by precipitation because the deep soil expression of plant productivity is only weakly (or not at all) responsive to changes in precipitation. Soil pCO2 for the minimally vegetated world is estimated to be 1700 ppm from reconstructions of Cambrian and early Ordovician atmospheric pCO2 levels (Yapp and Poths, 1992; Mora and others, 1996; Foster and others, 2017), though we note that there is evidence for soil pCO2 elevated above atmospheric levels in pre-vascular plant soils (Yapp and Poths, 1994). In contrast, for all plant scenarios, soil pCO2 scales with precipitation (following the above described modeled relationships), with higher soil pCO2 near the coast and lower inland (fig. 6C). We also show a minimally-vegetated world with soil pCO2 scaling with precipitation (red dashed lines in fig. 6).
Solute fluxes simulated by Model 2 using inputs from Model 1 (fig. 5). (A) Solute [HCO3−] with hydrologic feedbacks (effects a and b) only. Lines correspond to A, A', B' and C' in figure 7 (middle and left panels). Soil pCO2 and Cmax [HCO3−] are invariant when only hydrologic feedbacks are included. (B) Solute [HCO3−]; (C) soil pCO2; (D) Cmax [HCO3−] including both hydrologic and thermodynamic feedbacks. Lines correspond to A, A”, B”, and C” in figure 7 (right and left panels).
Bicarbonate Cmax values, which are scaled to precipitation, follow soil pCO2 (fig. 6D) and decrease inland, while [HCO3−] increases inland (fig. 6B). Increasing [HCO3−] inland results from decreasing q from the coast to the interior. Consequently, the concentration of weathering solutes approach Cmax values at low runoff because of less dilution (Maher and Chamberlain, 2014; Ibarra and others, 2016; Wymore and others, 2017). This increase is also smaller when only hydrologic feedbacks (effects a and b) are included (fig. 6A), accounting for 43 percent of the increase in both the gymnosperm and angiosperm cases (table 3). If total runoff across the domain is held constant, using constant flux to constant flux boundary conditions, the trends and magnitude of concentration change are comparable (see fig. A5, panel F).
Summary of model results and GEOCARBSULF parameters based on Berner (2006) and Royer and others (2014)
Sensitivity to Changes in ω and Dw
The results presented in figures 6 and 7 show the cumulative impacts of global hydrologic change and local thermodynamic effects on solute concentrations.
Contour plots of weathering concentration (flux weighted) for varying values of Dw and ω. All weathering concentrations are normalized to the minimum vegetation scenario (point A on left panel) with no plant recycling feedbacks (Dw = 0.03 m/yr; ω = 1.5, Cmax = 1700 ppm; see text for details). Points correspond to curves in figure 6: A – minimally vegetated, A' – minimally vegetated with recycling, B' – vascular plants with recycling, C' – angiosperms with recycling, A” – minimally vegetated with recycling and thermodynamic effects, B” vascular plants with recycling and thermodynamic effects, and C” angiosperms with recycling and thermodynamic effects. Black arrow trajectories show increases in both local effects and global effects of increased downstream vapor recycling.
The effect of increased distribution of precipitation by plants (that is, plant recycling feedbacks) on weathering solute concentrations is shown in figure 7 (represented by ω values on y-axis). Here, for any given bedrock weathering zone (represented by varying Dw), [HCO3−] increases from minimally vegetated to angiosperm worlds. The absolute value of this effect is strongly dependent upon the erosional regime, which is encapsulated by the Dw value (x-axis of fig. 7). All concentration values on both panels are normalized to the case of a minimally vegetated world with global average Damköhler coefficient values (ω = 1.5, Dw = 0.03 m/yr). In the left panel of figure 7, showing our results without plant recycling feedbacks, average soil weathering concentrations are primarily set by the Dw values (x-axis), with minimal sensitivity to ω values. Values increase slightly with increasing ω due to flattened runoff gradients (fig. 5D) leading to higher weathering concentrations closer to chemostatic concentrations.
Inclusion of plant recycling feedbacks dampens the response to Dw values and amplifies the model response to ω. Moving from point A to point A' in figure 7 represents addition of plant recycling feedbacks only (effects a and b with introduction of the recycling term, f, in eq. 7) and results in a 1.37-fold increase in weathering concentration (see also table 3). Moving from point A' to B' to C' on figure 7 (middle panel) represents increasingly efficient recycling as ω increases. Further, the effect of increasing recycling and precipitation rates on global-scale recycling feedbacks, along with local-scale feedbacks increasing soil pCO2 (fig. 6C), results in a slight decrease in the effective Dw values due to increases in Cmax (fig. 6D). This is because Cmax also appears in the denominator of Dw (see eqs. 5–7), which is represented by the trajectory of points A to A” to B” to C” in the right panel of figure 6, where [HCO3−] concentrations increase 2.39 (point B”) and 2.55 (point C”) fold. If the effect of Cmax on the denominator of Dw is not accounted for, this would lead to minimal overestimation (a few percent) of the weathering concentration increases (see contours of Dw and ω relative to points B” and C”). However, if plants also increased the length scale (deeper soils) represented by increased Dw values, these weathering concentration changes may be greater. For example, for modern mountain Dw values of 0.3 m/yr (Maher and Chamberlain, 2014), vascular plant-world weathering solute concentrations would have increased by a factor of ∼2.8, with angiosperm-world concentrations increasing to a factor of ∼3.0 (all relative to the minimally vegetated world; fig. 7, right panel). Thus, the magnitude of increased weathering solute concentrations depends upon whether these changes occur in low-relief cratonic areas versus highly erosional mountainous regimes, the paleogeography of those mountainous regions, and exposed continental area as influenced by sea level or continental freeboard changes. Increased inland vapor recycling increases the likelihood that precipitation reaches inland cratonic mountain ranges. Regardless of topography, each transition in plant evolution increased the thermodynamic potential for solute generation by increasing soil pCO2 (Cmax), increased vapor recycling (f), and possibly by deepening the soil weathering length scale (Dw).
Uncertainty of the Calculations
There is significant potential uncertainty in our calculations related to both changes in the Budyko landscape parameter (ω values) and potential effects due to temperature on the maximum weathering concentration (Cmax). The fits to the climate model fluxes shown in figure 2 have an uncertainty (standard error for ω values) of ± 0.02 to 0.09. The maximum uncertainty value of 0.09 propagates to a ∼10 percent uncertainty in the calculated concentration changes (points A', B' and C') when only effects a and b (the hydrologic effects) are accounted for (see blue dots in fig. A6).
Uncertainty due to temperature effects on the thermodynamics of the net weathering reaction are explored in figure 3A. While the slope does not change, the absolute value of the Cmax values increase with lower – or decrease with higher – temperatures, as recently shown using similar reactive transport simulations by Winnick and Maher (2018). Because the thermodynamic rate constants are non-linear for the net weathering reactions (see Maher and Chamberlain, 2014), the uncertainty is asymmetric. Propagating the range of Cmax values determined for 5 and 25 °C, relative to the GWB simulations at 15 °C, give an uncertainty of 31 percent and 15 percent, respectively, when accounting for all effects (both thermodynamic and dynamic) (see green dots in fig. A6; points A”, B”, and C”). As such, our calculations are more uncertain at colder temperatures (lower atmospheric CO2).
Analysis of the cumulative uncertainties associated with only hydrologic effects and both hydrologic and thermodynamic effects (fig. A6), shows that our vascular plant- and angiosperm-worlds are significantly different from a minimally vegetated world (even with the inclusion of all recycling effects). However, for changes due to only hydrologic effects (blue dots; points B' and C'), the change in weathering solute concentration is not significantly different between our vascular plant and angiosperm-worlds because the Budyko ω values are relatively similar (2.2 and 2.5, respectively), leading to similar runoff distributions (fig. 5D). Nonetheless, when all effects are considered (green dots; points B” and C”) changes in weathering concentration are distinct between our vascular plant- and angiosperm-world concentrations, with the lower bound of the angiosperm concentration (C”) approximately equal to the average value for the vascular plant concentration (B”).
DISCUSSION
Plant Evolution and Weathering
As illustrated by the modeling results presented in this study, the evolution and spread of land plants enhanced chemical weathering concentrations through two complementary mechanisms: (1) increasing the down-wind transport and recycling of precipitation; and, (2) increasing the thermodynamic maximum of solutes (Cmax) produced from weathering via an increase in soil pCO2. The ability of a continental surface to generate higher solute concentrations means that the weathering flux should be higher at any given runoff. However, the globally averaged silicate weathering flux must equal the volcanic and metamorphic flux of carbon on geologic (that is, 1 Ma) timescales (Berner and Caldeira, 1997; Broecker and Sanyal, 1998; Zeebe and Caldeira, 2008; Caves and others, 2016a). Consequently, a world populated with plants in which the efficiency of weathering increases requires a reduction in other factors which contribute to weathering, such as atmospheric CO2 concentrations, temperature, and runoff. Consequently, each step in plant evolution may permit a lower steady-state atmospheric pCO2, resulting in lower steady-state temperature and therefore runoff. Further, that weathering solute concentrations are higher in a world with plants indicates that sudden increases in runoff—driven perhaps by large carbon cycle perturbations—will produce correspondingly larger increases in the weathered solute flux and thus a more rapid reduction in temperature and runoff. Thus, we suggest that the strength of the silicate weathering feedback increased as plants colonized the continents and the largest increase is the result of deep rooting vascular plants in the late Devonian and Carboniferous, with a somewhat smaller increase during the expansion of angiosperms in the Cretaceous and Cenozoic.
The effect of plants on chemical weathering were qualitatively outlined by Berner (1992) and Drever (1994) and quantitatively incorporated into the GEOCARB family of models, most recently GEOCARBSULF (Berner, 2006; see also more detailed model description in Berner and Kothavala, 2001 and Royer and others, 2014, as well as incorporation into COPSE in Bergman and others, 2004). To incorporate the many possible mechanisms for enhanced weathering by plants (Berner, 2006), these models treat the effect of plants on weathering such that: 1) pre-Devonian weathering rates are 25 percent as sensitive to changes in pCO2 compared to modern weathering rates; 2) the effect of gymnosperms begins at 380 Ma and ramps up to full vascular plant effect at 350 Ma with weathering rates 87.5 percent as sensitive compared to modern; and 3) angiosperm weathering sensitivity is the same as modern weathering rate sensitivity and ramps up over the period from 130 to 80 Ma (see Berner and Kothavala, 2001; Royer and others, 2014; summarized in table 3). Of similar magnitude, steady-state mass balance calculations from the GEOCLIM model (Donnadieu and others, 2006, 2009; Le Hir and others, 2011) suggest increases in silicate weathering in the late Devonian decreased atmospheric pCO2 by a factor of 1.55. Additional effects of plant weathering have often been incorporated by a CO2 fertilization of productivity term that averages across the variability in response seen in the modern vegetation (Mooney and others, 1991; Berner, 1994; Banwart and others, 2009), although applicability to the non-angiosperm fossil record has been questioned (Boyce and Zwieniecki 2012). However, CO2 fertilization of productivity is not considered in the modeling performed here.
In a broad sense our results are similar to previous estimates, with some key differences. Before discussing these differences, it is important to note that our estimates of weathering rate changes depend upon: 1) where on the landscape weathering occurs (for example, mountains versus lowlands, which influences Dw); 2) the estimated Budyko dimensionless landscape parameter (ω) in equation 4, which determines the conversion efficiency of P to ET for vascular plant and angiosperm worlds; 3) whether water is recycled through the entire weathering zone or only part of it; and, 4) the spatial extent of deep rooting plant cover.
The enhanced movement of water across the landscape due solely to increased transpiration has minimal impact on weathering directly; the key is the repeated interactions with the substrate that this transpirative recycling affords (represented by moving vertically in the fig. 7 panels). Water vapor can be recycled such that water does not immediately enter runoff and instead can enter the weathering zone again through transpirative return of water to the atmosphere if rooting depth of plants occurs near (or below) the depth of the weathering front (for example, Stone and Kalisz, 1991). Our results support past recognition that inland transpiration recycling can even decrease total runoff amount overall as water is lost as vapor off the downwind side of the continent (Kleidon and others, 2000; Alkama and others, 2012) (fig. 5D). With our preferred boundary conditions (constant flux to no gradient) and calibrated ω values (fig. 2), runoff decreases by 17 percent and 23 percent for the vascular plant- and angiosperm-worlds, respectively, relative to the minimally vegetated scenario (fig. 5D, integral of green and blue curves relative to red curve).
Many individual plants have relatively shallow roots. As such, the focus here on water being entrained into transpiration only after it has already participated in weathering is an end-member case, but a reasonable one. The critical part of the weathering zone—where the most solute is generated, the most CO2 is converted to alkalinity, and primary minerals are depleted and secondary minerals precipitated—tends to be in the zone of seasonal water table fluctuations, particularly in low-lying areas (Stone and Kalisz, 1991; Ehleringer and Dawson, 1992; Jongmans and others, 1997; Landeweert and others, 2001; Lucas, 2001). Although plant roots tend to stop at the water table, they can draw substantial proportions of their water from below this interface without substantially penetrating it (Fan, 2015). Indeed, stable isotopic studies and studies of modern plants show that in many cases, trees tap deep-water sources that should be well below the weathering front (upper few meters of soil) (Stone and Kalisz, 1991; Ehleringer and Dawson, 1992; Jongmans and others, 1997; Kleidon and Heimann, 2000; Landeweert and others, 2001; Lucas, 2001). Thus, about a quarter of the modern global transpiration stream (Evaristo and McDonnell, 2017) is likely to come from below the zone of maximal weathering. The first vascular land plants of the Silurian and Early Devonian would have been exclusively shallow rooting, but the trees of the later Devonian and Carboniferous had deeper roots that would have been able to access the water table in many environments, with water table interactions directly recorded in some rooting fossils (Algeo and others, 2001; Pfefferkorn and Wang, 2009; Boyce and others, 2017).
Given these caveats, we can place some quantitative constraints on the role that plants play in increasing weathering fluxes through time (figs. 7 and 8). For the following estimates, we use our end-members for a no-transpiration world (minimally vegetated with no deep-rooted vegetation), vascular plant world, and angiosperm world, global average Dw values (0.03 m/yr) from modern constraints (Maher and Chamberlain, 2014), and also include both a precipitation – soil CO2 – Cmax relationship and the introduction of a recycling term into the weathering equations. We conclude that moving water vapor inland due to transpiration increases solute concentrations by a factor of 2.39 for vascular plants and by a factor of 2.55 for angiosperms (table 3). This effect was noted by Berner (1992) and Drever (1994) based on climate model results by Shukla and Mintz (1982), but is not modeled explicitly in GEOCARBSULF (although there is a function incorporated into GEOCARBSULF between increased runoff and increased CO2 levels; see Berner and Kothavala, 2001; Royer and others, 2014) or more specific process-based models of plant-enhanced weathering (for example, Banwart and others, 2009; Taylor and others, 2009, 2011, 2012; Quirk and others, 2012; Beerling and others, 2012; Thorley and others, 2015). Climate models comparing vegetated versus desert worlds have shown that plant transpiration can double precipitation over land (Kleidon and others, 2000).
Sensitivity of weathering fluxes (Fsil) to perturbations in climate, through changes in runoff (q), using framework presented in this paper (after Maher and Chamberlain, 2014; von Blanckenburg and others, 2015; Ibarra and others, 2016). All lines are plotted for the global average Damköhler coefficient (Dw = 0.03 m/yr) adjusted for increased Cmax (Dwcorr) as in figure 7 and plotted as normalized fluxes (see Appendix 3 for derivation of equations). The red lower curve is equivalent to point A in figure 6. The green and blue lines represent points B” and C” in figure 7 using flux averaged recycling factor (favg) and Cmax values, representing vascular plants and angiosperm curves respectively.
Using climate models, Lee and Boyce (2010) (see also Boyce and Lee, 2010, 2011, 2017) have shown that angiosperms—which have a greater transpirative capacity than other vascular-plants—cause an increase in precipitation over continents, although changes in the spatial pattern of precipitation are complex. Moreover, vegetation can increase precipitation over landmasses independent of transpiration via albedo effects (see desert world vs. no-transpiration world maps in Boyce and Lee, 2017). Thus, one would expect that both the advent of pigmented microbial soil communities, followed by the evolution of vascular plants and of angiosperms, played a significant role in increasing weatherability by altering climate. Although this effect is difficult to quantify even with climate models, any long-term carbon cycle model (for example models by, Berner and Kothavala, 2001; Bergman and others, 2004; Berner, 2006; Zeebe, 2012; Mills and others, 2014) must account for this effect.
Increases in precipitation over land due to increases in transpiration must have occurred at least twice in Earth's history: once at the advent of vascular plants, particularly deep-rooted trees in the later Devonian and Carboniferous, and again with the evolution of flowering plants with high transpiration capacities in the latest Cretaceous and early Cenozoic. Transpirative effects on weathering were not likely to be important before this time; for example, during the evolution of the first land plants in the Ordovician through early Devonian, bryophyte-grade and simple shallow-rooting vascular land plants had limited capacity to remobilize soil water in any significant manner (Boyce and Lee, 2017).
Finally, we conclude that angiosperms increased solute concentrations by up to 7 to 55 percent over other vascular plants (point B' or B” to C' or C” in Figure 7). The range is due to the unconstrained transition and trajectory between two likely end-members, a weathering system that is independent from plant biomass (represented by A'), to the current state (C”) where biomass and soil CO2 are correlated to precipitation. The first trajectory is that the transition between both scenarios occur early, linked to the appearance of roots, allowing for colonization of well drained soils and a root system that provides decaying biomass directly to the source (in this case B” over B' and a 7 percent change in solute concentration between B” and C”). Alternatively, a second trajectory is supported if the modern relationship between soil CO2 and precipitation is an angiosperm-only trait, providing a higher flux of decaying biomass to the system for a given precipitation and in that case a transition from B' (or more likely an intermediate value between B' and B”) to C” would occur with the associated 55 percent increase in solute concentration. When compared to other vascular and non-vascular plants, angiosperm leaves have shorter life spans, higher nutrient contents and are paired with higher leaf decay rates (Cornwell and others, 2008), though this difference in decomposability is not necessarily observed in wood or root traits (McCormack and others, 2012), and does not increase soil nitrogen availability (Mueller and others, 2012), leaving it hard to constrain whether the well supported differences in leaf traits play a large role in determining past soil CO2.
Nevertheless, our results are similar in range and magnitude compared to the output of GEOCARBSULF (Berner, 2006; see also Berner and Kothavala, 2001), where an increase of about 12.5 percent in weathering rates are observed with the transition to angiosperm dominance. This increase due to angiosperm dominance of ecosystems is unlikely to explain the long-term decrease in Cenozoic pCO2. Though Caves and others (2016a) demonstrated that the strength of the silicate weathering feedback increased during the mid to late Cenozoic, angiosperm ecosystems were well-established by the latest Paleocene, indicating that the reason for long-term declining Cenozoic pCO2 cannot be attributed to our estimated increase in weathering rates due to the rise of angiosperms. Given that the angiosperm proliferation was in the Late Cretaceous and early Cenozoic and they were already overwhelmingly dominant by the latest Paleogene, another explanation would be needed for further increases in the Neogene. The Neogene spread of grasslands, and associated potential albedo changes—a topic not modeled here—would be an obvious plant-related candidate (Freeman and Colarusso, 2001; Fox and Koch, 2004; Feakins and others, 2005; Edwards and others, 2010; Strömberg, 2011; Chamberlain and others, 2014), though other non-biological factors such as tectonics (for example, Misra and Froelich, 2012; Caves and others, 2016a) are possible.
Though plant evolution has caused a long-term, unidirectional change in the global weathering regime, superimposed on this biological effect are transient climatic (less than 10 Myr) disturbances, such as weathering-erosion perturbations hypothesized for the end-Permian and end-Cretaceous extinction events (see for example McElwain and Punyasena, 2007; Sheldon and others, 2014). Furthermore, the effects described here related to land plant evolution may be overprinted by the longer-term changes in the composition of subaerial continental crust, changes in exposed lithologies, and tectonically induced mountain building, which also influences changes in soil depth and weathering zone thickness (for example, Ronov and others, 1980; Bluth and Kump, 1991; Gibbs and others, 1999; Berner, 2006; Li and Elderfield, 2013).
Taken together, we suggest that weathering concentrations increased due to the rise of plants from a minimally vegetated world to the modern by more than 2-fold (point A to B” or C” in fig. 7) due to the combined effects of increases in movement of water vapor by transpiration, in soil CO2 concentrations, and in thickness of the weathering horizon. This result implies a similar but slightly weaker weathering sensitivity, than the average 4-fold increased used in GEOCARBSULF (Berner, 2006; see also Berner and Kothavala, 2001; Royer and others, 2014; table 3).
Climate Implications
Much of the above discussion focused on the mechanisms—including higher soil pCO2 and downwind water vapor recycling—by which plant evolution has altered silicate weathering. These mechanisms would have altered the effectiveness of silicate weathering in modulating climate. Here, we discuss how these mechanisms—coupled with the concept of a silicate weathering feedback—likely led to long-term, dramatic changes for the Earth surface.
Though a number of different formulations for the silicate weathering feedback have been proposed (Walker and others, 1981; Berner and others, 1983; Volk, 1987, 1989; Berner and Canfield, 1989; Berner and Kothavala, 2001; Bergman and others, 2004; Berner, 2006; Uchikawa and Zeebe, 2008; Caves and others, 2016a), central to all of these formulations is the idea that the silicate weathering flux is sensitive to climate. In these cases, climate usually refers to temperature, runoff, or pCO2, such that higher temperatures, greater runoff, or higher pCO2 results in a greater silicate weathering flux. Here, we have utilized the formulation of Maher and Chamberlain (2014) (eqs. 5 and 7), which directly links the concentration of silicate-derived solutes to climate through runoff.
The silicate weathering flux is calculated using eq. 1 (reproduced below):
where C is calculated using eq. 7 (reproduced below):
where we have proposed that f = 1 prior to plant evolution. The coupling between Fsil and climate (q) (that is, strength of the weathering feedback) is modified by the values of Cmax, f (recycling factor), and Dw. In figure 8, we plot a normalized silicate weathering feedback curve, as represented by Fsil as a function of q, prior to plant evolution (lowest red line). In this figure, the silicate weathering feedback is normalized to the value of q required to produce a normalized Fsil value of 1 (Fsil0). Because Fsil must balance the volcanic and metamorphic inputs of carbon on geologic timescales (> 1 Ma), the y-axis is equivalent to the normalized flux of volcanic and metamorphic carbon.
The evolution of vascular plants in the late Devonian increased Cmax and f; as a consequence, Fsil must have increased and then subsequently decreased to match the volcanic and metamorphic flux. Increased Fsil would remove CO2 from the atmosphere, resulting in declining pCO2 and hence cooling, which would result in lower q. Though the volcanic/solid Earth inputs of carbon may have changed through time (McKenzie and others, 2016), there is no evidence that this flux changed as a consequence of plant evolution. Thus, across the late Devonian, we assume the volcanic and metamorphic inputs remained constant and, thus, Fsil must eventually decline following the spread of vascular plants. This decline in Fsil is accomplished because q decreases in response to cooling even as Cmax and f remain higher than prior to the spread of vascular plants. Consequently, the silicate weathering feedback would be characterized by the middle green line in figure 8, which has a higher slope than the pre-vascular plant silicate weathering feedback. Moreover, following the spread of vascular plants, any carbon cycle perturbation would result in a larger increase in Fsil than prior to the evolution of plants, thus dampening the climatic impact of such a perturbation. In this sense, the strength of the silicate weathering feedback increases as plants evolve. A similar increase in the strength of the silicate weathering feedback would occur when the evolution of angiosperms further increased Cmax and f (highest blue line in fig. 8). Thus, the evolution of plants resulted in a long-term decline in the steady-state atmospheric pCO2, accompanied by decreases in surface temperatures and in runoff—a result also found by others who have modeled the impacts of plant evolution on climate (Berner, 2006; Le Hir and others, 2011).
A major assumption in this analysis is that both Cmax and f remain largely unchanged even as atmospheric pCO2 decreased. In reality, Cmax is influenced by soil pCO2, which is at least partially influenced by atmospheric pCO2 (Cerling, 1984). Thus, although soil respiration will be a larger influence after vascular plant evolution (Berner, 1992), Cmax is still likely to co-vary positively with atmospheric pCO2. Changes in f are harder to evaluate and ultimately would require a more detailed consideration of the plant biology and of soil depth changes that might have their own complex feedbacks (for example, Sellers and others, 1996; Betts and others, 1997; Swann and others, 2016). For instance, the fraction of precipitation that is partitioned to transpiration would be expected to increase with lower atmospheric pCO2, as more gas exchange would be required to acquire a given amount of carbon (for example, Ainsworth and Rogers, 2007; Dalling and others, 2016). Conversely, widespread observations of “greening” (for example, greater leaf-area; see also Skinner and others, 2017) with higher pCO2—both in the modern (Zhu and others, 2016) and in the geologic past (Caves and others, 2016b)—suggest that global transpiration fluxes may scale with increasing atmospheric pCO2. Quantifying these additional feedbacks on weathering is beyond the scope of the current study, though they represent unanswered, critical Earth system questions. We note, however, that the modeled changes in Cmax and f due to plant evolution are likely greater than those induced by feedbacks with changes to long-term atmospheric pCO2.
The evolution of plants would have thus wrought major changes to the Earth system simply through its impact on silicate weathering. In short, the evolution of plants may have dramatically cooled the climate and reduced runoff normalized to precipitation. At this point, the consequent impacts remain speculative, but likely range from affecting how rivers transported sediment to atmospheric circulation changes (for example, Davies and Gibling, 2003, 2010, 2010; Boyce and Lee, 2010, 2011). Additional quantitative constraints on the impact of plant evolution on silicate weathering, as well as accounting for a richer variety of processes than our simple model permits, will yield further insights into the evolution of our planet's carbon cycle.
CONCLUSION
In this contribution, we used modern empirical and modeling constraints on atmospheric vapor-transport over continents and terrestrial weathering fluxes. We found that the emergence of widespread transpiration and downstream vapor recycling, and subsequent increases in soil pCO2, would have more than doubled weathering concentrations when deep-rooted vascular plants (Devonian-Carboniferous) fully replaced a minimally vegetated (pre-Devonian) world (fig. 6; table 3). Further, we found that the domination of ecosystems by angiosperms (late Cretaceous and Cenozoic) would have further increased weathering concentrations by an additional 7 to 55 percent. Through this work, we show that incorporating inland vapor recycling—a factor that is not currently incorporated into existing 1-D carbon cycle models—results in a significant change in weathering solute concentrations. This inland vapor recycling increases weathering solute concentrations due to: (1) increased total precipitation, which increases soil pCO2 and hence theoretical maximum aqueous solute concentrations, and (2) decreased total runoff that drive the weathering system closer to chemostatic concentrations. This inland vapor recycling would have also increased the likelihood that precipitation reaches inland, cratonic mountain ranges. Partitioning of the vapor recycling effects and feedbacks on weathering thermodynamics indicate that 55 percent of the increase in weathering solute concentrations from a minimally vegetated world to an angiosperm-dominated world is attributed to the thermodynamic effects of increased soil pCO2, which decreases soil pH and increases equilibrium solute concentrations. While the absolute values of the modeling results carried out in this contribution certainly varied over the Phanerozoic due to long-term changes in erosive regime, volcanic degassing rates, super continent formation, and changes in the solar constant, the trends observed within this modeling framework are indicative of the directionality of change.
Our model results place previously hypothesized effects of plant evolution on the Earth system on a more mechanistic footing. We show that deep-rooted plants can act to modulate weathering fluxes through controlling the hydrologic cycle. To maintain carbon mass balance in vegetated worlds, the higher weatherability of the continental surface must be balanced by a globally averaged lower runoff value. Further, steady-state atmospheric pCO2 is lower for a more weatherable planet, matching model results following the rise of land plants (for example, Berner, 2006; Le Hir and others, 2011). Finally, the primary implication of this work is that the strength of the silicate weathering feedback increased as plants colonized the continents, in line with observations that carbon isotope excursions decrease in magnitude over the Phanerozoic (Bachan and others, 2017). Previous models of the impact of plant evolution on biogeochemical cycles have overlooked the importance of inland vapor recycling in setting the strength of this feedback. Thus, we propose that the largest increase in the silicate weathering feedback sensitivity was caused by the deep rooting vascular plants in the late Devonian and Carboniferous, with a smaller feedback sensitivity increase occurring during the expansion of angiosperms in the Cretaceous and Cenozoic.
ACKNOWLEDGMENTS
We thank Nathan D. Sheldon and Yves Godderis for thorough reviews and Associate Editor Christopher J. Poulsen for handling the manuscript. This manuscript benefited from critical discussions with Kate Maher, Matthew J. Winnick, Friedhelm von Blanckenburg, Isabel Montañez, and P. James Dennedy-Frank. This study was initiated with additional input from Thomas Boag, Danielle Moragne, Xiaowei Li, Maria Rugenstein, and Matthew Nelsen. Jeremy Caves Rugenstein is funded by an ETH Fellowship. Daniel Ibarra is funded by a Stanford EDGE-STEM Fellowship. Kimberly Lau is funded by an Agouron Geobiology Postdoctoral Fellowship and an ARCS Foundation Fellowship. Dana Thomas is funded by a Stanford Center for Carbon Storage Fellowship. This is EHGoS contribution #648.
APPENDIX 1: MODEL 1 DERIVATION
DERIVATION
The mass balance equations.
In the following we derive an advection-diffusion-reaction model simulating the transport of water vapor across a 1-dimensional continental domain. The advection term accounts for directional transport (prevailing winds), diffusion simulates turbulent mixing, and reaction terms account for moisture that is lost to precipitation and returned to the air mass via evapotranspiration. Within a domain of a given length we consider a control volume with dimensions of dx (width) by A (cross-sectional area). The mass balance for the control volume can be expressed as: the time rate of change of mass (TROCM) equals mass-rate-in (MRI) minus mass-rate-out (MRO) plus sources or sinks (S). Or:
where the time (t)-rate of change of water vapor (w) in a grid cell is the volume of the cell times its vapor content:
and the total water vapor flux per unit time into the control volume (MRI) is the sum of the advective (Qa) and diffusive fluxes (Qd):
Each flux is the product of the area-normalized flux and the cross-sectional area:
The inputs (at point x) are known but the outputs (at point x + dx) have to be estimated. We do so by expressing them as a first order Taylor series expansion of the input fluxes in the x direction.
Finally, the reaction terms (S/S), precipitation (P), evaporation (E), and transpiration (T), similarly are considered to scale with the volume of the control volume:
Putting these equations together and canceling terms:
Assuming that A and dx are constant across x we can divide through by Adx:
Where the diffusional flux is proportional to the water vapor gradient:
And the advective input flux is the transport velocity (u) times the water vapor content:
Substituting the flux relationships for qa and qd, and assuming that D and w are independent of x (allowing us to move them out of the differential and rearrange) we get a 2nd-order ODE:
Finally, we solve the system in steady state, so the time derivative is set to zero:
Boundary conditions.
A 2nd-order ODE requires two boundary conditions, one on each side. We assume that the source of moisture is always on the left side of the domain (west side of the continent) while the right side represents the interior or edge of the continent. Each boundary (left or right) can be defined by one of three types of conditions, making a total of 9 different boundary condition combinations. Not all combinations are physically reasonable.
First, we consider a constant flux condition on the left. In this case both the advective and diffusive fluxes sum up to a constant prescribed value so that there is constant net transport of material across the edge of the domain. On the right, all three boundary conditions make physical sense in combination with a constant flux condition on the left.
Possible combinations of constant flux, constant gradient, and constant gradient boundary conditions on the left and right sides. ‘x’ defines physically reasonable boundary conditions.
A constant flux condition on the right with a value of zero means that no moisture crosses the far edge of the domain. This combination of boundary conditions allows for testing of the accuracy of the code, since what comes in must come out. A value other than zero allows a prescribed amount of moisture to leave the domain. This boundary condition combination is shown in Figure A5.
A constant gradient condition on the right with a value of zero enforces a zero gradient condition across the edge of the domain, which effectively means zero diffusive flux. This condition corresponds to the scenario that a balance had been reached between precipitation, evaporation, and transpiration such that no additional gradient in water vapor content remains. Advection then sweeps the remaining moisture off the edge of the domain (that is, the continental mass).
Constant concentration condition with a zero value on the right enforces the condition of zero water vapor on the boundary and thus no advection off the edge.
2. Second, a constant concentration condition on the left with a positive value implies that our domain is connected to a much larger volume of air, so that upon removal of moisture the vapor content of the larger domain remains unchanged. This condition makes physical sense with a constant concentration or a constant gradient condition on the right. A constant flux condition only makes sense if diffusion is more dominant than advection (otherwise it leads to piling up of water vapor on the right side of the domain).
3. Lastly, a constant gradient condition on the left makes little physical sense.
NUMERICAL SOLUTION
Discretization
The advection-diffusion-reaction equation can be written as:
We can integrate over a single cell spanning xi–1/2 to xi+1/2 and obtain the mass balance for it:
Where the spatial integral over the time derivative, 
dx, on the left side is denoted by w̄1, and the integral over the sources and sinks for the cell, 
, by s̄1.
The expression simply states that the change in a mass of material in a cell is the difference between the sum of advective and diffusive fluxes coming in on the left minus those leaving the cell on the right:
For the advective flux, we estimate the value at the boundary by linear interpolation:
For diffusion, we estimate the derivative by a first order difference:
Substituting the quantities yields:
Equivalently one can use first-ordered centered in space discretization for both the first and second order terms:
And replace them in the equation steady state equation, which yields the same result:
To simplify we denote the constant coefficients with the letter C:
Where
Boundary Condition Implementation
Constant flux conditions.
In the finite volume method, constant flux boundary conditions that define a total flux are naturally resolved, and require no discretization. Specifically, the condition is that the total flux on the right side of the domain is given by:
So within the context of the finite volume balance:
Using the approximations for the advective and diffusive fluxes:
Similarly, on the right boundary of the domain at location ‘n’:
Constant concentration conditions.
The implementation of constant concentration condition is simple. In this case, the value of w in the first (or last) cell is known. So:
Thus, we set the corresponding coefficients such that all terms other than w1 vanish and add - winin our boundary conditions vector.
Constant gradient conditions.
In the case of constant gradient boundary conditions the diffusional flux across the edge of the domain is known:
Discretizing for the right edge:
Substituting coefficients:
Now discretizing for the left side:
Which using the coefficients goes to:
MATRIX SETUP AND INCORPORATION OF BUDYKO
As an example using the two constant flux boundary conditions, we discretize our equations (in this example with 6 cells), which gives the following matrix:
We use the Budyko relationship to define how precipitation relates to the sum of evaporation and transpiration. Evaporation and transpiration are not partitioned explicitly in this modeling framework. We define precipitation to be a linear function of the water vapor content:
And define the evapotranspiration (E+T=ET) as a function of both precipitation and potential evapotranspiration (Ep) following the formulation used by Greve and others (2015):
When ω=1, ET=0. As the source (ET) and sink (P) terms are non-linear functions of vapor content (w), we use a root-finding procedure to simultaneously solve for the unknowns in the above matrix (wi), rather than a matrix inversion method. We implement the model using the multiroot function in the R package ‘rootSolve’ (Soetaert, 2016), which uses a numerical Newton-Raphson method to find the roots of the above equations.
APPENDIX 2: COMPILATION OF MODERN TRANSPIRATION MODIFICATION STUDIES
Treatment sites from Bosch and Hewlett (1982) (see their Table 1), and Zhang and others (1999) used to test the effect of reduced canopy transpiration on Budyko relationships (fig. 2B) with 100% treatment (full clear-cutting), all P, q and Ep units are mm/yr. Ep values were calculated as estimated via spatial interpolation from Trabucco and Zomer (2009)
APPENDIX 3: DERIVATION OF NORMALIZED WEATHERING EQUATIONS
To show the sensitivity of the weathering equations, we normalize the equations from Maher and Chamberlain (2014) following von Blanckenburg and others (2015) (their Figure 3). The solute production model predicts concentration (C) as a function of runoff, with assumed values for Damköhler coefficient and maximum concentration (Cmax), with a scaling factor (τ = e2) (equation 4 from the main text):
Area normalized weathering flux (Fsil) is given by:
Normalizing this equation by a baseline weathering flux (Fsil0) at a given baseline runoff value (q0) gives:
Note that this assumes no change in Cmax, which cancels. If Cmax does not cancel the normalized equation is:
In figure 8, we use the above relationships to plot normalized fluxes for scenarios A, B” and C” of figure 7. For points B” and C” we also included the recycling term (f) as a multiplier before Dw (see text for details).
APPENDIX 4: SOIL CO2–PRECIPITATION RELATIONSHIP
Relationship between soil pCO2 and P, calculated as a linear fit to the datasets from Cotton and Sheldon (2012) (blue points) and Cotton and others (2013) (gray points).
APPENDIX 5: ALTERNATE BOUNDARY CONDITION EXAMPLE
Water and solute fluxes simulated using constant flux to constant flux boundary conditions (see text for explanation). Simulations are shown for ω values representing the simulations from Boyce and Lee (2010) and modern observations as in figure 2. Model 1: (A) Total water vapor content, (B) precipitation, (C) evapotranspiration, and (D) runoff. The flux out on the right (eastern) boundary condition was set to 40% of the incoming water vapor, comparable to the Angiosperm simulations shown in figure 5. Model 2: (E) soil pCO2; and (F) solute [HCO3−] calculated including both hydrologic and thermodynamic feedbacks.
APPENDIX 6: SENSITIVITY OF CALCULATIONS
Sensitivity of normalized [HCO3−] calculations for points discussed in the text (A'-C' and A”-C”). The dashed line corresponds to point A. Blue dots show uncertainty (10%) from the Budyko omega parameter propagation of hydrologic effects only (A', B' and C' in fig. 7). Green dots show total uncertainty (+31%/-15%) for propagation of thermodynamic temperature effects and hydrologic effects (A”, B” and C” in fig. 7).
REFERENCES
In this issue
Jump to section
- Article
- Abstract
- INTRODUCTION
- MODELING FRAMEWORK
- MODEL RESULTS
- DISCUSSION
- CONCLUSION
- ACKNOWLEDGMENTS
- APPENDIX 1: MODEL 1 DERIVATION
- APPENDIX 2: COMPILATION OF MODERN TRANSPIRATION MODIFICATION STUDIES
- APPENDIX 3: DERIVATION OF NORMALIZED WEATHERING EQUATIONS
- APPENDIX 4: SOIL CO2–PRECIPITATION RELATIONSHIP
- APPENDIX 5: ALTERNATE BOUNDARY CONDITION EXAMPLE
- APPENDIX 6: SENSITIVITY OF CALCULATIONS
- REFERENCES
- Figures & Data
- Info & Metrics
- References