Wetness controls on global chemical weathering

The formation of soils, the evolution of the biosphere, and the CO2 content in the atmosphere are strongly impacted by chemical weathering. Due to its manifold importance for the long-term stability of the Critical Zone, it is crucial to link weathering rates to the environmental conditions affecting it and develop accurate rate laws for landscape evolution and carbon cycle modeling. Here we use the π theorem of dimensional analysis to provide a theoretical framework to global datasets of weathering rates. As a result, a strong relation between chemical depletion, precipitation and potential evapotranspiration synthesizes the primary role of wetness. Based on this finding, we estimate the spatial distribution of chemical depletion fraction and find that, globally, soils are 50% chemically depleted, 61% of the land is in kinetic-limited conditions, while only 1% is supply-limited. The remaining 38% of the land is in a transitional regime and susceptible to changes in wetness.


Introduction
Chemical alteration of the parent material begins when tectonic and topographic stresses fracture the rock and expose minerals to the action of meteoric water (Clair et al 2015, Anderson 2015, Moon et al 2017, Eppes and Keanini 2017, Holbrook et al 2019. From that moment, chemical weathering becomes part of an intertwined system of feedbacks between the Critical Zone and the atmosphere (Gaillardet et al 1999, Kump et al 2000, White and Buss 2014, Pelak et al 2016, Braun et al 2016. As minerals dissolve, nutrients are released, clay minerals precipitate, and soils form, favoring the development of plants and microbial activity (Brady and Weil 2008, Amundson et al 2007, Chorover et al 2011. This in turn increases CO 2 production in the upper soil layers (Singh andGupta 1977, Hanson et al 2000), from which dissolved CO 2 is transported to the deeper soil layers (weathering zone) enhancing chemical weathering (Calabrese et al 2017). Weathering reactions then consume the incoming CO 2 , subtracting it from the atmosphere and mitigating the greenhouse gas effect. An increase in atmospheric temperature and CO 2 content in fact would enhance weathering, triggering the negative feedback that stabilizes atmospheric conditions (Berner 1994, Winnick andMaher 2018).
To understand large scale controls on weathering and find empirical laws relating its rate to precipication, runoff, temperature, and erosion, weathering rates have been measured globally from river chemistry (Gaillardet et al 1999, White and Blum 1995, Dessert et al 2003, Oliva et al 2003, West et al 2005 and from chemical depletion profiles of the regolith (April et al 1986, Brimhall et al 1992, Riebe et al 2003, Rasmussen et al 2011, Dere et al 2013, Ferrier et al 2016, Dixon et al 2016. Although good correlation between weathering rates and these factors is often observed, the lack of consensus on the functional form of the rate laws and the variability in their parameterization (e.g., apparent activation energy) ( Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Dixon et al 2009, Larsen et al 2014, Dixon et al 2016, for instance on the relationship between chemical weathering and erosion. Studies based on reactive transport modeling have proved useful for elucidating the mechanisms driving mineral dissolution at the Darcy scale (Maher 2010, Maher and Chamberlain 2014, Li et al 2014, Lebedeva and Brantley 2018, Maher and Navarre-Sitchler 2019, where the Damköhler number, which compares reaction and fluid timescales, controls the reaction. However, due to the complex heterogeneous structure of the subsurface, idealized conditions (e.g., homogeneous mineral distribution and constant vertical fluid velocity) are typically assumed, so that results remain mostly applicable to a laboratory rather than a field setting, or to a specific site where detailed information on the system is available. Furthermore, the upscaling of the reactive transport equations to soil column or watershed scales remains arduous, due to the mathematical complexity of the system of coupled partial differential equations describing the reactive transport problem (Lasaga 1984, Cederberg et al 1985, Yeh and Tripathi 1991, Saaltink et al 1998, Battiato et al 2009, Battiato and Tartakovsky 2011. To help resolve the complex interaction between chemical weathering and climate, here we provide a theoretical structure to interpret weathering rate data based on the π theorem of dimensional analysis. In complex problems involving a large number of variables, the π theorem is often of essential help to derive underlying physical law from experimental data (Barenblatt and Isaakovich 1996). It reduces the number of independent variables by grouping them into dimensionless groups and it allows to obtain scalable laws, which are applicable across sites with very different climatic and geologic contexts. By applying it we show that the chemical depletion fraction-a measure of weathering rate relative to the input of fresh minerals-is related to the ratio of long-term average of potential evapotranspiration over average rainfall (i.e., the dryness index), and not to precipitation or streamflow only. Analysis of chemical depletion data then shows that a strong relation exists between weathering and wetness at the global scale, whereas other factors may dominate local variability. Based on this relation, we estimate the spatial distribution of weathering rates. This allows us to determine the regions of the world under supply or kinetic limitations or more susceptible to wetness. This dimensionless framework can be used to study the relation between weathering rate and all other factors (e.g., surface and subsurface topography) that contribute to the large variability of weathering rates observed globally.

The governing dimensionless groups of chemical weathering
We begin by considering the balance equations for mineral and water in the regolith (see appendix A). Tectonic uplift and physical breakdown of the parent material by tectonic and topographic stresses continuously supply fresh minerals, at a rate referred to as denudation rate D, which become subject to the action of meteoric water (Anderson 2015). Then, these minerals partly dissolve at rate W, while the remaining are eroded away at rate E without being chemically altered (Riebe et al 2004, Dixon et al 2009. The weathering flux, and hence the partitioning of the denudation flux, is determined by all the factors that control the weathering reaction (Lebedeva and Brantley 2018). These include geochemical (e.g., intrinsic weathering rate), geologic (e.g., denudation rate), as well as hydrologic parameters (e.g., precipitation and potential evapotranspiration (Budyko 1974, Sposito 2017, Daly et al 2019).
The π theorem of dimensional analysis (Barenblatt and Isaakovich 1996) allows us to establish a relation between the dimensionless ratio of weathering to denudation rate W/D (also called chemical depletion fraction) and an array of dimensionless numbers, formed by grouping the parameters affecting the weathering rate according to their dimensions. Accordingly, as detailed in the appendix A, one can write the following general relationship where W max /D and PET/P are the two major groups governing W/D globally (see section 3 and appendix), W max and PET being the intrinsic reaction rate of the mineral (e.g., moles/year) and the potential evapotranspiration (e.g., mm/year), respectively. The parameter PET/P, commonly called dryness index D I (or aridity index), is a popular climatic index used to define the partitioning of precipitation into evapo-transpiration and runoff (Budyko 1974, Daly et al 2019) and a good indicator to which to relate plant productivity and microbial activity in response to water availability (Rodríguez-Iturbe and Porporato 2004). D I thus accounts not only for water availability but also for the supply of CO 2 produced by soil respiration (i.e., by vegetation and microbes). The chemical depletion fraction can also be expressed in a form independent of the type of mineral, ( ) / W D *, by dividing (W/D) by the maximum value reached for D I → 0 (appendix B). In this way, ( ) W D * reduces to a function of D I only, which can be studied through an analysis of available data. For each experimental site of the data collection, we calculated the dryness index from long-term averages of precipitation and potential evapotranspiration. Measurements of long-term precipitation were provided by all data collections, while potential evapotranspiration was retrieved from the Climate Research Unit (CRU) dataset, from which the temporal average was then calculated. For the sites considered by Chadwick et al (2003), we used the potential evapotranspiration values provided by the authors, because gridded datasets (i.e., CRU) likely do not capture the variability of PET across the large elevation gradient where these sites are located (Kohala mountain, Hawaii). Across all sites, mean annual temperature ranges from 2°C in the Santa Rosa Mountains (NV, USA) to 25°C in the Sonora Desert (Mexico), whereas mean annual precipitation ranges from 22cm/year in Southern Sierra (NV, USA) to 420 cm/year in the Rio Icacos (Luquillo forest, Puerto Rico). The potential evapotranspiration ranges from ≈1.7 mm/day in the McNabb Track (New Zealand) to ≈5.6 mm/day in the Sonora Desert (Mexico). The resulting dryness index for these sites (figure 1) goes from ≈0.02 in the McNabb Track (New Zealand) up to 12 in one of the sites in the Kohala mountain (Hawaii, USA), representing the wet and dry extremes, respectively.

Wetness as the dominant driver of chemical depletion
The dimensionless graph of figure 1 reveals the strong control exerted by climate. As the dryness index D I goes to 0 (extremely wet climate), ( ) W D * tends to 1, while as  ¥ D I (extremely dry climate), ( ) W D * approaches 0 following a slow decay. Between these two extreme climatic regimes, the points lie on a sigmoid curve representing an intermediate, transitional regime. Finding the relation between ( ) W D * and D I can be interpreted as a problem of finding a function that satisfies these limiting conditions and has a sigmoid shape that passes through the data points. A mathematical expression that satisfies these conditions is where the only empirical constant is the exponent α. As can be seen in figure 1, the data are in good agreement with the expression (2) above; the root-mean-square error (RMSE) is ≈0.19 and ≈70% of the variance is explained ( figure C2). The characteristic exponent of the dryness index is α≈2.5. This good agreement underlines the role of wetness in controlling weathering rates at the global scale, while the dependence on other dimensionless groups becomes important when focusing on individual locations or when comparing different sites with similar D I (figure 1).  (2) is approximately 0.03, root-mean-squared error is 0.19, while the fraction of variance explained is 0.7.

Equation
(2) also provides useful insight into how wetness affects the chemical depletion fraction. In particular, it emphasizes the strongly nonlinear relation between water availability and weathering rate. Starting from high D I values, ( ) / * W D increases only slightly up to D I ≈2, after which a steep increase begins before plateauing at D I <0.5. The dryness index of approximately 2, from which the strongly nonlinear enhancement of weathering initiates, corresponds to the establishment of grasslands, savannas, and shrublands ecosystems (Woodward et al 2004), emphasizing the important role of vegetation in acidifying the soil and in turn enhancing weathering rates. Interestingly, a similar nonlinear dependence on wetness was observed in the global distribution of organic carbon sorbed to reactive soil minerals (Kramer and Chadwick 2018).

Global patterns of weathering
The relationship between chemical depletion fraction and dryness index can be used to estimate the spatial distribution of chemical depletion fraction over the globe from long-term measurements of potential evapotranspiration and precipitation. We considered the CRU dataset with over a century of climatological data  at 0.5°by 0.5°spatial resolution, from which we retrieved potential evapotranspiration and precipitation, and computed a time average for each grid point (see appendix C). From the global pattern of D I (figure 2(a)), as expected the tropical regions, e.g., the Amazon, Central Africa, and Indonesia, have the lowest dryness index (D I approaching 0), mostly due to their high precipitation regime. High latitude regions also have very low dryness index values because of their low potential evapotranspiration. Both subtropical and temperate regions exhibit high levels of dryness, including Northern Africa and Western and Central Asia, as well as west USA, Chile and Southern Argentina, and most of Australia. The relation between D I and latitude is illustrated in the longitudinal average D I , clearly showing that the subtropics are the driest latitudes.
The resulting distribution of chemical depletion fraction follows an opposite pattern compared to D I , with ( ) W D * values that increase with wetness (i.e., lower D I ), but due to the strongly nonlinear relation between ( ) W D * and D I , equation (2), the distinction between regions of slightly and highly depleted soils is more pronounced. Many regions appear to have low values of ( ) W D *, while very high values are estimated only in  3(a)). This particular frequency distribution of ( ) W D * originates, by conservation of probability, from the frequency distribution of D I , according to which low values are more frequent, and from the nonlinear relation with ( ) W D * and the power law decay with D I . Overall, we found that the global average is ( ) = W D 0.5 * and that 50% of the land has ( ) < W D 0.5 * . Having an explicit relation between chemical depletion and climate, equation (2), it was possible to define uniquely the weathering regimes by identifying the thresholds between them, for instance by using the points of the curve with maximum curvature, which are located at D I approximately 0.23 and 0.96. For D I <0.23, conditions are favorable to the weathering reaction, so that ( ) » W D 1 * and limitations may arise from the supply of fresh minerals. On the contrary, for D I >0.96, ( ) < W D 0.6 * because conditions are not favorable to the kinetics of the weathering reaction. Weathering regimes could also be defined based on a fixed absolute value of the slope of (2), and essentially the same results would be obtained by fixing a value of the slope of 0.4.

Susceptibility to wetness
To identify those areas that potentially can be more impacted by shifts in wetness, we defined the susceptibility to wetness as the negative derivative of equation (2), * which represents a measure of the sensitivity of different climatic regions and their weathering rate to variability in wetness. It is important to bear in mind that equation (3) does not provide information about timescales, for example of a transition towards a different value of ( ) / W D * upon a change in D I . A change in wetness is likely to impact directly the weathering rate, but the extent to which the latter will have an effect on the chemical depletion fraction of a soil profile will also depend on geologic (e.g., mineralogy, denudation, erosion) and surface and subsurface topographic factors, as well as temperature (e.g., (Rempe andDietrich 2014, Riebe et al 2017 , Holbrook et al 2019, Richter et al 2020)). Therefore, changes in long term chemical depletion of soil profiles should be analyzed also with respect to these other factors. As can be seen in figure 4, the susceptibility tends to zero for low and high values of D I , while it peaks at intermediate values of D I (0.5<D I <1), showing that the transitional regime also corresponds to the highly susceptible one. For regions with low or high D I , corresponding to supply and kinetic limited areas, respectively, the susceptibility to changes in wetness is very low. In humid regions (low D I ) there is more water and acidity than needed to reach the maximum weathering rates, while in dry regions such as desert and semi-desert regions (high D I ) water availability does not impact considerably the production of soil CO 2 given the poor presence of vegetation and low soil organic matter. Areas in the transitional regime, with intermediate D I , are highly susceptible and could potentially enter the supply-or kinetic-limited regimes depending on whether D I decreased or increased, respectively.
According to the thresholds in dryness index, 0.23 and 0.96, we divided the land into regions of supplylimited, kinetic-limited or transitional regime ( figure 5). Globally, about 61% of land surface is in kinetic limited conditions, including a large portion of north America and Asia, only 1% is in supply limited conditions, and the remaining 38% is in the transitional regime and thus susceptible to changes in wetness (see inset in figure 5). Susceptible regions include the African rainforest, eastern north America, central and northern Europe, as well as continental south-east Asia.

Conclusions
In conclusion, the relation between weathering, denudation, and wetness (equation (2)) introduced in this analysis sheds light on the mechanisms through which climate, specifically water availability, governs chemical  weathering at the global scale (e.g., see figure 2) and underlines the important role of vegetation and soil CO 2 production. This, in turn, may help us improve the representation of chemical weathering in landscape evolution and geologic carbon models, to better forecast the evolution of the Critical Zone as well as of the longterm climate (Kump et al 2000, Minasny and McBratney 2001, Saltzman 2001, Engler et al 2017. Knowledge of the spatial distribution of chemical weathering at the global scale is extremely important for sustainable land management (Field et al 2015, Lambin and Meyfroidt 2011), for example to identify regions that are more suitable to mineral fertilization techniques aimed at CO 2 sequestration (e.g., Enhanced Weathering (Hartmann et al 2013)) or those that are more sensitive to shifts in wetness. In the latter regions, it would be beneficial to take into account the effects of land use and land cover on chemical weathering and the potential implications for the long-term soil sustainability (Richter and Markewitz 2001).
The process of chemical weathering along a soil profile can be described by the balance equations for mineral (Riebe et al 2004) and water (Budyko 1974, Sposito 2017, Daly et al 2019, where D is the denudation rate, E is the erosion rate, W is the weathering rate, P is the long-term average precipitation, ET is the long-term evapotranspiration, and R is the long-term runoff. For the long timescales under consideration here, a common assumption to solve the water balance is that the hydrologic partitioning is mostly driven by climatic factors (Budyko 1974), such as precipitation and potential evapotranspiration, PET, so that the long term evapotranspiration can be expressed as a function of P and PET, i.e., PET being the long-term average potential evapotranspiration. The weathering flux is a function of the denudation rate, which can be considered the input of fresh minerals to the weathering zone, the mineral reaction rate constant, temperature, and the hydrologic regime, which determines the contact time between water and minerals as well as the transport of the chemical species (e.g., (Dupré et al 2003, Dixon et al 2009 and references therein). The hydrologic regime is also an important driver of vegetation dynamics and microbial activity in the upper soil layers (Rodríguez-Iturbe and Porporato 2004), and hence controls soil respiration and the CO 2 enrichment of water as it flows through the soil column towards the weathering zone. To couple the weathering rate to the water balance, we can express the weathering flux as max where W max is the maximum weathering rate for the type of mineral. Equation (A.3) is also a function of other factors, such as temperature, activation energy and so on, but as it will be shown later, at the global scale these exert only a second order control on W, so we will not treat these factors explicitly. We can select two dimensionally independent variables in equation (A.3) and reduce it to a relation between fewer dimensionless variables. Selecting the denudation rate, D, and precipitation, P, equation (A.3) can be rewritten as where f is the scaling relationship between the ratio of weathering and denudation flux, W/D, referred to as the chemical depletion fraction (Riebe et al 2004), the 'potential' chemical depletion fraction, W D max , corresponding to the chemical depletion fraction obtained if the reaction proceeded at the maximum rate, and the ratio of potential evapotranspiration and precipitation, PET/P, referred to as the dryness index D I (Budyko 1974).
As for the expression of a chemical reaction rate, consisting of the product of the intrinsic rate constant and another term which is function of the saturation, here we expect the maximum reaction rate to be a multiplicative factor for the weathering rate. Accordingly, in this dimensionless framework equation (A.4) can be written as where, because W/D varies between 0 and 1, both functions f 1 and f 2 also vary between 0 and 1. The chemical depletion fraction W/D is equal to f 1 when  D 0 I , according to which  f 1 2 . For this limit, W/D in fact reduces to a function only of W D max (figure C1). For W max /D<1, denudation rates are high compared to the maximum weathering rate W max , so we have W=W max and W/D=W max /D. For denudation rate lower than the maximum weathering rate, W max /D>1, the input of minerals from regolith production could not fully support weathering, so that W=D and hence W/D=1.

Appendix B. Definition of specific chemical depletion fraction
To focus on the relationship between chemical weathering and climate, it is convenient to focus on the function f 2 . This is done simply by dividing W/D by f 1 , so as to obtain To study the relation between chemical depletion and D I from data, the chemical depletion fraction measurements for given element were divided by the value of the maximum chemical depletion fraction, corresponding to the average value of W/D for the sites whose  D 0 I (see figure C1). Specifically, for each element we considered sites with D I <0.5 and averaged the corresponding W/D values (for Na, we removed one outlier; W/D=0.25 at D I = 0.18). Consistently with observed mobility (White and Buss 2014), we obtained » 1 . Furthermore, f 1 was greater for basalt, reflecting the higher weatherability of basaltic lithologies (Dupré et al 2003). By dividing the W/D values by the corresponding f 1 , we obtained the specific chemical depletion fraction, which does not depend on the type of mineral. The explicit functional form of f 2 is shown in equation (2) in the main text.
As we mentioned above, the analysis explicitly accounts for W max /D and D I , as these are primary factors affecting W/D at the global scale. Arrhenius plots for D I -corrected ( ) W D * indeed do not show any strong relation with T (figure C3 available online at stacks.iop.org/ERC/2/085005/mmedia). Temperature however still indirectly affects water availability through PET. While, we do not have enough information on the field sites to explore the dependence of ( ) W D * on all possible factors, our hypothesis of the primary control of D I is confirmed by its large explanatory power, i.e., D I explains 70% of the variance of ( ) W D *.
Appendix C. Global distribution of D I and ( ) W D * To compute the dryness index for the globe, we retrieved precipitation and potential evapotranspiration data from the CRU dataset (Harris et al 2014), containing over a century of climatological data (1901-2018) at a spatial resolution of 0.5°by 0.5°. For each grid point, we computed the temporal average for both precipitation and potential evapotranspiration over the period 1901-2018. The dryness index was then computed for each grid point by dividing the average potential evapotranspiration by the corresponding average precipitation. For given latitude, the longitudinal average of D I was calculated as