Articles | Volume 17, issue 9
Development and technical paper
15 May 2024
Development and technical paper |  | 15 May 2024

A revised model of global silicate weathering considering the influence of vegetation cover on erosion rate

Haoyue Zuo, Yonggang Liu, Gaojun Li, Zhifang Xu, Liang Zhao, Zhengtang Guo, and Yongyun Hu

Silicate weathering, which is of great importance in regulating the global carbon cycle, has been found to be affected by complicated factors, including climate, tectonics and vegetation. However, the exact transfer function between these factors and the silicate weathering rate is still unclear, leading to large model–data discrepancies in the CO2 consumption associated with silicate weathering. Here we propose a simple parameterization for the influence of vegetation cover on erosion rate to improve the model–data comparison based on a state-of-the-art silicate weathering model. We found out that the current weathering model tends to overestimate the silicate weathering fluxes in the tropical region, which can hardly be explained by either the uncertainties in climate and geomorphological conditions or the optimization of model parameters. We show that such an overestimation of the tropical weathering rate can be rectified significantly by parameterizing the shielding effect of vegetation cover on soil erosion using the leaf area index (LAI), the high values of which are coincident with the distribution of leached soils. We propose that the heavy vegetation in the tropical region likely slows down the erosion rate, much more so than thought before, by reducing extreme streamflow in response to precipitation. The silicate weathering model thus revised gives a smaller global weathering flux which is arguably more consistent with the observed value and the recently reconstructed global outgassing, both of which are subject to uncertainties. The model is also easily applicable to the deep-time Earth to investigate the influence of land plants on the global biogeochemical cycle.

1 Introduction

On geological timescales, Earth's climate is primarily controlled by the atmospheric CO2 concentration (pCO2); the evolution of the Sun – its brightness increases with time – also plays an important role on the timescale of a hundred million years (100 Myr) but in a temporally smooth way (Li et al., 2023). However, how the sources and sinks of CO2 varied in Earth's history remains elusive (Y. Zhang et al., 2022; Mills et al., 2021), and large uncertainties exist even in the estimate of their present-day magnitudes (Hilton and West, 2020). Due to the small size of the ocean–atmosphere carbon reservoir (∼40 000 Pg; Lee et al., 2019; Berner, 2004; Canadell et al., 2023), a small imbalance between the carbon sources and sinks can lead to large variations in pCO2 in a relatively short time (Berner and Kothavala, 2001; Berner, 1991; Walker et al., 1981; Berner, 2004). Therefore, accurately determining the exact magnitude of carbon sources and sinks is crucial for comprehensively understanding the mechanisms behind Earth's climate variations.

One of the essential ways of determining the carbon sink is through numerical modeling, especially for that in the deep past. Numerical models not only provide the magnitude of the carbon sink but also allow us to study its sensitivity to various factors such as continental evolution and climate change. Our goal here is to improve the model calculation of the primary sink of CO2, i.e., the silicate weathering, with a focus on its present-day values for which the spatial distribution is relatively well constructed.

The rate of silicate weathering is affected by the composition and physical erosion of surface rocks, pCO2, surface temperature and terrestrial runoff (Gaillardet et al., 1999; Raymo and Ruddiman, 1992; Brantley et al., 2008; Maher, 2010; Maher and Chamberlain, 2014; Dessert et al., 2003; Ibarra et al., 2019; West et al., 2005). Seawater isotopes such as Sr, Os, Li or Be are often used to estimate the global silicate weathering flux in the past (Caves Rugenstein et al., 2019; Dellinger et al., 2015; Kalderon-Asael et al., 2021; Li et al., 2019). However, it is difficult to constrain the sensitivity of silicate weathering to certain factors (e.g., temperature) from such measurements, especially in local regions, due to both the uncertainties in their interpretation (Li et al., 2019; Dellinger et al., 2015) and their global nature. Simulating the weathering reactions in the lab can provide useful information for the factors that control the weathering rate, but lab conditions are generally much simpler than those in the natural field (Gruber et al., 2014; Calabrese et al., 2022; White and Brantley, 2003). Many other works focused on compiling the dissolved river loading to estimate the silicate weathering fluxes and rates in different regions for the present day (Bluth and Kump, 1994; Gibbs et al., 1999; Amiotte Suchet et al., 2003; Suchet and Probst, 2002). Despite the various uncertainties in these methods, they provide a basis for the development of numerical models.

Early zero-dimensional models (e.g., the Geologic Carbon Cycle (GEOCARB) family; Walker et al., 1981; Berner et al., 1983; Berner, 1991) and subsequent two-dimensional numerical models such as the Gibbs and Kump Weathering Model (GKWM) in 1994 (Bluth and Kump, 1994), the Global Erosion Model for CO2 fluxes (GEM-CO2) in 1995 (Suchet and Probst, 2002; Amiotte Suchet et al., 2003) and a model by Jens Hartmann in 2009 (Hartmann et al., 2009; Hartmann and Moosdorf, 2012; Hartmann et al., 2014) provided important understanding of the long-term carbon cycle. Studies using these models (Amiotte Suchet et al., 2003; Gibbs et al., 1999; Zhang et al., 2021) have identified the lithology and runoff as the strongest predictors of chemical weathering rates. However, basin- or catchment-scale compilation of weathering data (Gaillardet et al., 1999) indicates that the spatial variability of the weathering rate had to be explained through a combined effect of runoff, temperature and erosion rate. West et al. (2005) further showed that there were two-end-member schemes of the weathering – transport-limited and kinetically limited regimes.

Built on the work of West et al. (2005), Gabet and Mudd (2009) constructed a theoretical model (referred to as the GM09 model hereafter) that encompassed the continuum of these two weathering regimes for the first time. This model is probably the most sophisticated one to date in terms of global silicate weathering calculation and has been used in many works subsequently for both the present day and the past (West, 2012; Goddéris et al., 2017; Maffre et al., 2018; Park et al., 2020). However, the model contains a few unknown parameters, including cation abundance in the bedrock, the dissolution rate constant and its dependence on runoff and reaction time, and the regolith production rate, of which only rough ranges are given. Most previous works (Maffre et al., 2018; Park et al., 2020) using this model estimated these parameters through some fitting approach with the help of catchment-scale observations (Gaillardet et al., 1999).

The global total silicate weathering flux (Fw) of the present day given by Park et al. (2020) (referred to as Park20 hereafter) in terms of carbon is 4.5×1012 mol yr−1, which was thought to be consistent with the global outgassing rate estimated by Gerlach (2011). However, a few lines of evidence indicate that this flux may be overestimated. (1) The Fw estimated from the present-day observations is 2.5×1012 mol yr−1 (1.59×10122.75×1012 mol yr−1) (Gaillardet et al., 1999; Moon et al., 2014). (2) The global outgassing rate was re-estimated to be ∼23.3×1012 mol yr−1 by Müller et al. (2022). (3) The silicate weathering fluxes for individual river basins within the tropical region from the Park20 model were overall overestimated compared to the observations (Fig. 1b), which led to an overestimate of Fw (Fig. 1c). This overestimation over the tropical region by the Park20 model has also been argued to exist based on the observed 187Os/188Os (Caves Rugenstein et al., 2021).

Overestimation of the carbon sink by 100 % will lead to a dramatic decrease in pCO2 and an extreme icehouse climate in a few million years when the outgassing is fixed (Berner and Caldeira, 1997; D'Antonio et al., 2019) and thus should be dealt with properly. More important reasons may be that (1) the overestimation is not random among different sites but is systematic; the weathering fluxes over tropical river basins are much more likely overestimated than underestimated, whether in the original values (Fig. 1b) or in the logarithmic values (Fig. 1e). (2) The climate sensitivity of the silicate weathering, i.e., the ability of silicate weathering to stabilize the climate, may be mis-estimated due to this systematic error.

The lower-than-expected silicate weathering rate over the tropical region was noticed by Stallard as early as 1981 (Stallard and Edmond, 1981; Stallard, 1985; Stallard and Edmond, 1983). Godderis et al. (2008) and Hartmann et al. (2014) also found that considering only the effects of temperature and runoff would lead to a significant overestimation of weathering in the tropical region. They proposed the effect of soil shielding as a solution; i.e., the occurrence of leached soil in equatorial regions hinders deeper weathering. They then assumed a global soil shielding effect in regions with leached soil and improved their model performance. However, the soil shielding effect has already been considered in the GM09 model to some extent, where the physical erosion was parameterized. Therefore, the problem remains in this model, and our main goal in this paper is to find a simple solution to the problem and test whether it affects the sensitivity of global silicate weathering to climate change.

Figure 1The difference between the model calculated and observed silicate weathering fluxes for 81 large rivers (more details can be found in Sect. 2.2.5) over the world. The upper and lower panels show model-obs and log10(model)-log10(obs), respectively. The left, middle and right panels show rivers at the mid to high latitudes (if more than half of its river basin is located at or beyond 30° latitude), low latitudes (within 30° latitude) and over the whole globe, respectively. The model results were calculated using the GM09 model but with model parameters in Park20. The global total weathering flux is 4.5×1012 mol yr−1. The surface slope and all climate forcings are from Park20, in which the runoff used is the one denoted as “from Yves” in Park20 (more details in Sect. 2.2.1). A similar systematic upward bias in the tropical region appeared when the parameters as given in Maffre et al. (2022) were used (Fig. S1 in the Supplement).


Specifically, we will first test whether the historical climate data constructed by different institutes have any significant impact on the calculated silicate weathering rate using the GM09 model in the tropical region. Then, the influence of the magnitude of the seasonal cycle on physical weathering is tested next. It is then found that reducing physical erosion rates where leached soil is present works best in removing the systematic bias in the tropical region. In the end, we find a simple parameterization scheme related to vegetation that can attain a similar effect to that of leached soil but that is much more applicable to weathering calculation for other periods of Earth's history.

The rest of the paper is organized as follows. In Sect. 2, the GM09 model is briefly described, and the field observations used to validate the model and climate data used to calculate the weathering fluxes are also described. In Sect. 3, the results of various sensitivity tests and the parameterization for the vegetation effect are presented. The shortcomings as well as the consequences of the model revision are then discussed in Sect. 4, and a summary is provided in Sect. 5.

Figure 2Schematic diagram of the theoretical model of bedrock weathering and the simultaneous production of soil or regolith based on GM09. In stage 1, the unweathered bedrock moves vertically at a speed U due to tectonic movement, with weathering and erosion just to occur at the surface. In stage 2, soil is produced (Pr) at the surface of the bedrock and eroded (E) at the soil top, with silicate weathering occurring mostly within the soil. h represents the soil thickness, and h¯ and h̃ are the heights of the soil and the bedrock surface, respectively, relative to the reference plane (note that h¯ is nonzero in stage 2 but is not marked for the sake of esthetics). The part enclosed by dashed lines is eroded away. All variables evolve with time at this stage. In stage 3, a steady state is reached under continuous weathering such that the soil thickness and the weathering flux do not change with time anymore. The weathered material within the soil is carried away by water runoff into the oceans, with the weathering flux denoted as W.


2 Model and data

2.1 Theoretical model for silicate weathering

2.1.1 The weathering profile and weathering flux

For the convenience of the latter discussion, the model GM09 as presented in detail in Park20 and Maffre et al. (2022) is recapped here. The model includes an explicit simulation of a regolith layer, which extends from the soil surface to the unweathered bedrock (Fig. 2). The layer can be millimeters to tens of meters thick depending on the environment and is determined by

(1) d h d t = d h ̃ d t + d h ¯ d t = P r - U + U - E = P r - E ,

where h is the regolith thickness, Pr is the soil production rate and E is the erosion rate. The weathering rate J at depth z is proportional to the concentration of cations (e.g., Ca2+ and Mg2+) denoted as x and also depends on the temperature (T), runoff (q) and exposure time (τ) that the sample has experienced. The influences of T and q are generally considered using the Arrhenius equation and a linear or power-law relation (White and Blum, 1995; Dessert et al., 2003), respectively. When an exponential dependence of the weathering rate on runoff q is employed as in Park20, the weathering rate J is written as

(2) J ( z ) = K ( 1 - e - k w q ) e - E a R 1 T - 1 T 0 τ σ x ( z ) ,

where K is the dissolution constant, kw is the runoff sensitivity of the dissolution rate, Ea is the apparent activation energy at T0 for dissolution, R is the gas constant, and σ is an empirical constant.

It should be noted that the global total silicate weathering flux throughout this work pertains specifically to the weathering flux of Ca and Mg silicates. While Na and K silicates also participate in weathering, these are not traditionally regarded as carbon sinks on geological timescales (Berner et al., 1983) due to their inability to form carbonate minerals. However, the residence times of Na+ and K+ in the ocean are ∼80 and ∼10 Myr (Lécuyer, 2016; Emerson and Hedges, 2008; Olson et al., 2022; Berner and Berner, 2012; Hu et al., 2020), respectively. This long residence time means that the weathering of Na and K silicates could have an impact on the atmospheric CO2 on a 1 million year timescale. Moreover, Na+ and K+, when released into the soil through weathering reactions, may displace Ca2+ and Mg2+ through cation exchange with sediments or oceanic crust (France-Lanord and Derry, 1997), leading to carbonate deposition and carbon sinking indirectly. However, currently we are unable to quantify these aspects due to the intricacies of the Na and K cycles. Thus, we focus solely on the Ca2+ and Mg2+ silicate weathering flux in the current study.

The concentration of cations themselves changes with time according to

(3) x t = - U x z - K ( 1 - e - k w q ) e - E a R ( 1 T - 1 T 0 ) τ σ x .

In most cases, we do not need to track the evolution of surface topography, and it is as accurate, when calculating weathering flux, as just setting the reference plane to be at the regolith–bedrock interface. In that case, h̃0 and h¯h, and the uplifting speed in Eq. (3) can be replaced with Pr. The weathering profiles are often assumed to have reached a steady state; i.e., the soil production rate equals the erosion rate (Phillips, 2010). This assumption is appropriate if the lifetime of a weathering profile is much shorter than a few million years. The lifetime of a weathering profile may be estimated by using its typical thickness and the surface erosion rate. The global total erosion is ∼20 Gt yr−1 (Milliman and Farnsworth, 2011), which gives a global mean erosion rate of ∼133 t km−2 yr−1. If we use a relatively low value, say 60 t km−2 yr−1 (equivalent to 2×10-5 m yr−1), a typical weathering profile 10 m thick will require half a million years to completely renew. Thus, a weathering profile is near steady state if the environment changes slowly over a few million years. Such an assumption is not ideal but is necessary to make in order to study the long-term (hundreds of millions of years) evolution of silicate weathering at a reasonable cost.

Under such an assumption, neither the soil production rate Pr nor the cation concentration x changes with time (i.e., xt=0) as long as the tectonic setting and climate have not changed, and the exposure time τ is simply z/Pr. Equation (3) then becomes

(4) - P r x z - K ( 1 - e - k w q ) e - E a R 1 T - 1 T 0 z P r σ x = 0 .

The total weathering flux at the grid point is just the integration of J(z) through the regolith:

(5) W = 0 h J ( z ) d z = 0 h K ( 1 - e - k w q ) e - E a R 1 T - 1 T 0 z P r σ x d z .

There are still two undetermined variables in the formula above, i.e., h and Pr. The regolith thickness h can be calculated by assuming the balance between the soil production rate Pr and the surface erosion rate E. Next, we will describe how Pr and E are parameterized.

2.1.2 Soil production rate

Studies showed that the soil production rate could be controlled by temperature or water content (Heimsath et al., 1997, 2009; Dixon et al., 2009; Whipple et al., 2012; Carretier et al., 2014). Overall, the soil production rate has been found to decline exponentially with increasing depth of the regolith (h) due to the decrease in water percolation or biogenic disturbance (Dietrich et al., 1995; Heimsath et al., 1997, 1999; Riebe et al., 2004; Heimsath et al., 2009; Heimsath and Korup, 2012; Burke et al., 2007; Small et al., 1999). However, it has also been suggested that there is an optimum regolith thickness. Soil production also slows down when the regolith is too thin in certain environments (Anderson, 2002; Strudley et al., 2006). The soil production rate has thus been described by the so-called “humped” law:

(6) P r = k rp q e - E a R 1 T - 1 T 0 ( e - h d 0 - k 1 e - h d 1 ) ,

where the second exponential term in the parentheses is there to ensure that the soil production rate decreases when h is too small. Here we neglect this effect by setting k1 to 0, the same as in Park20. In Eq. (6), krp is the regolith production constant to be determined by fitting the observations and d0 is the attenuation depth and is set to 2.73 m, the same as those in Park20.

2.1.3 Erosion rate

The current estimation of the erosion rate is mainly from the suspended river loads (Milliman and Farnsworth, 2011) or in situ cosmogenic nuclides in river sediments (Wittmann et al., 2011, 2015, 2020; Blanckenburg et al., 2012; Dannhaus et al., 2017; Larsen et al., 2014). Supported by observations, modeling studies of erosion rates at a global scale have flourished, and several parameterization schemes are now available. For example, the model BQART, derived from a global database of 488 rivers, can estimate the erosion flux for the entire river basin with knowledge of water discharge, drainage area, basin relief, average temperature and anthropogenic influence (Syvitski and Milliman, 2007).

The river incision at the catchment scale is simulated using the classical empirical law – the stream power incision law (Davy and Crave, 2000; Howard, 1994) that has been widely used (Adams et al., 2020; Gasparini et al., 2007; Harel et al., 2016; Lague, 2014; Quye-Sawyer et al., 2020; Royden and Taylor Perron, 2013):

(7) E = k e B q m s n ,

where ke is the erodibility constant which is calibrated by setting the global total physical denudation flux to 20 Gt yr−1 and set to 0.0030713 m1−m/yr1−m in Park20. s is the surface slope. Exponents m and n are set to the values 0.5 and 1, respectively. A new parameter B is introduced herein to match the observed individual erosional fluxes in some of the tests performed herein, as will be explained in detail in Sect. 2.2.5. The BQART model is similar to Eq. (7), except that a temperature dependence is added (Syvitski and Milliman, 2007). This model was tested here, but results will not be shown because no improvement was achieved compared to the stream law model above.

Note that both the BQART and stream law models are not prepared for the grid-scale erosion rate but for the catchment scale. More explicit ways of representing the denudation are available (e.g., Carretier et al., 2018), which involve many detailed processes and hydrographic features. Such a method is not practical here since our purpose is to construct a model applicable to paleoclimate conditions for which limited information can be obtained.

2.1.4 The final solution for the weathering flux

The regolith thickness h in Eq. (5) can be calculated by equating the erosion rate E and soil production rate Pr under the steady-state assumption:

(8) h = P r - 1 ( E ) .

Since h, Pr and E are independent of z, the integration in Eq. (5) can be solved to get

(9) W = E x | z = 0 - x | z = 0 e - K 1 - e - k w q e - E a R 1 T - 1 T 0 σ + 1 h E σ + 1 ,

where x|z=0 is the concentration of relevant cations in the fresh rock and is dependent on the lithology. The second term in the large parentheses of Eq. (9) is actually the concentration of elements at the surface of the regolith layer (i.e., z=h) and will be represented by xs in what follows.

Five parameters (Table S1 in the Supplement) in this equation are unknown. Field measurements or laboratory experiments have provided reference ranges for some parameters (Rudnick and Gao, 2003; Heimsath et al., 1997; White and Brantley, 2003). Based on these reference ranges, previous studies estimated optimal values of these parameters by fitting the calculated weathering fluxes with the observed ones in various river catchments (Maffre et al., 2018, 2022; Park et al., 2020). We will use this theoretical model as a foundation and try to improve the model–data comparison by adding possible missing processes. The parameters in Eq. (9) are re-estimated when necessary.

2.2 Data

2.2.1 Climate data for the present day

The climate fields required in the model presented above are surface temperature and river runoff. To investigate the influence of these data on the comparison between the calculated and observed weathering fluxes, climate data from various sources are considered. The first one is the monthly 2 m temperature and runoff for 1950 to 2021 obtained from ERA5 (Muñoz Sabater, 2019). ERA5 is a re-analysis dataset obtained using a global climate model constrained by various observations from weather stations, ships and satellites. The dataset is gridded with a spatial resolution of 0.1° × 0.1°. Since Park20 has done elaborate work on testing parameters, we also used the temperature and runoff in their test. Their temperature was derived from CRU TS v.4.03 (Harris et al., 2014; denoted as T_CRU), while two runoff datasets were used: one was from UNH/GRDC Composite Runoff Fields V1.0 (Fekete et al., 2002; denoted as R_Park), and the other was from Yves as described in the runoff file provided by the data repository supplied along with Park20 (denoted as R_Yves). However, because the R_Park data are different from the runoff that we downloaded from UNH/GRDC Composite Runoff Fields V1.0 (, last access: 7 May 2022), this latter dataset was also tested and denoted as R_UNH herein. Other than these two datasets, an observation-based global gridded runoff dataset GRUN from 1902 to 2014 (Ghiggi et al., 2019) with a resolution of 0.5° × 0.5° was also used.

To account for the influence of global warming and human activities, we conducted tests using temperature and runoff averaged over three different periods. For temperature, the three time periods are 1950–1979, 1950–1997 and 1950–2021, denoted as T_ERA1, T_ERA2 and T_ERA3, respectively. For runoff, the three time periods are the same as those for the temperature for the ERA dataset but are 1902–1950, 1902–1996 and 1902–2014 for the GRUN dataset and denoted as R_GRUN1, R_GRUN2 and R_GRUN3, respectively. The distributions of temperature and runoff in different datasets and different time periods are shown in Figs. S2 and S3.

2.2.2 Climate data for the Last Glacial Maximum (LGM) and the future

To estimate the sensitivity of global silicate weathering (i.e., Fw) to the climate, data for both cold and warm climates are needed. For cold climates, the LGM was chosen and the data from M. Zhang et al. (2022) were used, denoted as T_LGM and R_LGM. For the warm climate, the abrupt quadruple-CO2 experiment carried out using CESM2 (Danabasoglu, 2019) was used, and data were downloaded from the CMIP6 data website (, last access: 1 July 2023), denoted as T_4CO2 and R_4CO2, respectively.

2.2.3 Surface topography

A key variable for calculating the erosion rate is the surface slope s. Global topography data from Scotese and Wright (2018) were used to calculate s according to the formula (Maffre et al., 2018)

(10) s = h x 2 + h y 2 .

The slope data from Park20 were also tested, whose topography field was from the Shuttle Radar Topography Mission (Farr et al., 2007). We denote the surface slopes calculated from Scotese and Wright (2018) and from Park20 as s1 and s2, respectively (Fig. S4).

2.2.4 Lithology

The spatial distribution of lithologies was obtained from the Global Lithologic Map (GliM) (Hartmann and Moosdorf, 2012). The original dataset includes 16 types of rock, and we grouped them into six categories, the same as done in Park20 (see their Fig. S1 and our Fig. S5). The concentrations of Ca and Mg cations in each type of rock can be estimated through the EarthChem library (, last access: 17 October 2022), in addition to rocks such as sedimentary and metamorphic rocks, whose characteristics are greatly dependent on protoliths. They may cause large uncertainty in the calculated silicate weathering flux, so Park20 treated the concentrations of these two types of rocks as fitting parameters in the model. This is also how it is done here.

2.2.5 Catchment measurements of weathering and erosional fluxes

For model validation, concentrations of cations such as Ca2+ and Mg2+ in the dissolved loading of river discharge from global catchments were collected from the literature. The weathering fluxes integrated over the corresponding river basins can be inferred from these catchment data. Cations in rivers have various origins, such as atmospheric input, carbonate weathering or silicate weathering (Moon et al., 2014). Since almost only Ca2+ and Mg2+ from silicate weathering can be considered a sink of atmospheric CO2 on a geological timescale, the elements from different sources have to be distinguished. Two standard methods have been widely used to differentiate silicate and non-silicate chemical sources. The forward method often uses the pre-assigned compositions for each element, which essentially relies on the knowledge of bedrock and environmental characteristics of the study area (Meybeck, 1987; Edmond et al., 1995; Galy and France-Lanord, 1999). In general, this approach is more easily applicable to small watersheds or watersheds with monolithic lithology than to large and complex watersheds. The inverse method starts from a priori ranges of elemental concentration ratios and determines the best a posteriori ratios based on the mass balance equation. This approach is useful when complete information on chemical compositions within the watershed is not available, such as in some large catchments (Gaillardet et al., 1999; Moon et al., 2014).

Since the silicate weathering model is mostly used for the geological past, where detailed information on surface topography, climate and lithology is not available, the spatial resolution of the model cannot be too high: usually around 0.5° × 0.5° or coarser. To ensure a comparable performance of the model for the past to the present day, the spatial resolution used herein is 0.5° × 0.5°. At such a coarse resolution, accurate identification of river routes is not possible, and data compiled for relatively large river basins are more reliable for model validation. Two such datasets are available (Gaillardet et al., 1999; Moon et al., 2014), and that compiled by Gaillardet for the 51 large river basins is the focus of our analysis (see Table S2 for the values and Fig. S6a for the definitions of the basins). Park20 removed the Brahmaputra watershed because it overlaps with the Ganges watershed. They also removed the Don watershed in their parameter exploration. Here we employed the modern river direction files contained within the Community Earth System Model (CESM) to refine the geographical delineation of rivers, ensuring that the Brahmaputra watershed was distinguished from the Ganges watershed. We also kept the Don watershed. As will be shown later, including the Brahmaputra and Don watersheds has little effect on the results.

Park20 also incorporated data from HYBAM, which consists of 32 small watersheds in the Amazon region (Moquet et al., 2011, 2016, 2018). The average weathering flux from the HYBAM Amazon basin data is approximately 0.07 mol m−2 yr−1, while the average weathering flux from Gaillardet et al. (1999) for the Amazon region is 0.02 mol m−2 yr−1. Due to this significant mismatch between the datasets, we used both the Gaillardet et al. (1999) data (denoted as Gaillardet) and the combination of Gaillardet and HYBAM data (denoted as Gaillardet+HYBAM) to validate the model.

The modeled erosion rates can also be validated to some extent by the observed suspended river loading, the so-called total suspended sediment (TSS). Different from the dissolved cations in the water, a significant portion of the suspended loading may have been deposited before they reached the catchment. Therefore, the suspended loading measured at the catchment may not represent the erosion rate over the river basin well. Nevertheless, we collected the river loading measurements from four sources (Table S3; Milliman and Syvitski, 1991; Milliman and Farnsworth, 2011; Milliman et al., 1995) and obtained the loading for each of the 51 large rivers mentioned above. Multiple measurements may be available at one river catchment; we prefer the older value in order to minimize the influence of human activities (see the details in Table S3).

Mean denudation rates are also available from cosmogenic nuclide analysis in sediment, like in situ cosmogenic 26Al and 10Be. In general, this represents a longer-term average erosion rate, typically on the scale of millions of years, unlike TSS, which represents the erosion over a short time period ( years). As a result, the denudation rates obtained through cosmogenic nuclide analysis may exclude the anthropogenic influence. Wittmann et al. (2020) compiled global denudation rates for more than 50 large rivers over a range of climatic and tectonic regimes in this way, but only 18 of the rivers overlap with our data. The final loading thus obtained is shown in Table S3. Figure 3 shows that the model-calculated erosion rates (Eq. 7 with B=1) deviate significantly from both TSS and isotope-derived erosion rates. Therefore, in some of the tests with the original Park20 model, the model-calculated erosion rates were scaled by tuning B such that the erosion of each basin was identical to the observed one. Note that, in these tests, B is a constant within each basin, but the erosion rate at each grid point varies within the basin. Moreover, if neither TSS nor the cosmogenic nuclide data are available for a river basin, B is set to 1 for this basin.

Figure 3The comparison between model calculated and observed erosional fluxes for individual river basins. The Park20 model is used with the forcing data R_Yves and s1 (defined in Sect. 2.2). Different colors represent the regions where the basins are located, and the sizes represent the areas of the basins. In panel (a), the observed erosion rates are from TSS data (the last column in Table S3), and the observed erosion rates in panel (b) are from cosmogenic nuclide analysis, from which data are available for only 18 rivers (the penultimate column in Table S3).


2.2.6 Vegetation

The primary vegetation data used herein are the areal fraction of different vegetation types and their associated leaf area index (LAI) provided by NCAR (Fig. S6c–d), which are derived by integrating observed land information (Lawrence and Chase, 2007). To test the performance of simulated vegetation, we also downloaded the preindustrial vegetation data simulated by the Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) dynamic vegetation model and the HadCM3 climate model (Allen et al., 2020). In the 4×CO2 experiment, the vegetation changed with climate and the data were downloaded from the CMIP6 home page, while the vegetation was assumed to be the same as in the present day except where the land was covered by ice sheets in the LGM experiment by M. Zhang et al. (2022).

2.2.7 Leached soil

The global soil distribution data are obtained from the Harmonized World Soil Database v1.2 (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012), which is provided by the Food and Agriculture Organization of the United Nations. Following Hartmann et al. (2014), we selected six specific soil types as leached soil, i.e., Ferralsols, Acrisols, Nitisols, Lixisols, Histosols and Gleysols. Figure S6b represents the proportion of leached soil within each grid cell, as determined according to the selected soil types.

2.3 Evaluation of model performance

The model–data discrepancy in silicate weathering flux is often measured by r2 (e.g., Park et al., 2020):

(11) r log 2 = 1 - ( log 10 M i - log 10 O i ) 2 ( log 10 O i - log 10 O ¯ ) 2 ,

where Mi and Oi are the model calculated and observed values, respectively, and the summation is over the index i. Since we are concerned with the global flux Fw and the weathering–climate sensitivity, Mi and Oi represent the catchment-weathering flux for river i rather than the weathering flux per unit area of the ith river basin. In the equation above, a logarithmic operation is taken to the values first before calculating the difference: a subscript “log” is thus added to differentiate it from the r2 calculated using the original values directly.

(12) r 2 = 1 - ( M i - O i ) 2 ( O i - O ¯ ) 2

Using rlog2 has the advantage of giving relatively balanced weights to both the very small and very large values, which is important because the weathering fluxes over different river basins differ a lot (Table S2). Park20 obtained their model parameters in Eq. (9) by maximizing rlog2. However, although there is a relatively small systematic bias in the logarithmic model–data errors (the data points distribute more symmetrically against the zero line in Fig. 1f), Fig. 1a–c show that there is an obvious systematic bias in the direct model–data errors. For similar magnitudes of the observational silicate weathering fluxes, the bias is much larger over the low-latitude regions (Fig. 1b) than over the high-latitude (Fig. 1a) regions. The bias in the direct errors in Fig. 1b will lead to an overestimation of the global weathering flux Fw (the global integral of W in Eq. 9) and may also be a misestimate of the weathering–climate sensitivity. Therefore, we argue that using the sum of rlog2 and r2 (denoted as “R2” hereafter) is better than using either of them as the criteria of model validation.

2.4 Experiments

In the first set of experiments, the original model of Park20 is tested for the influence of climate data and erosion rates from different sources or the same source but in different time periods. As described above, the temperature data come from two sources, ERA5 and CRU, and the data from ERA5 are organized into three different time periods; the runoff data come from five sources, ERA5, GRUN and UNH from Park20, UNH updated herein and Yves, where both the ERA5 and GRUN data are also organized into three different time periods. Slope data come from two sources: Scotese and Wright (2018) and Park20. The erosion rates are calculated in three different ways which all used Eq. (7), but the parameter B has different values: B=1, B tuned according to TSS data and B tuned according to both TSS data and the cosmogenic nuclide analysis. In the last case, the cosmogenic nuclide analysis supersedes TSS data if both of them are available for a basin. There are 4×9×2×3=216 experiments in total, which are summarized in Table S4.

In the second set of experiments, we try to improve the Park20 model by considering the effect of additional processes. In each of these experiments, rather than adopting the values from Park20, all the unknown model parameters (Table S1) are optimized again. Based on the results of the first set of experiments, only T_CRU is used for temperature, and R_Park, R_Yves and R_GRUN2 are used for runoff. Erosion correction is not applied (i.e., B=1), and both slope data s1 and s2 are tested. We will first show that changing the validation criteria (maximizing R2 rather than rlog2) is able to alleviate the systematic bias so that there is no overall overestimation, but the model–data discrepancy becomes even larger. In order to reduce this discrepancy, we try three different methods. The first method is to consider the influence of the seasonal cycle of temperature on the soil production rate, which will change the regolith thickness. The second and third methods consider the influence of leached soil and vegetation, respectively, on erosion rates. All three methods act to reduce the silicate weathering fluxes in the tropical region relative to those at the mid to high latitudes. All of these experiments are summarized in Table 1 below.

To consider the effect of vegetation, two different approaches have been tried, denoted by “m1” and “m2” in Table 1. m1 and m2 use the LAI of global vegetation from NCAR and simulated by the LPJ vegetation model, respectively. The global total erosional flux in m1 and m2 is reduced due to the shielding effect of vegetation. In order for the global total erosional flux to remain consistent with the observed value (20 Gt yr−1 in Park20), the erosion rate at every grid point is scaled uniformly (by changing ke in Eq. 7). For the sensitivity test, the LAI values of the global vegetation of the LGM (“m3”) and 4×CO2 (“m4”) are used (see Sect. 2.2.2).

Table 1Summary of the second set of experiments.

a m1–m4: the global LAI is used and the global erosional flux is fixed to 20 Gt yr−1 in all four cases, but the vegetation is from NCAR, the LPJ model, the LGM experiment and the 4×CO2 experiment, respectively. b The maximum R2 here is obtained with observed silicate weathering fluxes from Gaillardet and HYBAM data together.

Download Print Version | Download XLSX

3 Results

We will first show whether the overestimated weathering fluxes over the tropical river basins of the Park20 model were due to the uncertainty in climate data or error in the calculated erosion rates. Then, we will re-estimate model parameters by balancing rlog2 and r2, i.e., by maximizing R2 defined in Sect. 2.3. After that, we propose and test a few different parameterizations to see whether they are effective in further decreasing the model–data discrepancy measured by R2 (Table 1). Without specific indications, all the results described below are for the present day.

Figure 4The r2 (solid symbols) and rlog2 (hollow symbols) calculated using different temperature, runoff and slope data. In panel (a), all the observed catchment-weathering fluxes in Park20 are used, while in panel (b) only the 51 basins of Gaillardet et al. (1999) are used to calculate r2 and rlog2. The runoff datasets are denoted on the x axis, and their names can be found in Fig. S3. Circles and pentagrams denote results calculated using slope data s1 and s2, respectively. Blue and red denote results calculated using the temperature data T_CRU and T_ERA3, respectively. The results using the temperature data T_ERA1 and T_ERA2 are very similar to that using T_ERA3 and thus are not shown. Panels (c) and (d) are similar to panels (a) and (b), except that here the temperature is fixed to T_CRU, while colors mean different ways of revising the erosion rate: no change (blue), the erosion rate of each basin scaled according to TSS data (orange) or cosmogenic nuclide analysis (green).


3.1 Influence of the climate forcing and erosion rate in the original Park20 model

For this series of tests, everything is the same as the Park20 model, except that the temperature, runoff and surface slope from different sources or different time periods are used. Results show that climate and slope data do have some impact on rlog2 or r2, especially the latter (Fig. 4a, b). The runoff data have the largest impact, followed by slope, and the temperature data have the least impact, probably because the uncertainties in temperature are small (Fig. S2). For runoff, the data from different centers give quite different r2 values, while the data from the same center but different periods have a small effect. Although r2 can vary from −0.5 to −4.47 in different cases, all of them are below zero (Fig. 4a, b), meaning a large model–data discrepancy. For all cases, overestimation in the weathering fluxes over tropical river basins persists (not shown but largely the same as shown in Fig. 1b), and the total global weathering flux is similarly overestimated.

If the observed erosion rates are used, r2 is significantly improved, especially when the runoff datasets R_UNH and R_Park are used. The improvement is more significant when the erosion rates inferred from the cosmogenic nuclide analysis (Wittmann et al., 2020) are used. The tropical bias is also reduced but is still quite obvious (Figs. 5 and S7). Note that the results are improved even without tuning the empirical parameters in the Park20 model. This test shows us that the erosion rate may be a critical factor in alleviating the model bias. However, the erosion rates in either the past or the future are unknown and need to be parameterized if the model is to be applied to these time periods. Improving this parameterization is the major focus of our work herein and will be described in detail in what follows.

Figure 5The difference (model  observation) in silicate weathering fluxes for 81 large rivers (more details can be found in Sect. 2.2.5). The model results are from the experiment T_CRU_R_Yves_s1_Be (Table S4), and the observation data are Gaillardet+HYBAM. The left, middle and right panels show rivers at the mid to high latitudes (if more than half of the river basin is located at or beyond 30° latitude), low latitudes (within 30° latitude) and over the whole globe, respectively. rlog2 and r2 are 0.54 and 0.31, respectively. The global total weathering flux is 3.95×1012 mol yr−1.


3.2 Maximizing R2 – a new control model

Other than the inaccuracy in the erosion rate (Fig. 3), the systematic bias in the Park20 model (Fig. 1) may also be due to the model parameters being searched by maximizing rlog2. Here we check whether the bias can be alleviated by minimizing R2. Specifically, five parameters are searched with their searching ranges given in Table 2. Because the computational load of the model is relatively small, the searching is done by a forward calculation for all possible combinations. The total number of combinations is 240 240, and a full search takes 72 h on a desk computer and 1 h when 72 cores are used on a cluster. Only results for s2 are shown here, which aligns more closely with observations than those for s1 (not shown); results for s1 can be found in Table 1. Moreover, because of the relatively high sensitivity of the model results to runoff (Fig. 4), the model parameters are searched for three runoff datasets: R_Yves, R_Park and R_GRUN2. Only results for the former two are presented below, which is sufficient for demonstrating the effect of maximizing R2.

When calculating rlog2 and r2, two different sets of observed catchment-weathering fluxes have been used: Gailladet and Gailladet+HYBAM (see Sect. 2.2.5). The rlog2, r2 and R2 of all the parameter combinations are shown in Fig. 6, where each dot represents the result of a specific combination of model parameters and only the ones with values greater than 0 are shown. For Gailladet, the maximum rlog2 and r2 are −0.044 and 0.323, respectively, when R_Park is used (Fig. 6a; orange and light-blue dots, which represent rlog2 and R2, respectively, do not show up because all their values are smaller than 0). For Gailladet+HYBAM, the maximum rlog2 and r2 become 0.138 and 0.349, respectively (Fig. 6b). It can be seen that Fw tends to be overestimated if rlog2 is to be maximized (Fw=5.54×1012 mol yr−1 at the peak of the orange dot group in Fig. 6b) but underestimated if r2 is to be maximized (Fw=1.8×1012 mol yr−1 at the peak of the purple-red dot group in Fig. 6b). Although the maximum rlog2 (orange dots) and r2 (purple-red dots) are both greater than 0 in Fig. 6b, no parameter combination can give relatively high rlog2 and r2 simultaneously, so that light-blue dots can appear. Using R_Yves improves rlog2 significantly and thus R2; the maximum R2 value for Gailladet+HYBAM is 0.483 (Fig. 6d). It is notable that Fw is within the observational uncertainty range when R2 is maximized. The parameter combination associated with the maximum R2 is hence considered the new control model, and R_Yves is used in all the tests to be presented in what follows.

When R2 is maximized, either rlog2 or r2 or both are too small (Fig. 6). This means that errors for individual basins have increased overall, although the signs of errors are more balanced (Fig. 7) than before, so that the bias in Fw is small. However, inspection of the data points in Fig. 7 shows that the errors in the high-latitude region now have a negative bias compared to before (compare Figs. 7a and 1a), while the positive bias in the tropical region is somewhat reduced but remains (Fig. 7b). This redistribution of biases is clearly unsatisfying, and it may suggest that there is a missing process that distinguishes the tropical and extratropical regions.

Table 2Model parameters and their values to be searched.

a The data marked in italics are the additional values considered herein on top of those searched by Park20. The bold black values represent the optimal parameters selected by Park20. b Although the range of the cation concentration of metamorphic rocks overlaps with the sedimentary rocks, it is constrained so that the former must be larger than the latter during the search.

Download Print Version | Download XLSX

Figure 6The rlog2 (orange) and r2 (purple-red) together with their sums (light blue) for all possible combinations of the parameters in Table 2. Only the cases with values greater than zero are shown. The results in panels (a), (b), (c) and (d), respectively, are for experiments R_park_s2 and R_ Yves_s2 defined in Table 1 and differ only in the runoff data used. The Gailladet and Gailladet+HYBAM are used as observational data in the left and right panels, respectively. The black vertical line and grey zone show the observed global total weathering flux (i.e., Fw) and its uncertainty range.


Figure 7The difference (model  observation) in silicate weathering fluxes for 81 large rivers (more details can be found in Sect. 2.2.5). The model results are from the experiment R_Yves_s2 (Table S4), and the observation data are Gaillardet+HYBAM. The parameters are those that give the maximum R2 in Fig. 6d. The upper and lower panels show model-obs and log10(model)-log10(obs), respectively. The left (a, d), middle (b, e) and right (c, f) panels show rivers at the mid to high latitudes (if more than half of its river basin is located at or beyond 30° latitude), low latitudes (within 30° latitude) and over the whole globe, respectively. rlog2 and r2 are 0.33 and 0.16, respectively. The global total weathering flux Fw is 2.33×1012 mol yr−1.


3.3 Influence of the temperature-modulated soil production rate

Large seasonal changes in temperature can induce fractures in rocks and even form deep cracks in the surface soil layer (Liu et al., 2020), which may enhance the soil production rate at the base of the soil layer. Thus, the much weaker seasonal cycle in the tropical regions than at the higher latitudes (Fig. 6a) may be a factor to consider when calculating the erosion rate. To consider its influence, we assume that the soil production Pr is dependent on the amplitude of the seasonal cycle of surface temperature (defined as the difference between the maximum and minimum monthly temperatures), and the constant krp in Eq. (6) is now

(13) k rp = e ( A T - 24 ) / c - e - ( A T - 24 ) / c e ( A T - 24 ) / c + e - ( A T - 24 ) / c + b a ,

where AT is the amplitude of the seasonal cycle:

(14) A T = T max - T min ( K ) .

The constant 24 (K) in Eq. (13) is roughly the amplitude of the seasonal cycle at around 30° latitude (Fig. 8b). Across this critical amplitude, the soil production rate increases or decreases rapidly (Fig. 8a). Note that we have subjectively chosen to use a logistic function in Eq. (13) so as to make the soil production rate in the tropical region much lower than that in the extratropical region (Fig. 8b). This should be sufficient for the present purpose, which is to demonstrate whether the AT could have any significant impact on silicate weathering. The values of a, b and c determine the minimum values of krp and its variation with latitude, and a total of 12 combinations of a, b and c are tested (Table S5, Fig. 8a).

The forward calculation is repeated to search the parameter combinations (Table 2) that maximize R2 for all combinations of a, b and c. Results show that the best rlog2 and r2 are obtained when a, b and c are equal to 0.0244, 1.05 and 8 (yellow dotted line in Fig. 8a), respectively, when only considering the observation data from Gaillardet and 0.015, 2.3 and 1 (red solid line in Fig. 8a), respectively, when including HYBAM data (Fig. 8c–d). Both rlog2 and r2 are improved in terms of the bias in Fw. Compared to the new control model in Sect. 3.2, Fw, corresponding to the peaks of both rlog2 and r2, is slightly closer to the observational value (compare Fig. 8c–d to Fig. 6c–d). However, the values of rlog2 and r2, corresponding to the highest R2, are 0.201 and 0.204 when the HYBAM data are not included, remaining small. When the HYBAM data are included for model evaluation, the R2 value is higher but no better than that of the new control model (compare Fig. 8d to Fig. 6d). Nevertheless, this model is superior to the new control model in that the biases in both the tropical and extratropical regions are reduced this time (not shown).

Figure 8Possible effect of the seasonal cycle of surface temperature on the modeled silicate weathering flux. The soil production constant krp is assumed to depend on the amplitude of the seasonal cycle AT (b) according to the functions in panel (a). Yellow, red and blue colors in panel (a) correspond to the three columns of Table S5, respectively. Parameter c in Eq. (14) has values of 8, 4 and 2 for the dotted, dash-dotted and dashed lines, respectively. Panel (c) shows the rlog2 (orange), r2 (purple-red) and R2 (light blue) of all possible combinations of the parameters with the effect of AT considered. Instead of showing all the dots, only the envelopes (one for each color) are shown for the sake of clearness; the envelopes are obtained by curve fitting (cubic spline interpolation), and the data points used to do the fitting are still shown in the figure. Panel (d) is the same as panel (c), except that the HYBAM data are included in the observations. The results shown in panels (c) and (d) are from the experiment R_Yves_s2_td (Table 1). The black vertical line and grey zone show the observed Fw and its uncertainty range.

3.4 Implication of leached soil

Equation (9) tells us that local weathering flux is essentially the product of the erosion rate and the difference in the concentration of Ca and Mg cations between the bottom and top of the regolith. In the tropical regions, the cation concentration at the surface calculated by the Park20 model is near 0 where mountains are absent and is consistent with the distribution of leached soil (Figs. S6b and S8) defined in Sect. 2.2.7. The overestimation of tropical weathering fluxes (Fig. 1b) may thus indicate that the erosion rate in these regions is slower than that calculated by the model (Eq. 7). The cosmogenic nuclide analysis data do indicate lower erosion rates for the vast majority of rivers in the equatorial region than those from both TSS and the model calculation (Fig. 3b and Table S3). Figure 4 also shows that the results of the original Park20 model would be improved significantly if the observed erosion rates are used. Therefore, we think it is reasonable to slow down the erosion rate calculated by Eq. (7) when the areal fraction of leached soil in a grid box at mid to low latitudes (<30°) is greater than 20 %; the existence of such soil is an indication of slow erosion. The results are only slightly different if a different criterion is used because the areal fraction of leached soil is either very high or very low within the tropical region (Fig. S6b). Through a number of tests, it is found that the erosion rate by Eq. (7) should better be slowed down by an order of magnitude in these regions. The improvement of the erosion rate can be seen from Fig. 9b.

The model results are improved significantly with the simple change to the erosion rates above. The highest value of R2 reaches 0.73 and 0.93 when Gaillardet and Gaillardet+HYBAM, respectively, are used as observations (Fig. 10). Moreover, both rlog2 and r2 have high values (∼0.4) when R2 is at its maximum, higher than those obtained in Sect. 3.1 (Fig. 4c–d), where the model parameters were not optimized. Furthermore, the tropical bias is visibly reduced (compare Figs. 10 and 1b). These suggest that substantially slowing down the tropical erosion rates calculated by the Park20 model (Eq. 7) is an advisable choice. However, the appearance of leached soils is obviously a manifestation, not a reason, for the lower erosion rates. In addition, the distribution of leached soil is not available for the past or the future, just like the observed erosion rates tested in Sect. 3.1. Therefore, some other processes that are more fundamental and convenient than leached soil need to be found.

Figure 9Comparison between modeled and observed erosion rates. The observed erosion rates are from the cosmogenic nuclide analysis of Wittmann et al. (2020). Some data points (different in each panel but definitely less than 4) do not show up because the axis limits are set to relatively smaller values for the sake of clarity (compare this to Fig. 3b), but rlog2 and r2 as well as the linear correlation are calculated using all 18 data points. Colors represent regions where the basins are located (see also Fig. 3b), and the sizes represent the area of the basins. Panels (a)(e) differ in the way erosion is modeled. All the model calculations use R_Yves for runoff and s2 for the surface slope and use Eq. (7) with a different adjustment. (a) No adjustment. (b) The erosion rate is reduced by an order of magnitude where leached soils exist. (c) The erosion rate is reduced for a large LAI (from NCAR) according to Eq. (15). (d) Same as (c) except that the LAI is from the LPJ model. (e) Calculated by setting m in Eq. (7) to 0. In panels (c)(e), all the values are rescaled uniformly so that the global total erosion is 20 Gt yr−1.


Figure 10The effect of reduced erosion where leached soil is present in the modeled silicate weathering flux. Panels (a) and (b) show the rlog2 (orange), r2 (purple red) and R2 (light blue) of all possible combinations of the parameters. Similar to Fig. 8, only the envelopes and points used to fit the envelopes are shown. The results are from the experiment R_Yves_s2_soil (Table 1). The black vertical line and the grey zone show the observed Fw and its uncertainty range. Panel (c) shows the difference (model  observation) in silicate weathering fluxes for 81 large rivers, similar to Fig. 7a–c except that here the results corresponding to the highest R2 in panel (b) are shown. rlog2 and r2 are 0.56 and 0.37, respectively. The global total weathering flux is 2.29×1012 mol yr−1.


3.5 Influence of vegetation

It is observed that the distribution of leached soil (Fig. S6b) coincides with the flourishing of tropical vegetation (Fig. S6c, d) and may very well be a result of the latter. Although observations from arid regions indicate that the presence of vegetation significantly enhances mechanical erosion due to a rise in precipitation rates, mechanical erosion diminishes as vegetation cover increases in wet regions, owing to the dominant protective effects of vegetation (Mishra et al., 2019; Maffre et al., 2022). The presence of vegetation not only reduces the impact of raindrops on soil particles but also slows down the overland flow of water, decreasing the potential for soil detachment. Moreover, plant roots and organics contribute to soil cohesion and provide mechanical reinforcement (McMahon and Davies, 2018; Zeichner et al., 2021), thus reducing the overall likelihood of slope failures and landslides. Based on such thinking and the approximate coincidence between the distribution of leached soil and the region where the LAI is greater than 2, we design a way to modulate surface erosion with vegetation:

(15) E = E ( e - min 2 , LAI ) .

The basinal erosion rates calculated by the model in this way match those inferred from cosmogenic nuclide analysis better than when vegetation is not considered (Fig. 9c, d), substantiating the adjustment of erosion rate by vegetation. The erosion rates shown in Fig. 9c, d have been scaled up uniformly by changing ke in order for the global total erosion flux to retain a value of 20 Gt yr−1, although tests show that it only has a slight effect on the calculated silicate weathering fluxes. After considering the effect of vegetation, the maximum R2 can reach 0.84 (Fig. 11a, b), with the corresponding Fw being 2.8×1012 mol yr−1 (the experiment R_Yves_s2_LAI_global with the m1 method). Both rlog2 and r2 are also reasonably high (>0.3, Fig. S9a). Based on these results, here we propose that the suppression of erosion rates by vegetation was likely underestimated in previous studies on silicate weathering.

For the past or future, we will have to rely heavily on the model-simulated vegetation. However, the ability of current land models to simulate the vegetation and its response to climate change is still limited. Whether the effect of vegetation on silicate weathering can be properly considered is contingent upon how well the vegetation can be simulated. In one of the tests, the LAI simulated by the LPJ model (Fig. S10c, d; the experiment R_Yves_s2_LAI_old_global with the m2 method) was used. The results deteriorate substantially when the modeled global LAI was used; the maximum R2 merely attains a value of 0.64. This means that the defects in vegetation data cannot be compensated for by tuning other parameters in the weathering model. Therefore, getting better vegetation data by either reconstruction or model simulation is important for properly simulating silicate weathering of the past or future.

Figure 11Panels (a) and (b) show the envelopes of rlog2 (orange) and r2 (purple red) and their sums (light blue) of all possible combinations of the parameters. The effect of vegetation is considered by using the global LAI, and the global total erosion rate is scaled to 20 Gt yr−1. Only the cases with values greater than zero are shown. The black vertical line and the grey zone show the observed Fw and its uncertainty range. The results are from the experiment R_Yves_s2_LAI_global defined in Table 1.


3.6 Final parameters

A weathering model that adopts a parameterization for the effect of vegetation on erosion reduces the systematic error in the tropical region and is also easily applicable to other time periods. The Fw obtained by such a revised model is also closer to the most recently estimated global degassing flux (Müller et al., 2022). In this section, the optimal parameter set (that gives the highest R2) is provided for different combinations of runoff and surface slope (Table 3). The top five parameter sets with the highest R2 values (ranked based on the average R2 calculated for two sets of catchment-weathering flux measurements, R2a and R2b in Table 3) for each case are provided in Table 3. As can be seen, the parameter set highlighted in bold in Table 3 is amongst the best-performing parameter sets no matter which runoff or slope data are used. The weathering fluxes calculated using this set of parameters are much improved compared to those calculated using the original Park20 model in terms of both individual river basins (Fig. S11) and the global total (Fig. 11a, b). Subsequent calculations in this study are all based on this set of parameters unless otherwise stated.

Table 3Parameters chosen in the case of the global LAI. Bold fonts represent the selected best-performing parameters among different slope or runoff data.

a represents the fitting metrics with the observation of Gaillardet. b represents the observation data including the HYBAM network. c is the case that changes the erosion rate model by setting its sensitivity to the runoff to 0.

Download Print Version | Download XLSX

4 Discussion

4.1 Multiple effects of vegetation on silicate weathering

From the results presented in previous sections, we think that the silicate weathering fluxes calculated by previous models such as Park20 were systematically overestimated over the tropical region, and the overestimation was due at least in part to the overestimated erosion rate in this region. We thus propose that the overestimation in erosion was likely due to the underestimated effect of tropical vegetation on reducing erosion. The tests above show that this effect can be taken into account through a simple parameterization using the LAI, which can be obtained more easily by either reconstruction or model simulations (Binney et al., 2017; Krapp et al., 2021; Prentice et al., 2000; Prentice and Webb, 1998; Shao et al., 2018; Wang et al., 2008; Woillez et al., 2011; Yao et al., 2009; Andermann et al., 2022) for different time periods.

However, vegetation was generally thought to enhance silicate weathering by emitting organic acid (Caves Rugenstein et al., 2019; Berner, 2004, 1992), and the appearance of vegetation has been linked to the occurrence of a few ice ages in Earth's history (Lyla et al., 2011; Lenton et al., 2012). Here we are not arguing against such a mechanism and idea. Instead, we think that the ability of vegetation to enhance silicate weathering is universal and has been implicitly considered in model parameters such as the dissolution constant K in Eq. (2). In contrast, the effect of vegetation on soil protection could have been underestimated in silicate weathering models and could be geographically dependent. It is worth mentioning that Maffre et al. (2022) tested the effect of vegetation on slowing down soil erosion during the Devonian era, when vascular plants had just landed. Their work was more of a sensitivity study in that the observations (e.g., pCO2) could not provide vigorous constraints, which the basinal weathering fluxes used here do.

4.2 Influence of runoff

Some studies propose that the influence of runoff might have been overestimated in existing erosion rate frameworks. For instance, in a renowned model for erosion, BQART, sensitivity to runoff was adjusted downward from 0.5 to 0.31 (Syvitski and Milliman, 2007). Consequently, we did a simple test by assuming no correlation between erosion rate and runoff, i.e., setting m in Eq. (7) to 0. As can be seen from Fig. 9e, the model–data discrepancy is also reduced quite significantly by this method. For weathering calculation, the maximum R2 value obtained under this assumption is approximately 0.718 (Fig. S10), achieving a smaller improvement compared to the vegetation parameterization above but a notable one compared to other methodologies. The optimal parameter sets obtained for this test are provided in Table 3. A not unreasonable conjecture is that removing the dependence of erosion on runoff implicitly takes into account part of the influence of vegetation, since vegetation and runoff generally exhibit a positive correlation under contemporary conditions (see Figs. S3 and S6c). Because the factors affecting vegetation include not just precipitation (highly related to runoff) but also temperature, sunlight and pCO2, parameterizing erosion using vegetation is likely a superior way to using runoff.

4.3 Leached soil

It is worth emphasizing that the shielding effect on silicate weathering, although coincident with the distribution of leached soil, should be attributed to heavy vegetation as the former is likely a result of the retarding effect of the latter on erosion. This point can be partially inferred from the distribution of leached soils (Fig. 6b), LAI (Fig. 6c), surface slopes (Fig. S4a–b) and erosion rates (Fig. S4c–d). Leached soils should form in regions with relatively low erosion rates but no low silicate weathering rates. For example, regions with the lowest erosion rates are the desert regions where runoff is too small, but no leached soils form in these regions, probably due to the very low weathering rate. In some low-latitude regions where leached soils exist, the relatively low erosion rates (Fig. S4c–d) are consistent with the relatively small surface slope (i.e., flat terrain; Fig. S4a–b). However, in many other places (Fig. S6b), the erosion rates are relatively high (Fig. S4c–d) due to high runoff. These regions correspond more closely to the regions with a high LAI (Fig. 6c). Therefore, leached soils should be a result of a combination of vegetation development and relatively flat terrain, but the flat terrain is clearly not a necessary condition.

Godderis et al. (2008) also considered the effect of thick regolith cover on weathering by reducing the fluid that can reach the fresh bedrock, while we consider the effect herein by reducing the erosion rate. The two approaches agree with each other in that they both think that the weathering is transport-limited (i.e., fresh rocks are not exposed for weathering), but our approach is more direct and easier to apply to the paleo periods. This is because (1) the existence of the thick regolith cover is likely the result of weakened erosion (as described above) and (2) knowledge of regolith cover and thickness is unavailable for the past. Therefore, it seems better just to optimize the parameterization for the erosion rate directly, such as by considering the effect of land plants.

4.4 Sensitivity of global silicate weathering to climate change

The climate data from the 4×CO2 and LGM experiments (Sect. 2.2.2) are used to test the sensitivity of global silicate weathering to climate change. The land surface temperature increases from 278.4 K in the LGM to 286.6 K in the preindustrial and further to 301.1 K in the 4×CO2 experiment. Note that these changes are highly dependent on the climate model used but do not matter for the purpose here, which is to demonstrate the sensitivity of silicate weathering to climate changes between the Park20 model and the revised model in Sect. 3.6.

According to the Park20 model, the global silicate weathering flux Fw increases by 1.44 (46 %) from the LGM to PI and by 6.77 (149 %) from PI to the 4×CO2 situation (Table 4). For the revised model, Fw increases by 0.69 (32 %) from the LGM (experiment “RT_LGM_s2_LAI_global” with m3 in Table 1) to PI and by 4.38 (153 %) from PI to the 4×CO2 situation (experiment “RT_4CO2_s2_LAI_global” with m4 in Table 1). Thus, in terms of absolute values, the revised model is less sensitive to the climate, but in terms of relative values, the revised model is very similar to the original model. Because the relative change in the silicate weathering flux largely determines the relative change in pCO2 (see Eq. 2 of Goddéris et al., 2023), which determines the climate change, the weathering–climate sensitivity of the revised model is similar to that of the original model. However, due to the fact that Fw and its variations (in terms of absolute values) in the revised model are much smaller than before, other processes such as the burial of organic carbon may have been more important in Earth's carbon cycle than thought before.

Note that, although the LGM and 4×CO2 climates are used above to demonstrate the weathering–climate sensitivity, the timescales implied by these two experiments are only 10 000 and 100 years, respectively. These timescales are too short to be appropriate for the weathering models here, which assume that the weathering has reached a steady state; when the climate changes, vegetation may respond quickly (∼100 years), but the regolith layer and thus the weathering take a very long time to reach a new steady state.

Table 4Sensitivity of global silicate weathering to climate.

Download Print Version | Download XLSX

4.5 Caveats and future directions

The previously used measure for model–data discrepancy is rlog2, the maximization of which essentially optimizes the ratio between the model and the data. This measure has its advantages, but as we have shown above, such a measure cannot prevent the occurrence of a systematic error in the absolute difference between the model and the data (Fig. 1b). Optimizing r2, on the other hand, tends to underestimate Fw. We thus propose optimizing the sum of rlog2 and r2 (i.e., R2) so that Fw is nearest to the observation. It turns out that simply maximizing R2, while largely removing the systematic bias, would give very low values for both rlog2 and r2 (Fig. 7), meaning that changing the measure for model–data discrepancy alone cannot improve the model. To resolve the problem, certain physical processes have to be rectified, e.g., by invoking the influence of vegetation on erosion. A relatively satisfactory fit was finally obtained. However, R2 is still a subjective choice which may not be ideal. For example, R2 measures the overall degree of dispersion of the model-calculated fluxes around the observed fluxes, but it does not measure the correlation in spatial patterns. This may be one way to improve the measure for model–data discrepancy in the future.

The erosion rates derived from cosmogenic nuclides, as compared to those obtained from TSS, significantly alleviate the issue of overestimation of tropical weathering fluxes calculated by the model (see Fig. 4c, d). This improvement is likely due to the fact that the erosion rates derived from the cosmogenic nuclide analysis represent the erosion for long timescales, whereas TSS may have been substantially contaminated by human activities such as land use and deforestation (Hewawasam et al., 2003). Additionally, TSS could have been eroded primarily from near the river mouth, overestimating the erosion, or substantial deposition may have occurred on the way, underestimating the erosion (Wittmann et al., 2020). The major problem with the cosmogenic nuclide data is that they only cover a limited number of river basins. Both datasets also have the problem that neither of them can provide the spatial distribution of erosion distribution within river basins. Given the large area of many of the river basins (Fig. 3; the size of the circle represents the area of the river basins), uneven distribution of erosion within the river basin could have a great influence on the modeled weathering flux. More detailed observations of erosion rates and the related mechanisms are clearly needed in the future.

Although it seems that a simple parameterization for reducing the erosion rate by vegetation (Eq. 15) works well in improving model–data comparison, it must be noted that this may not be the sole or best resolution. The influence of vegetation on erosion may also depend on the local environment, which we have refrained from delving into further, primarily due to the plethora of uncertainties and insufficient constraints. Future observational evidence will be required to offer support for better parameterization. Another process that may be considered is the horizontal transport and deposition of materials. The current model is a one-dimensional one in which the regolith or soil comes from the bottom only, while in reality the soils can be eroded away easily and transported to another location, changing the local profile of the cation concentration.

5 Summary

A silicate weathering model that explicitly considers the regolith profile based on the formulation of GM09 and Park20 is studied in detail. This model has more than five underdetermined parameters which need to be constrained by the observed weathering fluxes for multiple river basins or watersheds over the globe. In doing so, the model–data discrepancy was normally measured by rlog2 (Eq. 11; larger values mean smaller discrepancies), and the parameter space was then searched to maximize rlog2. This method focuses more on minimizing the relative error (or discrepancy) than the absolute error. We demonstrate that the parameters determined this way tend to systematically overestimate the weathering fluxes over the tropical region, which leads to a significant overestimation of the global total flux Fw (Fig. 1). In addition, we show that such a problem is not due to uncertainties in the climate and surface slope data. We thus propose using R2 =rlog2+r2 as a new measure of model–data discrepancy, the maximization of which reduces both the relative and absolute errors in a more balanced way. By searching for the optimal parameters using this new measure, globally unbiased weathering fluxes are indeed obtained (Figs. 6 and 7c). However, the bias is removed by increasing the bias over the extratropical region (Fig. 7a) rather than by reducing the bias over the tropical region (Fig. 7b). Moreover, the model–data discrepancy is large; either rlog2 or r2 is small. Therefore, some other processes must be considered to reduce the bias over the tropical region and to reduce the model–data discrepancy.

The influence of the seasonal cycle of temperature on soil production is tested first based on the consideration that a stronger seasonal cycle can fracture and shatter rocks more easily. Little improvement can be achieved by such consideration (Fig. 8c, d). Next, the erosion rate is reduced in tropical regions where there are leached soils. It is found that the model–data discrepancy of silicate weathering fluxes is greatly reduced in this test (Fig. 10). Due to the fact that leached soil is the result, not the cause, of weakened erosion and the fact that the distribution of leached soil is almost coincident with that of forest, we propose that heavy vegetation is able to slow down erosion significantly. A simple parameterization is then put forward to consider the effect of vegetation on erosion by using the global LAI (Eq. 15). The LAI is used because it is relatively easy to obtain for other periods of Earth's history from Earth system model simulations. The Park20 model is revised to add this parameterization, and the model parameters are re-optimized using criterion R2 (Table 3). This revised model fits the observed weathering fluxes better than the original Park20 model (Fig. 11), and the modeled Fw is more consistent with both the observation and the most recently constructed global outgassing. Note that we are not against the idea of the evolution of land plants in Earth-enhanced silicate weathering: it is just that heavy vegetation could hinder silicate weathering by slowing down erosion over rainy regions. High precipitation and runoff in these regions would otherwise induce high erosion rates.

The revised model simulates a much smaller Fw than the original Park20 model. Correspondingly, the changes in Fw also become smaller under the same climate changes (Table 4), although the relative changes in Fw remain similar to the original model. If the revised model is reliable, it implies that the variations of other sinks of carbon such as organic carbon could have played a more important role than previously thought in the models. It will be interesting to see how the reconstruction of the Phanerozoic carbon cycle using models (e.g., Berner and Kothavala, 2001) will be impacted when the shielding effect of vegetation on silicate weathering as proposed herein is considered.

Code and data availability

The main code and data can be found at (Zuo, 2023).


The supplement related to this article is available online at:

Author contributions

HZ developed the model and ran all the simulations. HZ played the pivotal role in improving the model, to which YL contributed. YL designed the project and wrote the first draft together with HZ. All the authors contributed to the analyses and editing of the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We are grateful to Wenjing Liu for providing us with guidance on the source of the observational data and for revising the article.

Financial support

This work is supported by the National Natural Science Foundation of China (grant no. 42488201) and the National Key Research and Development Program of China (grant no. 2022YFF0800200). Zhifang Xu is supported by the National Key Research and Development 415 Program of China (grant no. 2020YFA0607700).

Review statement

This paper was edited by Lele Shu and reviewed by Pierre Maffre and two anonymous referees.


Adams, B. A., Whipple, K. X., Forte, A. M., Heimsath, A. M., and Hodges, K. V.: Climate controls on erosion in tectonically active landscapes, Sci. Adv., 6, eaaz3166,, 2020. 

Allen, J., Forrest, M., Hickler, T., Singarayer, J., Valdes, P., and Huntley, B.: Global vegetation patterns of the past 140,000 years, J. Biogeogr., 47, 2073–2090,, 2020. 

Amiotte Suchet, P., Probst, J.-L., and Ludwig, W.: Worldwide distribution of continental rock lithology: Implications for the atmospheric/soil CO2 uptake by continental weathering and alkalinity river transport to the oceans, Global Biogeochem. Cy., 17, 1038,, 2003. 

Andermann, T., Strömberg, C. A. E., Antonelli, A., and Silvestro, D.: The origin and evolution of open habitats in North America inferred by Bayesian deep learning models, Nat. Commun., 13, 4833,, 2022. 

Anderson, R.: Modeling the tor-dotted crests, bedrock edges, and parabolic profiles of high alpine surfaces of the Wind River Range, Wyoming, Geomorphology, 46, 35–58,, 2002. 

Berner, E. K. and Berner, R. A.: Global Environment: Water, Air, and Geochemical Cycles – Second Edition, 2, Princeton University Press,, 2012. 

Berner, R.: The Phanerozoic Carbon Cycle: CO2 and O2, Oxford Academic,, 2004. 

Berner, R., Lasaga, A., and Garrells, R.: The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years, Am. J. Sci, 283, 641–683,, 1983. 

Berner, R. A.: A model for atmospheric CO2 over Phanerozoic time, Am. J. Sci., 291, 339,, 1991. 

Berner, R. A.: Weathering, plants, and the long-term carbon cycle, Geochim. Cosmochim. Ac., 56, 3225–3231,, 1992. 

Berner, R. A. and Caldeira, K.: The need for mass balance and feedback in the geochemical carbon cycle, Geology, 25, 955–956,<0955:TNFMBA>2.3.CO;2, 1997. 

Berner, R. A. and Kothavala, Z.: GEOCARB III: A revised model of atmospheric CO2 over phanerozoic time, Am. J. Sci., 301, 182–204,, 2001. 

Binney, H., Edwards, M., Macias-Fauria, M., Lozhkin, A., Anderson, P., Kaplan, J. O., Andreev, A., Bezrukova, E., Blyakharchuk, T., Jankovska, V., Khazina, I., Krivonogov, S., Kremenetski, K., Nield, J., Novenko, E., Ryabogina, N., Solovieva, N., Willis, K., and Zernitskaya, V.: Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns, Quaternary Sci. Rev., 157, 80–97,, 2017. 

Blanckenburg, F., Bouchez, J., and Wittmann, H.: Earth surface erosion and weathering from the 10Be (meteoric)/9Be ratio, Earth Planet. Sc. Lett., 351–352, 295–305,, 2012. 

Bluth, G. and Kump, L.: Lithologic and climatologic controls of river chemistry, Geochim. Cosmochim. Ac., 58, 2341–2359,, 1994. 

Brantley, S. L., Bandstra, J., Moore, J., and White, A. F.: Modelling chemical depletion profiles in regolith, Geoderma, 145, 494–504,, 2008. 

Burke, B., Heimsath, A., and White, A.: Coupling chemical weathering with soil production across soil-landscapes, Earth Surf. Proc. Land., 32, 853–873,, 2007. 

Calabrese, S., Wild, B., Bertagni, M. B., Bourg, I. C., White, C., Aburto, F., Cipolla, G., Noto, L. V., and Porporato, A.: Nano- to global-scale uncertainties in terrestrial enhanced weathering, Environ. Sci. Technol., 56, 15261–15272,, 2022. 

Canadell, J. G., Monteiro, P. M. S., Costa, M. H., Cotrim da Cunha, L., Cox, P. M., Eliseev, A. V., Henson, S., Ishii, M., Jaccard, S., Koven, C., Lohila, A., Patra, P. K., Piao, S., Rogelj, J., Syampungani, S., Zaehle, S., and Zickfeld, K.: Global Carbon and Other Biogeochemical Cycles and Feedbacks, in: Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Intergovernmental Panel on Climate, C., Cambridge University Press, Cambridge, 673–816,, 2023. 

Carretier, S., Goddéris, Y., Delannoy, T., and Rouby, D.: Mean bedrock-to-saprolite conversion and erosion rates during mountain growth and decline, Geomorphology, 209, 39–52,, 2014. 

Carretier, S., Goddéris, Y., Martinez, J., Reich, M., and Martinod, P.: Colluvial deposits as a possible weathering reservoir in uplifting mountains, Earth Surf. Dynam., 6, 217–237,, 2018. 

Caves Rugenstein, J., Ibarra, D., Zhang, S., Planavsky, N., and Blanckenburg, F.: Isotope mass-balance constraints preclude that mafic weathering drove Neogene cooling, P. Natl. Acad. Sci. USA, 118, e2026345118,, 2021. 

Caves Rugenstein, J. K., Ibarra, D. E., and von Blanckenburg, F.: Neogene cooling driven by land surface reactivity rather than increased weathering fluxes, Nature, 571, 99–102,, 2019. 

Danabasoglu, G.: NCAR CESM2 model output prepared for CMIP6 CMIP, WCRP [data set],, 2019. 

Dannhaus, N., Wittmann, H., Krám, P., Christl, M., and Blanckenburg, F.: Catchment-wide weathering and erosion rates of mafic, ultramafic, and granitic rock from cosmogenic meteoric 10 Be/ 9 Be ratios, Geochim. Cosmochim. Ac., 222, 618–641,, 2017. 

D'Antonio, M., Ibarra, D., and Boyce, C.: Land plant evolution decreased, rather than increased, weathering rates, Geology, 48, 29–33,, 2019. 

Davy, P. and Crave, A.: Upscaling local-scale transport processes in large-scale relief dynamics, Phys. Chem. Earth., 25, 533–541,, 2000. 

Dellinger, M., Gaillardet, J., Bouchez, J., Calmels, D., Louvat, P., Dosseto, A., Gorge, C., Alanoca, L., and Maurice, L.: Riverine Li isotope fractionation in the Amazon River basin controlled by the weathering regimes, Geochim. Cosmochim. Ac., 164, 71–93,, 2015. 

Dessert, C., Dupré, B., Gaillardet, J., François, L. M., and Allègre, C. J.: Basalt weathering laws and the impact of basalt weathering on the global carbon cycle, Chem. Geol., 202, 257–273,, 2003. 

Dietrich, W., Reiss, R., Hsu, M.-L., and Montgomery, D.: A process-based model for colluvial soil depth and shallow landsliding using digital elevation data, Hydrol. Process., 9, 383–400,, 1995. 

Dixon, J., Heimsath, A., and Amundson, R.: Critical role of climate and saprolite weathering in landscape evolution, Earth Surf. Proc. Land., 34, 1507–1521,, 2009. 

Edmond, J. M., Palmer, M. R., Measures, C. I., Grant, B., and Stallard, R. F.: The fluvial geochemistry and denudation rate of the Guayana Shield in Venezuela, Colombia, and Brazil, Geochim. Cosmochim. Ac., 59, 3301–3325,, 1995. 

Emerson, S. and Hedges, J.: Chemical Oceanography and the Marine Carbon Cycle, Cambridge University Press, Cambridge,, 2008. 

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D.: The shuttle radar topography mission, Rev. Geophys., 45, RG2004,, 2007. 

Fekete, B., Vörösmarty, C. J., and Grabs, W.: Highresolution fields of global runoff combining river discharge and simulated water balances, Global Biogeochem. Cy., 16, 15-1–15-10,, 2002. 

FAO/IIASA/ISRIC/ISSCAS/JRC: Harmonized World Soil Database (version 1.2), FAO, Rome, Italy and IIASA, Laxenburg, Austria [data set], (last access: 1 December 2021), 2012. 

France-Lanord, C. and Derry, L. A.: Organic carbon burial forcing of the carbon cycle from Himalayan erosion, Nature, 390, 65–67,, 1997. 

Gabet, E. J. and Mudd, S. M.: A theoretical model coupling chemical weathering rates with denudation rates, Geology, 37, 151–154,, 2009. 

Gaillardet, J., Dupré, B., Louvat, P., and Allègre, C. J.: Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers, Chem. Geol., 159, 3–30,, 1999. 

Galy, A. and France-Lanord, C.: Weathering processes in the Ganges–Brahmaputra basin and the riverine alkalinity budget, Chem. Geol., 159, 31–60,, 1999. 

Gasparini, N., Whipple, K., and Bras, R.: Predictions of steady state and transient landscape morphology using sediment-flux-dependent river incision models, J. Geophys. Res., 112, F03S09,, 2007. 

Gerlach, T.: Volcanic versus anthropogenic carbon dioxide, Eos Trans. Agu, 92, 201–202,, 2011. 

Ghiggi, G., Humphrey, V., Seneviratne, S. I., and Gudmundsson, L.: GRUN: an observation-based global gridded runoff dataset from 1902 to 2014, Earth Syst. Sci. Data, 11, 1655–1674,, 2019. 

Gibbs, M., Bluth, G., Fawcett, P., and Kump, L.: Global chemical erosion over the last 250 MY: Variations due to changes in paleogeography, paleoclimate, and paleogeology, Am. J. Sci, 299, 611–651,, 1999. 

Godderis, Y., Donnadieu, Y., Tombozafy, M., and Dessert, C.: Shield effect on continental weathering: Implication for climatic evolution of the Earth at the geological timescale, Geoderma, 145, 439–448,, 2008. 

Goddéris, Y., Donnadieu, Y., Carretier, S., Aretz, M., Dera, G., Macouin, M., and Regard, V.: Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering, Nat. Geosci., 10, 382–386,, 2017. 

Goddéris, Y., Donnadieu, Y., and Mills, B. J. W.: What models tell us about the evolution of carbon sources and sinks over the Phanerozoic, Annu. Rev. Earth Pl. Sc., 51, 471–492,, 2023. 

Gruber, C., Zhu, C., Georg, R. B., Zakon, Y., and Ganor, J.: Resolving the gap between laboratory and field rates of feldspar weathering, Geochim. Cosmochim. Ac., 147, 90–106,, 2014. 

Harel, M. A., Mudd, S. M., and Attal, M.: Global analysis of the stream power law parameters based on worldwide 10Be denudation rates, Geomorphology, 268, 184–196,, 2016. 

Harris, I., Jones, P., Osborn, T., and Lister, D.: Updated high-resolution grids of monthly climatic observations – The CRU TS3.10 Dataset, Int. J. Climatol., 34, 623–642,, 2014. 

Hartmann, J. and Moosdorf, N.: The new global lithological map database GLiM: A representation of rock properties at the Earth surface, Geochem. Geophy. Geosy., 13, Q12004,, 2012. 

Hartmann, J., Jansen, N., Dürr, H. H., Kempe, S., and Köhler, P.: Global CO2-consumption by chemical weathering: What is the contribution of highly active weathering regions?, Global Planet. Change, 69, 185–194,, 2009. 

Hartmann, J., Moosdorf, N., Lauerwald, R., Hinderer, M., and West, A. J.: Global chemical weathering and associated P-release – The role of lithology, temperature and soil properties, Chem. Geol., 363, 145–163,, 2014. 

Heimsath, A. and Korup, O.: Quantifying rates and processes of landscape evolution, Earth Surf. Proc. Land., 37, 249–251,, 2012. 

Heimsath, A., Dietrich, W., Nishiizumi, K., and Finkel, R.: Cosmogenic nuclides, topography, and the spatial variation of soil depth, Geomorphology, 27, 151–172,, 1999. 

Heimsath, A., Fink, D., and Hancock, G.: The “humped” soil production function: Eroding Arnhem Land, Australia, Earth Surf. Proc. Land., 34, 1674–1684,, 2009. 

Heimsath, A. M., Dietrich, W. E., Nishiizumi, K., and Finkel, R. C.: The soil production function and landscape equilibrium, Nature, 388, 358–361,, 1997. 

Hewawasam, T., von Blanckenburg, F., Schaller, M., and Kubik, P.: Increase of human over natural erosion rates in tropical highlands constrained by cosmogenic nuclides, Geology, 31, 597–600,<0597:IOHONE>2.0.CO;2, 2003. 

Hilton, R. G. and West, A. J.: Mountains, erosion and the carbon cycle, Nat. Rev. Earth Env., 1, 284–299,, 2020. 

Howard, A.: A detachment-limited model of drainage-basin Evolution, Water Resour. Res., 30, 2261–2285,, 1994. 

Hu, Y., Teng, F.-Z., Plank, T., and Chauvel, C.: Potassium isotopic heterogeneity in subducting oceanic plates, Sci. Adv., 6, eabb2472,, 2020. 

Ibarra, D. E., Rugenstein, J. K. C., Bachan, A., Baresch, A., Lau, K. V., Thomas, D. L., Lee, J.-E., Boyce, C. K., and Chamberlain, C. P.: Modeling the consequences of land plant evolution on silicate weathering, Am. J. Sci, 319, 1–43,, 2019. 

Kalderon-Asael, B., Katchinoff, J., Planavsky, N., Hood, A., Dellinger, M., Bellefroid, E., Jones, D., Hofmann, A., Ossa, F., Macdonald, F., Wang, C., Isson, T., Murphy, J., Higgins, J., West, A. J., Wallace, M., Asael, D., and Pogge von Strandmann, P.: A lithium-isotope perspective on the evolution of carbon and silicon cycles, Nature, 595, 394–398,, 2021. 

Krapp, M., Beyer, R. M., Edmundson, S. L., Valdes, P. J., and Manica, A.: A statistics-based reconstruction of high-resolution global terrestrial climate for the last 800,000 years, Sci. Data, 8, 228,, 2021. 

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61,, 2014. 

Larsen, I. J., Almond, P. C., Eger, A., Stone, J. O., Montgomery, D. R., and Malcolm, B.: Rapid soil production and weathering in the Southern Alps, New Zealand, Science, 343, 637–640,, 2014. 

Lawrence, P. and Chase, T.: Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0), J. Geophys. Res., 112, G01023,, 2007. 

Lécuyer, C.: Seawater residence times of some elements of geochemical interest and the salinity of the oceans, Bulletin de la Société Géologique de France, 187, 245–260,, 2016. 

Lee, C.-T. A., Jiang, H., Dasgupta, R., and Torres, M.: A framework for understanding whole-Earth carbon cycling, in: Deep Carbon, 313–357,, 2019. 

Lenton, T. M., Crouch, M., Johnson, M., Pires, N., and Dolan, L.: First plants cooled the Ordovician, Nat. Geosci., 5, 86–89,, 2012. 

Li, S., Li, W., Beard, B. L., Raymo, M. E., Wang, X., Chen, Y., and Chen, J.: K isotopes as a tracer for continental weathering and geological K cycling, P. Natl. Acad. Sci. USA, 116, 8740–8745,, 2019. 

Li, X., Hu, Y., Yang, J., Wei, M., Guo, J., Lan, J., Lin, Q., Yuan, S., Zhang, J., Wei, Q., Liu, Y., Nie, J., Xia, Y., and Hu, S.: Climate variations in the past 250 million years and contributing factors, Paleoceanogr. Paleoclimatol., 38, e2022PA004503,, 2023. 

Liu, Y., Yang, J., Bao, H., Shen, B., and Hu, Y.: Large equatorial seasonal cycle during Marinoan snowball Earth, Sci. Adv., 6, eaay2471,, 2020. 

Lyla, T., Steve, B., Jonathan, L., and David, J. B.: Modeling the evolutionary rise of ectomycorrhiza on sub-surface weathering environments and the geochemical carbon cycle, Am. J. Sci., 311, 369,, 2011. 

Maffre, P., Ladant, J.-B., Moquet, J.-S., Carretier, S., Labat, D., and Goddéris, Y.: Mountain ranges, climate and weathering. Do orogens strengthen or weaken the silicate weathering carbon sink?, Earth Planet. Sc. Lett., 493, 174–185,, 2018. 

Maffre, P., Godderis, Y., Pohl, A., Donnadieu, Y., Carretier, S., and Hir, G.: The complex response of continental silicate rock weathering to the colonization of the continents by vascular plants in the Devonian, Am. J. Sci, 322, 461–492,, 2022. 

Maher, K.: The dependence of chemical weathering rates on fluid residence time, Earth Planet. Sc. Lett., 294, 101–110,, 2010. 

Maher, K. and Chamberlain, C. P.: Hydrologic regulation of chemical weathering and the geologic, Science, 343, 1502–1504,, 2014. 

McMahon, W. J. and Davies, N. S.: Evolution of alluvial mudrock forced by early land plants, Science, 359, 1022–1024,, 2018. 

Meybeck, M.: Global chemical weathering of surficial rocks estimated from river dissolved loads, Am. J. Sci., 287, 401–428,, 1987. 

Milliman, J. and Farnsworth, K.: River Discharge to the Coastal Ocean – A Global Synthesis, Cambridge University Press,, 2011. 

Milliman, J. and Syvitski, J.: Geomorphic tectonic control of sediment discharge to ocean – The importance of small mountainous rivers, J. Geol., 100, 525–544,, 1991. 

Milliman, J. D., Rutkowski, C., and Meybeck, M.: River discharge to the sea; a global river index (GLORI), loicz reports & studies, no. 2, (last access: 15 June 2022), 1995. 

Mills, B. J. W., Donnadieu, Y., and Goddéris, Y.: Spatial continuous integration of Phanerozoic global biogeochemistry and climate, Gondwana Res., 100, 73–86,, 2021. 

Mishra, A., Placzek, C., and Jones, R.: Coupled influence of precipitation and vegetation on millennial-scale erosion rates derived from 10Be, PLOS ONE, 14, e0211325,, 2019. 

Moon, S., Chamberlain, C. P., and Hilley, G. E.: New estimates of silicate weathering rates and their uncertainties in global rivers, Geochim. Cosmochim. Ac., 134, 257–274,, 2014. 

Moquet, J.-S., Crave, A., Viers, J., Seyler, P., Armijos, E., Bourrel, L., Chavarri, E., Lagane, C., Laraque, A., Lavado, W., Pombosa, R., Noriega, L., Vera, A., and Guyot, J.-L.: Chemical weathering and atmospheric/soil CO2 uptake in the Andean and Foreland Amazon basins, Chem. Geol., 287, 1–26,, 2011. 

Moquet, J.-S., Guyot, J.-L., Crave, A., Viers, J., Filizola Jr, N., Martinez, J., Oliveira, T., Hidalgo Sánchez, L., Lagane, C., Lavado, W., Noriega, L., and Pombosa, R.: Amazon River dissolved load: temporal dynamics and annual budget from the Andes to the ocean, Environ. Sci. Pollut. R., 23, 11405–11429,, 2016. 

Moquet, J.-S., Guyot, J.-L., Morera, S., Crave, A., Rau, P., Vauchel, P., Lagane, C., Sondag, F., Lavado, W., Pombosa, R., and Martinez, J.: Temporal variability and annual budget of inorganic dissolved matter in Andean Pacific Rivers located along a climate gradient from northern Ecuador to southern Peru, Cr. Geosci., 350, 76–87,, 2018. 

Müller, R. D., Mather, B., Dutkiewicz, A., Keller, T., Merdith, A., Gonzalez, C. M., Gorczyk, W., and Zahirovic, S.: Evolution of Earth's tectonic carbon conveyor belt, Nature, 605, 629–639,, 2022. 

Muñoz Sabater, J.: ERA5-Land monthly averaged data from 1950 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set],, 2019. 

Olson, S., Jansen, M. F., Abbot, D. S., Halevy, I., and Goldblatt, C.: The effect of ocean salinity on climate and its implications for Earth's habitability, Geophys. Res. Lett., 49, e2021GL095748,, 2022. 

Park, Y., Maffre, P., Godderis, Y., Macdonald, F., Anttila, E., and Swanson-Hysell, N.: Emergence of the Southeast Asian islands as a driver for Neogene cooling, P. Natl. Acad. Sci. USA, 117, 25319–25326,, 2020. 

Phillips, J.: The convenient fiction of steady-state soil thickness, Geoderma, 156, 389–398,, 2010. 

Prentice, I. C. and Webb III, T.: BIOME 6000: reconstructing global mid-Holocene vegetation patterns from palaeoecological records, J. Biogeogr., 25, 997–1005,, 1998. 

Prentice, I. C., Jolly, D., and BIOME 6000 participants: Mid-Holocene and glacial-maximum vegetation geography of the northern continents and Africa, J. Biogeogr., 27, 507–519,, 2000. 

Quye-Sawyer, J., Whittaker, A. C., and Roberts, G. G.: Calibrating fluvial erosion laws and quantifying river response to faulting in Sardinia, Italy, Geomorphology, 370, 107388,, 2020. 

Raymo, M. E. and Ruddiman, W. F.: Tectonic forcing of late Cenozoic climate, Nature, 359, 117–122, 1992. 

Riebe, C. S., Kirchner, J. W., and Finkel, R. C.: Erosional and climatic effects on long-term chemical weathering rates in granitic landscapes spanning diverse climate regimes, Earth Planet. Sc. Lett., 224, 547–562,, 2004. 

Royden, L. and Taylor Perron, J.: Solutions of the stream power equation and application to the evolution of river longitudinal profiles, J. Geophys. Res-Earth., 118, 497–518,, 2013. 

Rudnick, R. and Gao, S.: Composition of the Continental Crust, Treatise on Geochemistry, 1–64 pp.,, 2003. 

Scotese, C. R. and Wright, N.: PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic, PALEOMAP Project [data set], (last access: 20 April 2019), 2018. 

Shao, Y., Anhäuser, A., Ludwig, P., Schlüter, P., and Williams, E.: Statistical reconstruction of global vegetation for the last glacial maximum, Global Planet. Change, 168, 67–77,, 2018. 

Small, E., Anderson, R., and Hancock, G.: Estimates of the rate of regolith production using 10Be and 26Al from an alpine hillslope, Geomorphology, 27, 131–150,, 1999. 

Stallard, R. F.: River Chemistry, Geology, Geomorphology, and Soils in the Amazon and Orinoco Basins, in: The Chemistry of Weathering, edited by: Drever, J. I., Springer Netherlands, Dordrecht, 293–316,, 1985. 

Stallard, R. F. and Edmond, J. M.: Geochemistry of the Amazon: 1. Precipitation chemistry and the marine contribution to the dissolved load at the time of peak discharge, J. Geophys. Res-Oceans., 86, 9844–9858,, 1981. 

Stallard, R. F. and Edmond, J. M.: Geochemistry of the Amazon: 2. The influence of geology and weathering environment on the dissolved load, J. Geophys. Res-Oceans., 88, 9671–9688,, 1983. 

Strudley, M., Murray, A. B., and Haff, P.: Emergence of pediments, tors, and piedmont junctions from a bedrock weathering-regolith thickness feedback, Geology, 34, 805–808,, 2006. 

Suchet, P. and Probst, J.-L.: A global model for present-day atmospheric/soil CO2 consumption by chemical erosion of continental rocks (GEM-CO2), Tellus B., 47, 273–280,, 2002. 

Syvitski, J. and Milliman, J.: Geology, geography, and humans battle for dominance over the delivery of fluvial sediment to the coastal ocean, J. Geol., 115, 1–19,, 2007. 

Walker, J. C. G., Hays, P. B., and Kasting, J. F.: A negative feedback mechanism for the long-term stabilization of Earth's surface temperature, J. Geophys. Res.-Oceans, 86, 9776–9782,, 1981. 

Wang, G., Feng, X., Han, J., Zhou, L., Tan, W., and Su, F.: Paleovegetation reconstruction using δ13C of Soil Organic Matter, Biogeosciences, 5, 1325–1337,, 2008. 

West, A. J.: Thickness of the chemical weathering zone and implications for erosional and climatic drivers of weathering and for carbon-cycle feedbacks, Geology, 40, 811–814,, 2012. 

West, A. J., Galy, A., and Bickle, M.: Tectonic and climatic controls on silicate weathering, Earth Planet. Sc. Lett., 235, 211–228,, 2005. 

Whipple, K., Heimsath, A., and DiBiase, R.: Soil production limits and the transition to bedrock-dominated landscapes, Nat. Geosci., 5, 210–214,, 2012. 

White, A. F. and Blum, A. E.: Effects of climate on chemical weathering in watersheds, Geochim. Cosmochim. Ac., 59, 1729–1747,, 1995. 

White, A. F. and Brantley, S. L.: The effect of time on the weathering of silicate minerals: Why do weathering rates differ in the laboratory and field?, Chem. Geol., 202, 479–506,, 2003. 

Wittmann, H., Oelze, M., Gaillardet, J., Garzanti, E., and Blanckenburg, F.: A global rate of denudation from cosmogenic nuclides in the Earth's largest rivers, Earth-Sci. Rev., 204, 103147,, 2020. 

Wittmann, H., Blanckenburg, F., Bourgoin, L., Guyot, J.-L., Filizola Jr., N., and Kubick, P. W.: Sediment production and delivery in the Amazon River basin quantified by in situ produced cosmogenic nuclides and recent river loads, Geol. Soc. Am. Bull., 123, 934–950,, 2011. 

Wittmann, H., Blanckenburg, F., Dannhaus, N., Bouchez, J., Gaillardet, J., Guyot, J.-L., Bourgoin, L., Roig, H., Filizola Jr, N., and Christl, M.: A test of the cosmogenic 10 Be(meteoric)/9 Be proxy for simultaneously determining basin-wide erosion rates, denudation rates, and the degree of weathering in the Amazon basin, J. Geophys. Res-Earth., 120, 2498–2528,, 2015. 

Woillez, M.-N., Kageyama, M., Krinner, G., de Noblet-Ducoudré, N., Viovy, N., and Mancip, M.: Impact of CO2 and climate on the Last Glacial Maximum vegetation: results from the ORCHIDEE/IPSL models, Clim. Past, 7, 557–577,, 2011. 

Yao, Y.-F., Bera, S., Ferguson, D. K., Mosbrugger, V., Paudayal, K. N., Jin, J.-H., and Li, C.-S.: Reconstruction of paleovegetation and paleoclimate in the Early and Middle Eocene, Hainan Island, China, Clim. Change, 92, 169–189,, 2009. 

Zeichner, S. S., Nghiem, J., Lamb, M. P., Takashima, N., de Leeuw, J., Ganti, V., and Fischer, W. W.: Early plant organics increased global terrestrial mud deposition through enhanced flocculation, Science, 371, 526–529,, 2021.  

Zhang, M., Liu, Y., Zhu, J., Wang, Z., and Liu, Z.: Impact of dust on climate and AMOC during the Last Glacial Maximum simulated by CESM1.2, Geophys. Res. Lett., 49, e2021GL096672,, 2022. 

Zhang, S., Bai, X., Zhao, C., Tan, Q., Yun, L., Wang, J., Li, L., Wu, L., Chen, F., Li, C., Deng, Y., Yang, Y., and Xi, H.: Global CO2 consumption by silicate rock chemical weathering: Its past and future, Earths Future, 9, e2020EF001938,, 2021. 

Zhang, Y., Mills, B., Yang, T., He, T., and Zhu, M.: Simulating the long-term carbon cycle in the Phanerozoic: current status and future developments, Chinese Journal, 68, 1580–1592,, 2022. 

Zuo, H.: zuohaoyue1/Silicate-weathering-model: Silicate weathering model through fitting parameters (Version v1), Zenodo [code],, 2023. 

Short summary
Compared to the silicate weathering fluxes measured at large river basins, the current models tend to systematically overestimate the fluxes over the tropical region, which leads to an overestimation of the global total weathering flux. The most possible cause of such bias is found to be the overestimation of tropical surface erosion, which indicates that the tropical vegetation likely slows down physical erosion significantly. We propose a way of taking this effect into account in models.