A greening Earth has reversed the trend of decreasing carbonate weathering under a warming climate

Carbonate weathering simulations

The research area in this study is located in Southwest China (102°25′–114°25′E, 21°59′–30°25′N), including eight provinces (Fig. 1). Southwest China is noted for its extensive carbonate outcrops (total of ~9.7 × 10⁵km², mainly limestone) and diverse karst landscapes³². This area is dominated by a subtropical monsoon climate and rainfall events are concentrated from April to September⁴⁹. In the last century, deforestation, and damaging cultivation practices caused widespread rock desertification and ecosystem degradation in the carbonate areas of Southwest China⁵⁰. Hence, China has invested nearly 19 billion USD to implement several national ERPs in this region since the 1980s³⁴. These ERPs have resulted in a substantial increase in vegetation cover and have led to Southwest China becoming one of the most striking greening regions worldwide over the past few decades^33,35 (Fig. 1a, b). In this study, we focused on the carbonate areas in Southwest China and simulated the variations in [HCO₃⁻] and WS_atm-CO2 during this greening period. Additionally, we simulated [HCO₃⁻] in the global carbonate areas during the same period. Details of the biogeochemical models (NPP–pCO₂ and MPD), calculations, and the multi-source database used in the simulations are given below.

Carbonate weathering intensity estimation

In previous studies, we applied an MPD to calculate the equilibrium concentration of [HCO₃⁻] as a function of temperature and soil pCO₂¹⁰. In this study we used the results of this carbonate equilibrium model (simulated [HCO₃⁻]) to represent the actual [HCO₃⁻] at the global scale. This model can also separately estimate the individual influences from thermodynamics (temperature) and soil pCO₂ variations. Equation 1 is a simplified equation to represent how carbonate weathering captures atmospheric CO₂. Actually, for a pure calcite dissolution system, the weathering of carbonate is controlled by a series of chemical equations^2,3,13:

These equations are governed by different mass action laws, which are controlled by the temperature-dependent reaction constants K₁, K_C, K_H, and K₂, and the temperature, T (°C)¹³. K₁ and K₂ are the first and second dissociation constants for carbonic acid in water (Eq. 4 and Eq. 3), and K_c is the equilibrium constant for calcite dissolution (Eq. 2). K_H is the Henry’s Law constant for CO₂ gas in water (Eq. 5)¹⁴. These reactions are characterized by rapid kinetics and the drainage water can reach an equilibrium state within a couple of hours¹³. The pCO₂ (atm) at the water–rock interface (atmosphere or soil) is an additional determinant of the magnitude of HCO₃⁻ production¹³.

In nature, the water interacting with CaCO₃ is impure and the activity coefficient should be introduced to relate concentrations to activities. We can denote \({\gamma }_{1}\) and \({\gamma }_{2}\) as the activity coefficients for single-charged and double-charged species, respectively, and the Davies equation can then be applied to calculate the activity coefficients³:

Here, A is a constant based on temperature, z is the charge of the chemical species, and I is the ionic strength of the solution. Assuming that in a natural water body (pH ~8), the solutes include several common single-charged and double-charged species (Ca²⁺, Mg²⁺, Na⁺, K⁺, HCO₃⁻, CO₃²⁻, SO₄²⁻, NO₃⁻, and Cl⁻), then the electro-neutrality of its composition implies that:

We can then group the species with no acid-base behavior and neglect H⁺, OH⁻ and CO₃²⁻, and Eq. 7 can be modified to:

with \({RA}=2[{{{{\rm{SO}}}}}_{4}^{2-}]+[{{{{\rm{NO}}}}}_{3}^{-}]+[{{{{\rm{Cl}}}}}^{-}]-2[{{{{\rm{Mg}}}}}^{2+}]-[{{{{\rm{Na}}}}}^{+}]-[{{{{\rm{K}}}}}^{+}]\) being the reduced alkalinity³. In this study, we mainly consider the influences of greening and warming on soil pCO₂ and the thermodynamics (temperature-dependent reaction constants). Thus, an ideal pure calcite dissolution–controlled system is assumed in global carbonate areas. Under such conditions RA can set to 0 at equilibrium (the solutions only contain the ionic species Ca²⁺, HCO₃⁻, CO₃²⁻, H⁺, and OH⁻), and then the equilibrium HCO₃⁻ concentration ([HCO₃⁻]_eq) can be calculated by refs. ^3,13,36:

Here, all the K constants can be estimated by temperature. pCO₂ (atm) is the partial pressure of CO₂ at the water–rock interface, which can be estimated by soil pCO₂ models. γCa²⁺ and γHCO₃⁻ are the activity coefficients for Ca²⁺ and HCO₃⁻ at equilibrium (Ca²⁺ and HCO₃⁻ are assumed to be the only single-charged and double-charged species in the pure calcite dissolution system), which can be calculated using Eq. 6. T is the temperature in °C. To save computing time for long-term and large-scale spatial evaluations, we applied an approximate equation which summarizes all the equations (Eqs. 2–8) for ideal calcite dissolution at equilibrium:

This equation requires only temperature and pCO₂ and is valid for pCO₂ > 3 × 10^−4[ 51.

Soil pCO₂ parameterizations

Soil pCO₂ is a critical parameter for [HCO₃⁻] simulation. In previous studies, soil pCO₂ was estimated by different ecological models based on empirical relationships and processes. T, ET, NPP, and SM are the common environmental parameters applied in these models for soil pCO₂ simulation. These parameters are also controlled by vegetation cover variations. In this study, we selected a net primary productivity-based soil pCO₂ models for estimating long-term carbonate weathering intensity. This model can generate accurate simulations across different land-use types at a global scale. Details of these pCO₂ models are given below³:

NPP–based soil pCO₂ (NPP–pCO₂) is a semi-mechanistic model which calculates the pCO₂ profile across soil layers. In this model, CO₂ production in the root zone is assumed to be 75% of NPP ⁵². This ratio is based on a large-scale statistical analysis of in-situ aboveground and belowground carbon measurements for 79 global sites⁵³. Soil pCO₂ can then be expressed as a function of the atmospheric CO₂ concentration, temperature, and the NPP of the vegetation ecosystem³:

where A is the set of conversion unit constants. In this study, the A values for annual and seasonal pCO₂ were set to 1.03 × 10⁶ and 1.236 × 10⁷, pCO_2atm is the atmospheric CO₂ pressure (ppmv) NPP in grams of dry matter per m² for a given year or month (g m⁻² yr⁻¹/month⁻¹), T is the temperature (K), and pCO_2(soil) is the maximum CO₂ pressure reached beneath the root zone (at the soil–rock interface) (ppmv). In previous studies, we compared the results of the NPP–pCO₂ model with field measurements from different global individual studies^5,10. We found that the model produced very accurate results and was widely applicable in different land cover types. Hence, we used NPP–pCO₂ to simulate seasonal [HCO₃⁻] variations in Southwest China and for global carbonate areas. Additionally, to evaluate the influences of vegetation production and warming on future global [HCO₃⁻], we constructed an inverse modelling framework based on NPP–pCO₂ and NCEAS model (NPP = 2540/(1+exp(1.584–0.0622 T)))³⁶. We applied the future temperature data from six CMIP6 models to estimate the required NPP, which can sustain an undiminished global carbonate weathering intensity.

Weathering related atmospheric CO₂ sink calculation

To calculate the total atmospheric CO₂ capture produced by minerals weathering, we used a weathering-related atmospheric CO₂ sink equation (WS_atm-CO2) [t C km⁻² yr⁻¹]²:

Here, WS_atm-CO2 (t C km⁻² yr⁻¹) is the weathering sink for atmospheric CO₂. RF (mm) is the RF, [HCO₃⁻] (mmol L⁻¹) is the bicarbonate concentration (simulated by an NPP–pCO₂ and MPD) driven by carbonates and silicates weathering, and m_r (g mol⁻¹) is the molar mass of carbon. \({f}_{{{{\rm{c}}}}}\) and \({f}_{{{{\rm{s}}}}}\) are the fractions of HCO₃⁻ produced by carbonate and silicate weathering, respectively. The partitioning factor n is 1 for silicate weathering, and 0.5 for carbonate weathering, as in the latter case one part of [HCO₃⁻] is from atmospheric CO₂ and the other part is from carbonate rock. In this study, we only calculate the WS_atm-CO2 driven by carbonate weathering.

Model sensitivity testing and trend uncertainty evaluation

The sensitivity testing of input parameters to models and the trend uncertainty evaluation are significant to validate the model simulation results. In this study, we applied a Bayesian statistical model⁵⁴ by using the PyMC3 package (version 3.10) in a Python 3 environment to evaluate the sensitivity of input parameters on [HCO₃⁻] and WS_atm-CO2 simulations. PyMC3 is an open-source probabilistic programming framework with an intuitive and readable, yet powerful, syntax that is close to the natural syntax statisticians use to describe models. It features next-generation Markov chain Monte Carlo sampling algorithms. We used PyMC3 to construct a Bayesian multivariate linear regression and posterior analysis by using modeling [HCO₃⁻] results (or WS_atm-CO2) and its drivers (model input parameters). Then, the sensitivity of input parameters to the model outcomes (their contributions and significances) was quantified. In addition, to further evaluate the significance of long-term trends in simulated [HCO₃⁻] (or WS_atm-CO2) and its drivers, we also conducted trend uncertainty tests using Chi-Square Test for Linear-by-Linear Association (SPSS 20 software, Crosstabs). The P value of this test is applied to validate the credibility of the trends of [HCO₃⁻] (or WS_atm-CO2) (P 

Materials

To estimate the long-term weathering intensity in the global carbonate areas, we utilized digital global carbonate area boundary data and different high-resolution databases (meteorological parameters, vegetation index, runoff, and NPP from 1982 to 2018). The world karst aquifer map (WOKAM) produced by the Product Center of the Federal Institute for Geosciences and Natural Resources was used to extract the simulated [HCO₃⁻] values in the global carbonate areas. Due to uncertainties in precisely distinguishing limestone from the globally less common dolostone on the geological map, we calculated the carbonate weathering intensity by assuming that all carbonates are mainly calcite. The carbonate areas in WOKAM are divided into “continuous” and “discontinuous” amounts of carbonate minerals mixed with silicate minerals. Due to the very fast kinetics of carbonate weathering compared to silicates^55,56, we assumed that the weathering processes in these carbonate terrains are dominated by a pure calcite dissolution system. The NOAA CDR AVHRR NDVI V5 (5 km) was used to evaluate the spatio-temporal greening trend within the carbonate areas during 1982–2018. The original grid data were aggregated to monthly NDVI using the maximum value composite method in the Google Earth engine and then summarized to annual mean NDVI. To more accurately represent the NDVI changes around our experimental site (Puding), we extracted high-resolution MODIS monthly NDVI data near the Puding Comprehensive Karst Research and Experimental Station using MOD13A1 products. The runoff data in the carbonate areas of Southwest China were from the China Natural Runoff Dataset version 1.0, with a spatial resolution of 0.25° × 0.25° (CNRD v1.0, tpdc.ac.cn). The runoff data of major river catchments in Southwest China were obtained from national hydrological stations at the outlet of each catchment. NPP was the key input parameter for the three soil pCO₂ models used in this study. We obtained 8-day NPP data from the MUSES NPP database (MUltiscale Satellite remote Sensing products, The 8-day NPP data from MUSES were recalculated to monthly NPP values. The temperature used for [HCO₃⁻] calculations and the precipitation for RC estimations were from the CRUTEM 5.0 database ( with a spatial resolution of 0.5° × 0.5°. To estimate the RNPP for sustaining stable future [HCO₃⁻], we used the future near-surface air temperature of two scenarios (2015–2100, SSP126-Sustainable and SSP585-Development paths based on fossil fuels) from six CMIP6 models (ACCESS-ESM1-5, BCC-CS M2-MR, CanESM5, CNRM-CM6-1, IPSL-CM6A-LR, and MPI-ESM1−2-LR) outcomes from the DKRZ archive To optimize the calculation time, all the above-mentioned spatial data were re-projected onto a grid with a resolution of 0.25° × 0.25°.

Field experiment

To verify and quantify the individual influences from vegetation cover change and warming on carbonate weathering, we constructed a simulation test site (Puding Site) in the central area of Southwest China. This site is located in the Puding Comprehensive Karst Research and Experimental Station, Puding County, Guizhou Province, Southwest China (Figs. 3e and 7). This area experiences a humid sub-tropical monsoon climate with an annual mean air temperature of 15.1°C; and the mean annual precipitation is 1315 mm, 80% of which falls in the wet season from May to September. In January 2014, five concrete tanks with water pipe outlets, S1–S5, were built at the Puding experimental station, each simulating a carbonate catchment with a groundwater outlet under a single land-use scenario (Fig. 7a).

Each tank is 20-m long, 5-m wide, and 3-m deep, coated with epoxy resin (to avoid the influence of possible concrete erosion on tank hydrochemistry), filled with 2 m of dolomitic limestone rubble in the lower part, and topped with 0.5 m of soil. One water pipe was set at the bottom of each tank, which can be regarded as a groundwater outlet. Tank S1 was left with bare carbonate rocks, to provide a scenario without soil and plants. Tank S2 simulates land with bare soil and underlying carbonate rubble. Tank S3 was planted with corn to simulate intensive cultivation. Alfalfa was sown in the soil of Tank S4 to simulate grassland land use. The soil of tank S5 was first sown with a small proportion of Rosa roxburghii (a native shrub species in this area) in 2014 (Fig. 7a). To simulate land with vegetation restoration, we then allowed the natural growth of grasses and shrubs. During 2014–2021 the surface of S5 was fully covered by natural vegetation (Fig. 3f). A lateral drainage hole in each tank simulated a natural carbonate spring. Each tank system thus simulated a single land-use catchment with precisely known, identical dimensions, and with the same climatic and geological conditions^37,57. The water flow of each simulated catchment was controlled and the water temperature and [HCO₃⁻] concentration were measured monthly from 2015 to 2021. Thus, the bare rock-dominated catchment (S1) had no vegetation cover during 2015–2021, which simulates the scenario of carbonate weathering only controlled by warming. By contrast, the restored land-dominated catchment (S5) changed from a sparse shrub cover to full vegetation cover during the same period, which simulates the scenario of carbonate weathering controlled by greening and warming. The individual influences of greening and warming on [HCO₃⁻] could thus be further verified and quantified for this site.