Response to reviewers’ comments

Response to Reviewer #1:

[Comment 1] In this paper, the authors hypothesize that the mismatch between modelled and measured global silicate-weathering rates are due to uncertainties in the underlying erosion models. The authors therefore present a global erosion-rate model trained against cosmogenic-nuclide derived erosion rates using a random forest approach. They show that this model can reproduce erosion rates much better than a basic stream-power-erosion approach. The authors then incorporate their predicted erosion patterns into an existing steady-state silicate-weathering model and demonstrate that the new product can predict both the total global weathering flux and the spatial distribution of weathering rates better than previous models.

I found the study to be very well motivated, well written, and convincingly argued. The analysis is sound and provides a significant contribution for modelling global silicate-weathering fluxes. The figures and text are of a high quality. In particular, I appreciate the explanations of complex terms and approaches (e.g. the SHAP work) for the unfamiliar reader. I only have a few relatively minor comments that may be interesting to consider before publication.

Response: Thank you very much for your thoughtful and constructive review of our manuscript. We are gratified by your positive assessment of the study's motivation and the clarity of our explanations. Your comments have helped us identify areas where our terminology was ambiguous or where the data presentation required more rigorous quantification.

[Comment 2] L478 – 496: You write that Fig. 9 “clearly illustrates the trade-off”, that the value in Fig. 9a is “far to the right of the observed range” and that its value within the square-band is “low”. To my eyes, the peak is just to the right of the optimum and there are still reasonable values within the grey band. Next, you write that in Fig. 9b, “peaks of all three metric curves are now located squarely within the grey observational band”. For the red and blue curves that is true, but the peak of the yellow curve is clearly not in the square band and actually looks quite close to the location in Fig 9b. Can you revise this text so that the text matches the data more clearly?

Response: We sincerely thank the reviewer for this rigorous and insightful evaluation of Figure 9. We appreciate the reviewer pointing out that our descriptions of the curves' positions relative to the observational grey band were not sufficiently precise.

We agree that the term "far to the right" in L478-496 is subjective. In Figure 9a (SPIM), the R_log^2 peak (representing the optimal fit to the spatial pattern of observations evaluated with the original metric Park et al. used in the paper) occurs at a global flux of ~4.2 × 10^12  mol 〖yr〗^(-1) which is approximately 1.7 times of the observation-based central estimate of 2.5 × 10^12  mol 〖yr〗^(-1). While some R_log^2 values within the grey band are technically "reasonable" (0.41 compared to its optimal of 0.53), the model's best spatial performance (e.g., the peak) cannot be achieved without overestimating the global total flux.

The reviewer is correct that the peak of the yellow curve (R_log^2) in Figure 9b sits at ~3.5 × 10^12  mol 〖yr〗^(-1), which is indeed outside the 2.0-3.0 × 10^12  mol 〖yr〗^(-1)  grey band. Our previous statement that "all three metric curves" were "squarely within" the band was inaccurate.

We have revised the text to clarify that the significant improvement in MErSiM (Fig. 9b) is the convergence of the metrics. Specifically, the peak of the joint metric (R^2+R_log^2, blue curve) and the standard R^2 (red curve) now fall within the observational range. The final parameter set we used for simulation is those who optimize R^2+R_log^2.

[Comment 3] Perhaps, you can also explicitly quantify (1) the global weathering fluxes predicted by the two different models (MErSiM versu Park20) including an estimate of uncertainty for both, and (2) the difference between each one of these predictions and the measured total weathering fluxes.

Response: In the original manuscript, we have explicitly quantified the predicted flux in Section 3.5.1 and in Table 3 where the sensitivity test is presented, but the reviewer is right that this quantification could be mentioned in this section. Following the reviewer’s suggestion, we have added a text summary explicitly comparing the model predictions to observations. To account for uncertainty in parameter selection, we consider all parameter sets whose metric value is within 0.05 of the maximum. The prediction uncertainty is then quantified by the variability of model outputs across these acceptable parameter sets.

Measured target: 2.5±0.5 × 10^12  mol 〖yr〗^(-1)

Park20 (GM09-SPIM): 4.54±0.4 × 10^12  mol 〖yr〗^(-1) (Error: 2.04 × 10^12  mol 〖yr〗^(-1), 82% bias).

MErSiM (Revised): 3.11±0.3 × 10^12  mol  〖yr〗^(-1) (Error: 0.61 × 10^12  mol 〖yr〗^(-1), 24% bias).

[Comment 4] Finally, in the section L478 – 496 you claim that these results speak to the trade-off between matching total weathering fluxes versus the pattern of weathering fluxes. I did not understand how Figure 9 relates to (or informs) the spatial pattern of silicate weathering fluxes.

Response: We appreciate the reviewer's request for clarification on the relationship between the metrics in Figure 9 and the "spatial pattern" of weathering. In our evaluation framework, the performance metrics (R^2 and R_log^2) shown in Figure 9 are not calculated for a single location, but represent the agreement between modeled and observed weathering fluxes across 81 major river basins distributed globally. Because these 81 basins encompass a vast range of tectonic settings (from active orogens like the Andes to stable cratons like the Amazon) and climatic zones (from the tropics to high latitudes), a high R^2 value indicates that the model successfully reproduces the observed spatial heterogeneity and relative differences in weathering rates among these diverse regions. Therefore, the "model performance" plotted on the y-axis of Figure 9 is a direct proxy for how well the model captures the global spatial pattern of silicate weathering, while the value on the x-axis indicates the performance in matching the total weathering flux. A tradeoff is thus obtained if the maximum y value is obtained for an x value within the grey band. We have clarified this definition in the revised manuscript to avoid ambiguity.

[Comment 5] L497 – 504: Similar to above, the text here seems a bit more pushed than warranted based on Fig 10. I agree that the errors are clearly smaller in the RF model, but there are still quite large errors, in particular in the tropical regions. When you write things like “large” errors and is “much less” severe, can you explain quantitatively what you understand by a large and small error. How much were errors reduced on average etc.?

Response: We thank the reviewer for this constructive critique. We agree that using qualitative descriptors like "large" or "much less severe" without clear benchmarks can be subjective. In the revised manuscript, we have explicitly defined "large errors" as those exceeding 5 × 10^10  mol 〖yr〗^(-1).

To provide a clearer comparison, we have calculated the following metrics across the 81 major basins. First, MErSiM reduced the Mean Absolute Error (MAE) across all 81 basins by approximately 35% compared to the original GM09-SPIM (1.77 × 10^10  mol 〖yr〗^(-1) compared to 2.72 × 10^10  mol 〖yr〗^(-1)). Second, the number of basins exhibiting "large" positive errors (>5 × 10^10  mol 〖yr〗^(-1)) dropped from 10 in GM09-SPIM to 4 in MErSiM. Third, specifically in the tropics (basins located between 23S-23N), where GM09-SPIM errors often exceeded 7 × 10^10  mol 〖yr〗^(-1), MErSiM reduced the average basin-scale bias by 90% (0.41 × 10^10  mol 〖yr〗^(-1) compared to 4.12 × 10^10  mol 〖yr〗^(-1)). We have added these quantitative description in the explanation paragraph of Fig. 10.

We appreciate the reviewer for pointing out that significant errors remain in the tropics. We have added a paragraph noting that these persistent biases likely stem from uncertainties in the cation concentrations of sedimentary and metamorphic rocks, which remain as adjustable bulk parameters in the current framework.

[Comment 6] I was confused by the discussion about the effect of CO2 increases. You explain your models small response to the CO2 increase by widespread chemostasis and a resulting runoff insensitivity of the weathering flux. My understanding is exactly the opposite: To me, chemostasis means that concentrations remain constant with changes in runoff – contrary to what is suggested in L568 (e.g. Godsey et al., 2019; Godsey et al., 2009). Under these conditions (common in active mountains (Godsey et al., 2019)) the weathering fluxes should strongly increase with runoff. Hence, the system should have a high sensitivity to increased CO2 and runoff. Something is off here, and maybe it is just about the term chemostasis? I guess your model also has a negative temperature response, does that matter here?

Response: We are very grateful to the reviewer for pointing out the potential confusion regarding the term "chemostasis." We agree that in the hydrological community (e.g., Godsey et al., 2009, 2019), chemostasis typically refers to a state where solute concentrations remain nearly constant despite changes in discharge, which implies that weathering fluxes should scale linearly with runoff.

In the original manuscript, we used "chemostatic" in the context of the Damköhler number framework (Maher and Chamberlain, 2014) to describe a state of weathering flux saturation. We acknowledge that this usage is ambiguous and contradicts the definition favored by the reviewer.

In our revised model, under greenhouse conditions (), the system enters a regime where the fluid travel time (τ_f) is much shorter than the reaction equilibration time (τ_eq). This leads to a "sluggish" response, causing the total weathering flux to plateau despite in higher temperature scenario.

The reviewer correctly intuited the role of the negative temperature response. In our data-driven erosion module, warming in high-latitude or high-altitude regions can move the environment out of the optimal "frost-cracking window (Anderson and Anderson, 2010; Zhao et al., 2026)," leading to a localized decrease in erosion rates. This reduction in mineral supply effectively decouples the weathering flux from the kinetic acceleration expected from rising temperatures.

We have revised the text to replace "chemostatic" with "kinetically limited" and added a discussion on the coupling between the negative temperature-erosion response and the attenuated weathering feedback.

[Comment 7] Your paper demonstrates that an erosion model based on 7 input parameters performs equally well to a model based on 14 parameters. That makes me wonder: would a model with even fewer parameters also work? For example, if I understand Fig. 6 correctly, the lithology may play a minor role. Would it be feasible to progressively eliminate parameters from your model (starting with the least important) and see when the model starts breaking down?

Response: We thank the reviewer for this insightful suggestion regarding model robustness. Achieving a balance between model simplicity and predictive accuracy is particularly important for deep-time applications where boundary conditions are less certain.

While our SHAP analysis (Fig. 6) suggests that lithology (ROCK) has a lower GLOBAL importance compared to other predictors, it remains a fundamental factor in geomorphology as it dictates rock strength and initial cation concentrations (Bluth and Kump, 1994; Bufe et al., 2022). When modeling LOCAL erosion rates, lithology sometimes serves as the only other predictor other than slope (e.g., Zondervan et al., 2023). Removing this variable could lead to systematic spatial biases, especially in regions with distinct lithological controls.

To address this, we have previously conducted experiments by progressively removing variables. When ROCK was eliminated (leaving SLOPE, ELEVATION, MAT, LAI, RUNOFF, and PWET), the erosion model's validation R^2 decreased to 0.78 (from 0.81 in the baseline), and the weathering model's maximum R^2+R_log^2 dropped to 0.72 (from 0.86 in the 7-parameter version). Although the performance decline is not "catastrophic," the 7-parameter model consistently provides a superior fit to both erosion and weathering datasets. So are the results of experiments eliminating other 6 predictors.

To further evaluate the necessity of these 7 variables, we have performed a Permutation Importance Test. Unlike SHAP, which measures local contributions, this test assesses global importance by randomly shuffling a single feature's values and measuring the resulting change in the model's Mean Absolute Error (MAE).

In this test, all seven variables yielded a positive MAE change, meaning the exclusion of any single variable increases the model's prediction error (Figure R1). In terms of its impact on global error (MAE change), ROCK was found to be more significant than both RUNOFF and LAI. This confirms that while lithology may not be the primary driver of local variation in all cases, it is essential for the model's overall stability and global accuracy. Therefore, we believe the 7-parameter configuration represents the "optimal" set of predictors for MErSiM v1.0.

 

Figure R1 Mean Absolute Error (MAE) change of randomly shuffling a single feature's value

[Comment 8] L10: Which “widely used weathering model”- can you be specific? There are several out there.

Response: We referred to the GM09 framework. The rephrased L10 is “The widely-used weathering model framework of Gabet and Mudd (2009), when driven by stream power erosion laws…”

[Comment 9] L51: became

Response: Fixed to "became." Thanks.

[Comment 10] L60: How can a model “bias [be] supported by […] data”? Wasn’t your point that the bias is a bias because it doesn’t match global fluxes?

Response: Thanks. We have rephrased. The "bias" is supported by data in the sense that independent proxies (like Osmium isotopes) also suggest the model's spatial distribution of weathering is incorrect.  The rephrased L60 is “It stems from a systematic overestimation of weathering fluxes specifically in tropical regions, a spatial mismatch corroborated by osmium isotope data (Rugenstein et al., 2021).”

[Comment 11] L102 – 136: Please systematically define all variables (even if they are obvious); e.g. missing t in equation 1 (definition comes only for equation 2) and many of the parameters in the Arrhenius relationship (equation 6).

Response: We have defined all parameters in the manuscript. The rephrased explanation is “The core of the GM09 model describes the change as a function of time (t, in years) in regolith thickness (h, in meters) as a balance between…” for equation (1) and “Where E_a is the activation energy (J/mol), R is the ideal gas constant (8.314 J/(mol·K)), T_0 is the reference temperature (K), and k_rp is a proportionality constant” for equation (6).

[Comment 12] L140: I guess you could cite West (2012) here who explores that steady-state expression that you show

Response: Thanks. We have added the citation here.

[Comment 13] L149: You have the variable name R already in the Arrhenius relationship. It would be useful to choose a different name for runoff – for example qw.

Response: We have renamed the runoff variable from R to q to avoid confusion with the universal gas constant.

[Comment 14] L263: By “this compilation” you refer to the Gaillardet et al. (1999) compilation, right? Maybe specify to clarify sentence

Response: We will specify this compilation as “the compilation of Gaillardet et al. (1999) and HYBAM”

[Comment 15] Figure 3: Please define the abbreviations in the figure caption

Response: We have added definitions for all abbreviations (e.g., PWET, LAI) and units (e.g., mm/kyr, °C) to the caption. Thanks for bringing this to our attention.

[Comment 16] Figure 8: Can you indicate the units of these variables?

Response: We have added the units in the caption.

Reference

Anderson, R.S., Anderson, S.P., 2010. Geomorphology : The Mechanics and Chemistry of Landscapes.

Bluth, G., Kump, L., 1994. Lithologic and Climatologic Controls of River Chemistry. Geochimica Et Cosmochimica Acta 58, 2341–2359. https://doi.org/10.1016/0016-7037(94)90015-9

Bufe, A., Cook, K.L., Galy, A., Wittmann, H., Hovius, N., 2022. The effect of lithology on the relationship between denudation rate and chemical weathering pathways – evidence from the eastern Tibetan Plateau. Earth Surface Dynamics 10, 513–530. https://doi.org/10.5194/esurf-10-513-2022

Coogan, L.A., Caves Rugenstein, J.K., 2025. Regulation of the carbon cycle on geological timescales, in: Treatise on Geochemistry. Elsevier, pp. 419–465. https://doi.org/10.1016/B978-0-323-99762-1.00060-7

Dessert, C., Dupré, B., Gaillardet, J., François, L.M., Allègre, C.J., 2003. Basalt weathering laws and the impact of basalt weathering on the global carbon cycle. Chemical Geology, Controls on Chemical Weathering 202, 257–273. https://doi.org/10.1016/j.chemgeo.2002.10.001

Gabet, E.J., Mudd, S.M., 2009. A theoretical model coupling chemical weathering rates with denudation rates. Geology 37, 151–154. https://doi.org/10.1130/G25270A.1

Gaillardet, J., Dupré, B., Louvat, P., Allègre, C.J., 1999. Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers. Chemical Geology 159, 3–30. https://doi.org/10.1016/S0009-2541(99)00031-5

Maher, K., Chamberlain, C.P., 2014. Hydrologic Regulation of Chemical Weathering and the Geologic Carbon Cycle. Science 343, 1502–1504. https://doi.org/10.1126/science.1250770

Moon, S., Chamberlain, C.P., Hilley, G.E., 2014. New estimates of silicate weathering rates and their uncertainties in global rivers. Geochimica et Cosmochimica Acta 134, 257–274. https://doi.org/10.1016/j.gca.2014.02.033

Rugenstein, J.K.C., Ibarra, D.E., Zhang, S., Planavsky, N.J., Von Blanckenburg, F., 2021. Isotope mass-balance constraints preclude that mafic weathering drove Neogene cooling. Proc. Natl. Acad. Sci. U.S.A. 118, e2026345118. https://doi.org/10.1073/pnas.2026345118

Zhao, J., Liu, Y., Li, G., Zuo, H., 2026. Thresholds in the controls of denudation rates: A global analysis of tectonic, climatic and biological factors based on machine learning. Earth and Planetary Science Letters 674, 119750. https://doi.org/10.1016/j.epsl.2025.119750

Zondervan, J.R., Hilton, R.G., Dellinger, M., Clubb, F.J., Roylands, T., Ogrič, M., 2023. Rock organic carbon oxidation CO2 release offsets silicate weathering sink. Nature 1–5. https://doi.org/10.1038/s41586-023-06581-9

Response to reviewers’ comments

Response to Reviewer #2:

[Comment 1] Zhao and co-authors develop a new data-driven model of erosion and weathering and compare this model to similar output presented in Park et al. (2020), which uses a stream-power model coupled to Gabet and Mudd (2009; GM09). They find that their new model does a better job capturing global erosion rates, lowers estimated silicate weathering fluxes, and makes silicate weathering much less sensitive to changes in climate.

I found this an interesting paper to read and the figures are well-done. I especially appreciated the explanations of how the model achieved certain results, and the random-forest technique is especially well-suited to understand the controls on erosion. I also liked the sensitivity test of weathering to changes in CO2, since this is an important parameter in Earth system science that is poorly constrained. However, I have a number of concerns that should be addressed, so that readers understand the foundation on which this model was built.

Response: We sincerely thank Dr. Rugenstein for his constructive and encouraging assessment of our work. We are gratified by his positive evaluation of our data-driven modeling approach and the clarity of our figures and explanations. We have carefully considered and addressed each of the concerns raised to strengthen the foundation and clarity of the MErSiM framework, as detailed in our point-by-point responses below.

[Comment 2] First, the paper refers repeatedly to the “original GM09” model as including a stream-power component. GM09 does not include anything remotely related to stream-power and is perhaps better described as a reactive-transport-inspired model that attempts to predict weathering on hillslopes. I don’t know which version of GEOCLIM included a stream-power model piggy-backed on to GM09 (perhaps an earlier paper by Pierre Maffre? (Maffre et al., 2018)), but this is used in Park et al. (2020), though I would definitely not characterize this as the original GM09 model.

Response: We sincerely thank the reviewer for this important correction and clarification. We agree that the original Gabet and Mudd (2009) (GM09) model is a reactive-transport-inspired framework specifically designed to predict regolith production and weathering on hillslopes, and it does not include a stream-power erosion component.

In our manuscript, we used "original GM09" as a shorthand to refer to the specific global modeling configuration implemented by Park et al. (2020), which coupled the GM09 weathering physics with a Stream Power Incision Model (SPIM) to represent the erosion driver. We recognize that this terminology is imprecise and potentially misleading to readers familiar with the original GM09 study.

To resolve this, we have revised the manuscript to distinguish between the core weathering physics (GM09) and the erosion driver (SPIM). We have added a paragraph to the Introduction section describing the combination of GM09 and SPIM in GEOCLIM models: “The GM09 model has thus became a popular tool, notably being integrated into larger Earth System Models like GEOCLIM to explore long-term climate-tectonic interactions, where the erosion rate is provided by a Stream Power Incision Model (SPIM) (Maffre et al., 2018; Park et al., 2020). We refer to this combined configuration as GM09-SPIM.”

We have also replaced all subsequent references to “original GM09” with the more accurate term “GM09-SPIM” throughout the text.

[Comment 3] Second, the authors confuse the meaning of chemostatis (Godsey et al., 2009). Chemostasis refers to constant concentrations as runoff changes. This, in turn, means the silicate weathering fluxes scale linearly with runoff. I believe the authors mean some sort of kinetic-limit, such that, as runoff increases, concentrations decline, and silicate weathering fluxes are—as a consequence—flat.

Response: We sincerely thank the reviewer for this crucial clarification regarding the term "chemostasis." We agree with the reviewer’s definition based on Godsey et al. (2009, 2019): chemostasis implies that solute concentrations remain nearly constant as runoff changes.

The reviewer is correct in intuiting our intended meaning. In the original manuscript, we used "chemostatic" to describe a state where the weathering flux reaches a plateau and becomes insensitive to further increases in runoff. As the reviewer suggested, this is indeed a manifestation of a kinetic limit or regime. Our model results, consistent with the framework of Maher and Chamberlain (2014), indicate that under high-runoff conditions, the system is pushed toward this kinetic boundary because fluid travel times become too short for reactive equilibration.

We recognize that our previous terminology was imprecise and have revised the manuscript to replace "chemostatic" with "kinetically limited" throughout the text.

[Comment 4] More critically, the authors motivate much of the paper by arguing that Park et al. (2020) overestimate the silicate weathering flux. Maybe they do, maybe they don’t, but the model presented herein (MErSiM) produces a very low silicate weathering flux. It may match data from Müller et al. (2022), but the uncertainty on the volcanic flux is enormous (over an order of magnitude) (Coogan and Rugenstein, 2025) (see their Figure 13). To the extent that basically the silicate weathering flux must match the volcanic flux, I don’t think optimizing a model to the low-end of this estimate is particularly useful. What would happen if the model was optimized for a higher-end estimate? A higher-end estimate is, anyway, what Moon et al. (2014) estimate; the authors state that their model matches Moon’s estimates, though Moon estimates a global silicate weathering carbon flux nearly an order of magnitude higher than MErSiM. Either way, the authors need to address what would happen if the model was optimized to a different global silicate weathering flux.

Response: We thank the reviewer for the challenging comment regarding the global flux magnitude. We particularly appreciate the opportunity to clarify the geochemical benchmarks used in our study.

We have carefully re-examined the results in Moon et al. (2014) and believe there may be a misunderstanding regarding the specific flux metrics reported. Moon et al. (2014) as well as Gaillardet et al. (1999) distinguished between the Ca+Mg silicate weathering flux and the total CO2 consumption rate. In Moon et al. (2014), the global SWRCa+Mg—which is the appropriate benchmark for long-term alkalinity and carbon balance models like ours—is estimated at 2.17 × 10¹² mol/yr (2σ range: 1.59–2.75 × 10¹² mol/yr). Our model’s target 2.5 × 10¹² mol/yr is from the estimation of Gaillardet et al. (1999) and is aligned with this "Ca+Mg" benchmark SWRCa+Mg.

The higher figure mentioned by the reviewer (~11 × 10¹² mol/yr) likely refers to Moon et al.'s estimate of the total global CO2 consumption rate, where 7.85 × 10¹² mol/yr is from silicates (including SWRCa+Mg and Na/K silicates), and 4.08× 10¹² mol/yr is from volcanic provinces estimated by Dessert et al. (2003). 

In the context of long-term carbon cycle modeling (such as the GEOCLIM/MErSiM framework), the Ca+Mg silicate weathering flux is the critical parameter because only the weathering of Ca and Mg silicates provides the alkalinity required for marine carbonate precipitation and long-term CO2 sequestration. Weathering of Na and K silicates consumes CO2 in the short term but does not lead to carbonate burial. Our model’s predicted magnitude of 3.11 × 10¹² mol/yr is, therefore, not a "low-end" estimate but is aligned with—and even slightly higher than—the specific Ca+Mg benchmark (2,17 × 10¹² mol/yr) reported by Moon et al. (2014).

Regardless of the chosen global target, the primary contribution of MErSiM is its improved spatial consistency. As shown in our sensitivity analysis (Fig. 9), the MErSiM framework achieves optimal spatial performance with high R^2+R_log^2 for 81 watersheds around the world. This emergent property suggests that previous models' overestimations were indeed a symptom of incorrect spatial erosion distribution, rather than just an adjustable global bias. We have revised Section 3.5.1 to explicitly clarify these geochemical definitions and ensure the comparison with Moon et al. (2014) is transparent.

We agree that the uncertainty in volcanic degassing remains high as demonstrated in (Coogan and Caves Rugenstein, 2025). We believe that further work is needed to estimate weathering fluxes that include Na, K, and volcanic provinces, and to compare these with new estimates of degassing fluxes.

[Comment 5] Relatedly, I don’t follow the argument about why erosion is the particular parameter that is most underconstrained and therefore produces a wrong silicate weathering flux. Maybe I just don’t understand the argument, but the fact that we don’t even know what the silicate weathering flux should be means it’s difficult to assign to a specific parameter why a certain model doesn’t match certain estimates of silicate weathering. I like the authors’ approach to building a better erosion model; however, I don’t see this statement as a particularly helpful motivation for building such a model.

Response: We thank the reviewer for this critical reflection on our model's motivation. We agree that there are too many possible causes for an off estimate of silicate weathering, erosion rate does not have to be the one. Our motivation was from the previous attempt (Zuo et al., 2024) to reduce the systematic bias in the tropical region in the GM09-SPIM model. In that paper, we exhausted ways of reducing the bias by using different climate forcings and fitting parameters but failed all of them, even when chemical kinetic parameter (k_d) was tuned to their lower limits. The only way we could find was to substantially reduce the erosion rate within the tropical region, which motivated us to build a better erosion model in Zhao et al. (2026). However, before doing the test in the current work, we do not know whether such a data-driven erosion model would also improve the silicate weathering flux since we did not tune the model specifically to reduce the erosion rate within the tropical region; the difference in erosion rates between this data-driven model and the subjectively tuned model of Zuo et al. (2024) is large (as shown in Fig. 5e of this manuscript). The fact that the MErSiM produces better results than the GM09-SPIM model then confirms that improving the erosion model is critical to improving the silicate weathering model.

Regarding the reviewer's statement "the fact that we don’t even know what the silicate weathering flux should be", we agree that the inherent uncertainties in weathering fluxes of Na+K silicate and volcanic provinces make the global CO2 consumption difficult to estimate. Before any obviously improved estimate becomes available, we think it is advisable to choose to trust the results of Gaillardet et al. (1999) and Moon et al. (2014), who gave similar estimation for Ca+Mg silicate weathering flux, with a discrepancy of less than 12%.

Importantly, other than the global total, spatial pattern of weathering fluxes is also improved by MErSiM. Specifically, MErSiM reduced the Mean Absolute Error (MAE) across all 81 basins by approximately 35% compared to the original GM09-SPIM (1.77 × 10^10  mol 〖yr〗^(-1) compared to 2.72 × 10^10  mol 〖yr〗^(-1)). Second, the number of basins exhibiting "large" positive errors (5 × 10^10  mol 〖yr〗^(-1)) dropped from 10 in GM09-SPIM to 4 in MErSiM. Third, in the tropics (basins located between 23S-23N), where GM09-SPIM errors often exceeded 7 × 10^10  mol 〖yr〗^(-1), MErSiM reduced the average basin-scale bias by 90% (0.41 × 10^10  mol 〖yr〗^(-1) compared to 4.12 × 10^10  mol 〖yr〗^(-1)). Therefore, we are not simply adjusting a global scalar to match an uncertain target, but rather resolving a structural discrepancy in how mineral supply is distributed globally. 

[Comment 6] A few other caveats on the erosion model would be valuable. For example, I presume that there is no glacial erosion processes, meaning that calculating an LGM erosion flux assumes that glaciers only negligibly modify erosion and weathering processes. This is, of course, a major controversy in the field, and the lack of such a glacial erosion process is a major caveat. Similarly, at high CO2 levels, one robust expectation of a warming climate is a change in the timing and intensity of storms, which is likely to modify erosion. Based upon the parameters in MErSiM, I don’t think such a shift in precipitation timing/intensity and its effect on erosion is captured in MErSiM. Another caveat that should be mentioned.

Response:  We sincerely thank the reviewer for identifying these two critical caveats. We agree that the absence of explicit glacial processes and high-frequency intense precipitation represents a significant boundary for our model's interpretations, especially for the LGM and 4×CO₂ scenarios.

The reviewer is correct that MErSiM is primarily trained on modern cosmogenic-nuclide derived erosion rates. Since these training basins are located in non-glaciated or minimally glaciated regions, the model likely represents an "ice-free" erosional response to climate. We acknowledge that the lack of glacial grinding—which can significantly increase reactive surface area and modulate weathering—is a major caveat for our LGM simulations. We have added a dedicated paragraph to Section 3.6 (Scope of Applicability and Model Limitations) to address this. 

We also agree that MErSiM's use of monthly/annual climatic averages (e.g., MAT, PWET, and Runoff) does not resolve changes in the timing and intensity of individual storms. While runoff acts as a macro-scale proxy for water flux, it cannot capture the erosional pulses driven by "flashier" precipitation regimes expected in a warmer world. The problem is not only in the erosion model, but also in the climate models where the storms are difficult to be simulated due to the low spatial resolution. One well-known example is the simulation of tropical cyclones, which have strong influence on the precipitation in coastal regions but require horizontal resolution of at least 0.25°×0.25° (Chang et al., 2020), far higher than normally used (coarser than 2°×2°) for the deep paleoclimate (e.g., Goddéris et al., 2012; Li et al., 2022; Valdes et al., 2021). This limitation is now explicitly discussed as a caveat in the context of our Section 3.6. We believe these additions significantly improve the transparency and foundational understanding of the model's scope.

[Comment 7] I found the sensitivity experiment to be particularly interesting, and I appreciated the explanation of why the model responds the way it does. However, it should be noted that in, for example, the 4x CO2 experiment, the rise in silicate weathering flux is best understand as the instantaneous transient adjustment to a change in climate. Given the need to maintain mass balance in the carbon cycle (Berner and Caldeira, 1997; Caves et al., 2016; Zeebe and Caldeira, 2008), the silicate weathering flux will ultimately have to rise by exactly the same amount as whatever input flux of CO2 caused the increase in atmospheric CO2. This would presumably involve changes in erosion beyond what is predicted by MErSiM (or including other components of the Earth system, such as weathering on marine shelves (Trapp-Müller et al., 2025) or in seafloor basalts (Coogan and Gillis, 2018)). However, the authors state that MErSiM is best thought of as a long-term (and not short-term, transient) model. Some acknowledgement that this is still a barrier to be resolved would help to place these sensitivity experiments in context.

Response: We sincerely thank the reviewer for this insightful comment regarding the fundamental requirement of carbon cycle mass balance (Berner and Caldeira, 1997). We agree that in a long-term steady state, the global carbon sink must equal the CO₂ input. 

In the revised manuscript, we clarify that our 4×CO₂ simulations using CMIP6 datasets are designed as sensitivity experiments rather than transient simulations to a new steady state. They represent an instantaneous response of the current terrestrial landscape to a prescribed climate forcing. The result reveals that the terrestrial Ca+Mg silicate weathering thermostat in MErSiM is less efficient under extreme greenhouse conditions than previously models. However, we acknowledge that components of the Earth system may provide more rapid transient feedbacks, including:

1) as shown by Li et al. (2016), basalt weathering exhibits a much stronger temperature dependence in natural settings compared to global averages. These reactive terrains (e.g., volcanic islands and large igneous provinces) can deliver a rapid surge in weathering flux during the early stages of a greenhouse warming event, effectively acting as a high-frequency regulator before the more extensive continental surfaces reach equilibrium.

2) while long-term carbonate-silicate cycle models focus on Ca+Mg flux, the weathering of alkali silicates (Na+K) provides a significant transient CO₂ sink that responds rapidly to changes in runoff and temperature. This transient drawdown can be crucial in damping climate excursions on shorter timescales, even if it does not lead to permanent carbon burial.

3) as the reviewer correctly noted, other important sinks such as marine shelf weathering (Trapp-Müller et al., 2025) and seafloor basalt alteration (Coogan and Gillis, 2018) may also adjust to fulfill the mass balance requirement when the terrestrial continental sink approaches kinetic or supply limits. 

We have revised the manuscript Section 3.5.2 to acknowledge that the weak sensitivity of the "continental thermostat" in MErSiM implies a greater reliance on these highly reactive lithologies and marine sinks to maintain global mass balance. 

[Comment 8] Lastly, equation 10 (similar to West (2012)) demonstrates that erosion will be the predominant variable impacting the estimate silicate weathering flux (Equation 10 needs to be derived…it’s not clear to me how the authors reach Equation 10 from the previous equations). This is in contrast to Maher and Chamberlain (2014) who parameterize silicate weathering as being predominantly impacted by runoff and the equilibrium concentration (equation 3 and Fwsil = q*C, where q is runoff and C is concentration). I bring this up to point out that I don’t think this paper demonstrates that their formulation is indeed the best formulation to understand silicate weathering. It simply does a better job than a similar formulation that also places a heavy emphasis on erosion. One would need to do a similar analysis using Maher and Chamberlain (2014) (or another model) to demonstrate which model produces a better estimate of the silicate weathering flux. Even this is difficult, since our constraints on the weathering flux are poor. Thus, MErSiM ends up predicting that catchments are near the kinetic boundary and will be sensitive to changes in erosion, but it’s not clear to me that this is supported by the data, particularly if another model is used to interpret catchment solute data. Perhaps this is beyond the scope of this paper, but it is an overall caveat in how one conceptualizes the modeling presented herein.

Again, thanks for an interesting paper, and I hope these comments are helpful in revising the paper and addressing any outstanding questions.

Response: We thank the reviewer for pointing out that the derivation of Equation 10 was not sufficiently transparent. We have added a detailed derivation in the manuscript (also in the response below of [Comment 13]). 

We appreciate the reviewer’s insightful comparison with the framework of Maher and Chamberlain (2014). We agree that their model conceptualizes weathering through a different lens. The fundamental difference between two branches of models lies in the leading variable—whereas Maher and Chamberlain (2014) focus on the competition between fluid travel time and equilibrium time, MErSiM (rooted in GM09/West 2012) focuses on the competition between mineral supply (erosion) and reaction kinetics within a regolith column.

However, the two frameworks are not mutually exclusive. Maher and Chamberlain (2014) also acknowledge that the mineral supply rate affects the Damköhler coefficient and thus the sensitivity of the system.

The reviewer suggests that our prediction of catchments being near the "kinetic boundary" (where fluxes are sensitive to erosion) might be a consequence of our model's structure rather than the data itself. We emphasize that this conclusion is an emergent property of our objective global optimization. In our calibration across 81 major basins, the model was free to select parameters that could have placed these basins in a regime more indifferent to erosion rate, that is, a large enough base dissolution rate constant k_d. But when MErSiM achieved a spatial and global flux estimation more consistent with data, it used k_d=1×10^(-4). Therefore, during this parameter calibration process, data “chose” a model more sensitive to erosion.

We accept the reviewer's point that we have not definitively proven MErSiM to be the "best" possible formulation for silicate weathering. Our study demonstrates that MErSiM is an advancement over previous erosion-driven frameworks (like Park et al. 2020) because it resolves longstanding spatial biases. We have tempered our language and added a paragraph to Section 3.6 to acknowledge that different conceptual frameworks may lead to different interpretations of catchment data. We agree that a formal inter-model comparison using the same global datasets would be a valuable direction for future research.

[Comment 9] Line 18: What is the systematic tropical overestimation?

Response: By "systematic tropical overestimation," we refer to the phenomenon where GM09-SPIM predict weathering flux in tropical basins (e.g., the Amazon and Congo) that are significantly higher than observed values of Gaillardet et al. (1999).

[Comment 10] Line 48: I’m not sure the global degassing flux is that well-constrained. Müller et al. (2022) may claim it is, but recent compilations suggest nearly an order of magnitude uncertainty (Coogan and Rugenstein, 2025) (see their Figure 13 and references therein).

Reponse: We thank the reviewer for pointing out this recent and important work. We have added a citation to Coogan and Rugenstein (2025) and revised the text to acknowledge that the global degassing flux is subject to nearly an order of magnitude of uncertainty, rather than being strictly constrained by any single study.

[Comment 11] Line 67: I would not characterize GM09 as using a simplified stream power incision model. If anything, they use a sort of reactive-transport-inspired model that looks only at hillslopes.

Response: We have revised the manuscript to clarify that the original GM09 (Gabet and Mudd, 2009) does not include a stream-power component.

[Comment 12] Equation 8: This equation is not in the original GM09 model. I believe it is a modification in GEOCLIM. In the original GM09 model, erosion is simply an independent variable.

Response: We have revised the manuscript to clarify that this equation is from GM09-SPIM framework.

[Comment 13] Equation 10: Please provide the derivation for this equation.

Response: Thanks. We have added detailed derivation for Equation 10 in the manuscript.

“While GM09 is a transient model, for long-term geological applications it is typically run to a steady-state solution where the regolith thickness, mineral concentration and exposure time are all constant (dh/dt= 0, ∂x/∂t= 0, ∂τ/∂t= 0). Under this assumption, the regolith production rate equals the erosion rate (P_r=E) according to Eq. 1. The particle residence time at height z above the bedrock therefore simplifies to τ =z/E according to Eq.3. Substituting these into the mineral mass balance (Eq. 2) yields K(z/E)^σ x=-E (∂x/∂z). Integrating this rate throughout the entire regolith column (z=0 to z=h, Eq. 4) yields the total flux W = E(├ x┤|_(z=0)-├ x┤|_(z=h) ), where ├ x┤|_(z=0) and ├ x┤|_(z=h) are the mineral concentrations at the bedrock and surface, respectively. Vertical integration of K(z/E)^σ x=-E (∂x/∂z) gives ∫_(├ x┤|_(z=0))^(├ x┤|_(z=h))▒〖1/x dx=-K/E^(σ+1) 〗 ∫_0^h▒〖z^σ dz 〗, from which we have the solution ├ x┤|_(z=h)=├ x┤|_(z=0) e^((-K)/(σ+1) (h/E)^(σ+1) ). Substituting the expression of K (Eq. 9), the weathering flux W simplifies to an expression:

█(W =E(├ x┤|_(z=0)-├ x┤|_(z=0) e^((-k_d (1 - e^(-k_w  q) )  e^[(E_a/R)  (1/T_0 -1/T)] )/(σ+1) (h/E)^(σ+1) ) )  #(10) )”

[Comment 14] Lines 185-9: Only 2 purposes are listed, rather than 3.

Response: Thanks. The manuscript has been corrected.

[Comment 15] Line 490: I wouldn’t say that the peak of the yellow curve is strictly within the gray bar…it is still somewhat shifted to the right.

Response: The reviewer is correct that the peak of the yellow curve (R_log^2) in Figure 9b sits at ~3.5 × 10^12  mol 〖yr〗^(-1), which is indeed outside the 2.0-3.0 × 10^12  mol 〖yr〗^(-1)  grey band.

We have revised the text to clarify that the peak of the joint metric (R^2+R_log^2, blue curve) and the standard R^2 (red curve) now fall within the observational range instead of all three metrics. 

[Comment 16] Line 515: The Moon et al. (2014) estimate is substantially higher when considering all fluxes (ie, ~11 x 10^13 mols C/yr)

Response: We thank the reviewer for prompting us to clarify our comparison with the estimates in Moon et al. (2014). The reviewer correctly points out that Moon et al. (2014) report a higher flux when considering total CO2 consumption. However, the global Ca + Mg silicate weathering flux in that paper is 2.17 × 10¹² mol/yr. In long-term carbon cycle modeling (including MErSiM and the GEOCLIM framework), this Ca+Mg flux is the exact target metric, as it represents the alkalinity flux that ultimately leads to long-term carbon sequestration via marine carbonate precipitation (excluding Na and K weathering). Our target flux of ~2.5 × 10¹² mol C/yr aligns with this specific metric.

The higher estimate mentioned by the reviewer corresponds to Moon et al.'s "global CO2 consumption rates from silicates" when including volcanic provinces, which they report as 11.93 × 10¹² mol/yr (we note that the "11 × 10¹³" in the comment appears to be a minor typographical error in the exponent). Because MErSiM specifically simulates the long-term stabilizing carbon flux (Ca+Mg), it is geochemically consistent to benchmark against their 2.17 × 10¹² mol/yr estimate.

To avoid any confusion for future readers, we have revised the manuscript to explicitly clarify that our benchmark corresponds specifically to the Ca+Mg silicate weathering flux reported by Moon et al. (2014), rather than the total CO2 consumption rate.

[Comment 17] Line 567: “Chemostatic” has the opposite meaning; that is, solute concentrations remain the same as runoff changes, which has been used to suggest that the catchment has reached reaction equilibrium. I think you mean to say, “kinetically limited”.

Line 570: Again, I think you mean close to the kinetic boundary

Line 618: Again, wrong use of chemostatic

Response: We are very grateful to the reviewer for pointing out the incorrect use of the term "chemostasis." We agree that in the hydrological community (e.g., Godsey et al., 2009, 2019), chemostasis typically refers to a state where solute concentrations remain nearly constant despite changes in discharge.

In the original manuscript, we used "chemostatic" in the context of the Damköhler number framework (Maher and Chamberlain, 2014) to describe a state of weathering flux saturation. We acknowledge that this usage is ambiguous and contradicts the definition favored by the reviewer.

We have revised the text to replace "chemostatic" with "kinetically limited" and added a discussion on the coupling between the negative temperature-erosion response and the attenuated weathering feedback.

Reference:

Berner, R.A., Caldeira, K., 1997. The need for mass balance and feedback in the geochemical carbon cycle. Geology 25, 955–956. https://doi.org/10.1130/0091-7613(1997)025%253C0955:TNFMBA%253E2.3.CO;2

Chang, P., Zhang, S., Danabasoglu, G., Yeager, S.G., Fu, H., Wang, H., Castruccio, F.S., Chen, Y., Edwards, J., Fu, D., Jia, Y., Laurindo, L.C., Liu, X., Rosenbloom, N., Small, R.J., Xu, G., Zeng, Y., Zhang, Q., Bacmeister, J., Bailey, D.A., Duan, X., DuVivier, A.K., Li, D., Li, Y., Neale, R., Stössel, A., Wang, L., Zhuang, Y., Baker, A., Bates, S., Dennis, J., Diao, X., Gan, B., Gopal, A., Jia, D., Jing, Z., Ma, X., Saravanan, R., Strand, W.G., Tao, J., Yang, H., Wang, X., Wei, Z., Wu, L., 2020. An Unprecedented Set of High‐Resolution Earth System Simulations for Understanding Multiscale Interactions in Climate Variability and Change. J Adv Model Earth Syst 12, e2020MS002298. https://doi.org/10.1029/2020MS002298

Coogan, L.A., Caves Rugenstein, J.K., 2025. Regulation of the carbon cycle on geological timescales, in: Treatise on Geochemistry. Elsevier, pp. 419–465. https://doi.org/10.1016/B978-0-323-99762-1.00060-7

Coogan, L.A., Gillis, K.M., 2018. Low-Temperature Alteration of the Seafloor: Impacts on Ocean Chemistry. Annu. Rev. Earth Planet. Sci. 46, 21–45. https://doi.org/10.1146/annurev-earth-082517-010027

Dessert, C., Dupré, B., Gaillardet, J., François, L.M., Allègre, C.J., 2003. Basalt weathering laws and the impact of basalt weathering on the global carbon cycle. Chemical Geology, Controls on Chemical Weathering 202, 257–273. https://doi.org/10.1016/j.chemgeo.2002.10.001

Gabet, E.J., Mudd, S.M., 2009. A theoretical model coupling chemical weathering rates with denudation rates. Geology 37, 151–154. https://doi.org/10.1130/G25270A.1

Gaillardet, J., Dupré, B., Louvat, P., Allègre, C.J., 1999. Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers. Chemical Geology 159, 3–30. https://doi.org/10.1016/S0009-2541(99)00031-5

Goddéris, Y., Donnadieu, Y., Lefebvre, V., Le Hir, G., Nardin, E., 2012. Tectonic control of continental weathering, atmospheric CO2, and climate over Phanerozoic times. Comptes Rendus. Géoscience 344, 652–662. https://doi.org/10.1016/j.crte.2012.08.009

Li, Gaojun, Hartmann, J., Derry, L.A., West, A.J., You, C.-F., Long, X., Zhan, T., Li, L., Li, Gen, Qiu, W., Li, T., Liu, L., Chen, Y., Ji, J., Zhao, L., Chen, J., 2016. Temperature dependence of basalt weathering. Earth and Planetary Science Letters 443, 59–69. https://doi.org/10.1016/j.epsl.2016.03.015

Li, X., Hu, Y., Guo, J., Lan, J., Lin, Q., Bao, X., Yuan, S., Wei, M., Li, Z., Man, K., Yin, Z., Han, J., Zhang, J., Zhu, C., Zhao, Z., Liu, Y., Yang, J., Nie, J., 2022. A high-resolution climate simulation dataset for the past 540 million years. Sci Data 9, 371. https://doi.org/10.1038/s41597-022-01490-4

Maher, K., Chamberlain, C.P., 2014. Hydrologic Regulation of Chemical Weathering and the Geologic Carbon Cycle. Science 343, 1502–1504. https://doi.org/10.1126/science.1250770

Moon, S., Chamberlain, C.P., Hilley, G.E., 2014. New estimates of silicate weathering rates and their uncertainties in global rivers. Geochimica et Cosmochimica Acta 134, 257–274. https://doi.org/10.1016/j.gca.2014.02.033

Trapp-Müller, G., Caves Rugenstein, J., Conley, D.J., Geilert, S., Hagens, M., Hong, W.-L., Jeandel, C., Longman, J., Mason, P.R.D., Middelburg, J.J., Milliken, K.L., Navarre-Sitchler, A., Planavsky, N.J., Reichart, G.-J., Slomp, C.P., Sluijs, A., Van Hinsbergen, D.J.J., Zhang, X.Y., 2025. Earth’s silicate weathering continuum. Nat. Geosci. 18, 691–701. https://doi.org/10.1038/s41561-025-01743-y

Valdes, P.J., Scotese, C.R., Lunt, D.J., 2021. Deep ocean temperatures through time. Clim. Past 17, 1483–1506. https://doi.org/10.5194/cp-17-1483-2021

Zhao, J., Liu, Y., Li, G., Zuo, H., 2026. Thresholds in the controls of denudation rates: A global analysis of tectonic, climatic and biological factors based on machine learning. Earth and Planetary Science Letters 674, 119750. https://doi.org/10.1016/j.epsl.2025.119750

Zuo, H., Liu, Y., Li, G., Xu, Z., Zhao, L., Guo, Z., Hu, Y., 2024. A revised model of global silicate weathering considering the influence of vegetation cover on erosion rate. Geosci. Model Dev. 17, 3949–3974. https://doi.org/10.5194/gmd-17-3949-2024