close
close

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 × 105 km2, mainly limestone) and diverse karst landscapes32. This area is dominated by a subtropical monsoon climate and rainfall events are concentrated from April to September49. In the last century, deforestation, and damaging cultivation practices caused widespread rock desertification and ecosystem degradation in the carbonate areas of Southwest China50. Hence, China has invested nearly 19 billion USD to implement several national ERPs in this region since the 1980s34. 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 decades33,35 (Fig. 1a, b). In this study, we focused on the carbonate areas in Southwest China and simulated the variations in [HCO3] and WSatm-CO2 during this greening period. Additionally, we simulated [HCO3] in the global carbonate areas during the same period. Details of the biogeochemical models (NPP–pCO2 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 [HCO3] as a function of temperature and soil pCO210. In this study we used the results of this carbonate equilibrium model (simulated [HCO3]) to represent the actual [HCO3] at the global scale. This model can also separately estimate the individual influences from thermodynamics (temperature) and soil pCO2 variations. Equation 1 is a simplified equation to represent how carbonate weathering captures atmospheric CO2. Actually, for a pure calcite dissolution system, the weathering of carbonate is controlled by a series of chemical equations2,3,13:

$${{CaCO}}_{3(s)}={{Ca}}_{({aq})}^{2+}+{{CO}}_{3({aq})}^{2-}$$

(2)

$${\,{HCO}}_{3({aq})}^{-}={{CO}}_{3({aq})}^{2-}+{H}_{({aq})}^{+}$$

(3)

$${\,{CO}}_{2({aq})}^{*}+{H}_{2}{O}_{(1)}={{HCO}}_{3({aq})}^{-}+{H}_{({aq})}^{+}$$

(4)

$${{CO}}_{2(g)}^{*}={{CO}}_{2({aq})}^{*}$$

(5)

$${H}_{2}{O}_{(1)}=\,{H}_{({aq})}^{+}+{{OH}}_{({aq})}^{-}$$

(6)

These equations are governed by different mass action laws, which are controlled by the temperature-dependent reaction constants K1, KC, KH, and K2, and the temperature, T (°C)13. K1 and K2 are the first and second dissociation constants for carbonic acid in water (Eq. 4 and Eq. 3), and Kc is the equilibrium constant for calcite dissolution (Eq. 2). KH is the Henry’s Law constant for CO2 gas in water (Eq. 5)14. These reactions are characterized by rapid kinetics and the drainage water can reach an equilibrium state within a couple of hours13. The pCO2 (atm) at the water–rock interface (atmosphere or soil) is an additional determinant of the magnitude of HCO3 production13.

In nature, the water interacting with CaCO3 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 coefficients3:

$$\log (\gamma )=-{{{\rm{A}}}}{z}^{2}\left(\frac{\sqrt{I}}{1+\sqrt{I}}-0.3\sqrt{I}\right)$$

(7)

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 (Ca2+, Mg2+, Na+, K+, HCO3, CO32−, SO42−, NO3, and Cl), then the electro-neutrality of its composition implies that:

$$\begin{array}{c}2[{{{{\rm{Ca}}}}}^{2+}]+2[{{{{\rm{Mg}}}}}^{2+}]+[{{{{\rm{Na}}}}}^{+}]+[{{{{\rm{K}}}}}^{+}]+[{{{{\rm{H}}}}}^{+}]\\=2[{{{{\rm{CO}}}}}_{3}^{-}]+2[{{{{\rm{SO}}}}}_{4}^{2-}]+[{{{{\rm{HCO}}}}}_{3}^{-}]+[{{{{\rm{NO}}}}}_{3}^{-}]+[{{{{\rm{Cl}}}}}^{-}]+[{{{{\rm{OH}}}}}^{-}]\end{array}$$

(8)

We can then group the species with no acid-base behavior and neglect H+, OH and CO32−, and Eq. 7 can be modified to:

$$2[{{{{\rm{Ca}}}}}^{2+}]=[{{{{\rm{HCO}}}}}_{3}^{-}]+{RA}$$

(9)

with \({RA}=2[{{{{\rm{SO}}}}}_{4}^{2-}]+[{{{{\rm{NO}}}}}_{3}^{-}]+[{{{{\rm{Cl}}}}}^{-}]-2[{{{{\rm{Mg}}}}}^{2+}]-[{{{{\rm{Na}}}}}^{+}]-[{{{{\rm{K}}}}}^{+}]\) being the reduced alkalinity3. In this study, we mainly consider the influences of greening and warming on soil pCO2 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 Ca2+, HCO3, CO32−, H+, and OH), and then the equilibrium HCO3 concentration ([HCO3]eq) can be calculated by refs. 3,13,36:

$${\left[{{{{\rm{HCO}}}}}_{3}^{-}\right]}_{{{{\rm{eq}}}}}\,={\left(\frac{{2K}_{1}\left(T\right){K}_{C}\left(T\right){K}_{H}\left(T\right)}{{K}_{2}\left(T\right){\gamma }_{{{{{\rm{Ca}}}}}^{2+}}{\gamma }_{{{{{\rm{HCO}}}}}_{3}^{-}}^{2}}{p{{{\rm{CO}}}}}_{2}\right)}^{1/3}$$

(10)

Here, all the K constants can be estimated by temperature. pCO2 (atm) is the partial pressure of CO2 at the water–rock interface, which can be estimated by soil pCO2 models. γCa2+ and γHCO3 are the activity coefficients for Ca2+ and HCO3 at equilibrium (Ca2+ and HCO3 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:

$${\left[{{{{\rm{HCO}}}}}_{3}^{-}\right]}_{{{{\rm{eq}}}}}\approx 20.75\cdot \left(1-0.0139T\right) {\root 3 \of {{p{{{\rm{CO}}}}}_{2}}}$$

(11)

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

Soil pCO2 parameterizations

Soil pCO2 is a critical parameter for [HCO3] simulation. In previous studies, soil pCO2 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 pCO2 simulation. These parameters are also controlled by vegetation cover variations. In this study, we selected a net primary productivity-based soil pCO2 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 pCO2 models are given below3:

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

$$\,{p{{{\rm{CO}}}}}_{2({soil})}=\,{p{{{\rm{CO}}}}}_{2\left({atm}\right)}+\frac{A\cdot 0.75\cdot {NPP}}{{T}^{2}}$$

(12)

where A is the set of conversion unit constants. In this study, the A values for annual and seasonal pCO2 were set to 1.03 × 106 and 1.236 × 107, pCO2atm is the atmospheric CO2 pressure (ppmv) NPP in grams of dry matter per m2 for a given year or month (g m−2 yr−1/month−1), T is the temperature (K), and pCO2(soil) is the maximum CO2 pressure reached beneath the root zone (at the soil–rock interface) (ppmv). In previous studies, we compared the results of the NPP–pCO2 model with field measurements from different global individual studies5,10. We found that the model produced very accurate results and was widely applicable in different land cover types. Hence, we used NPP–pCO2 to simulate seasonal [HCO3] variations in Southwest China and for global carbonate areas. Additionally, to evaluate the influences of vegetation production and warming on future global [HCO3], we constructed an inverse modelling framework based on NPP–pCO2 and NCEAS model (NPP = 2540/(1+exp(1.584–0.0622 T)))36. 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 CO2 sink calculation

To calculate the total atmospheric CO2 capture produced by minerals weathering, we used a weathering-related atmospheric CO2 sink equation (WSatm-CO2) [t C km−2 yr−1]2:

$$\,{{WS}}_{{{{\rm{atm}}}}-{{{\rm{CO}}}}2}={10}^{-3} \cdot ({{{\rm{n}}}}{f}_{{{{\rm{c}}}}}+{{{\rm{n}}}}{f}_{{{{\rm{s}}}}}) \cdot {[{{{{\rm{HCO}}}}}_{3}^{-}] \cdot {{{\rm{m}}}}}_{{{{\rm{r}}}}} \cdot {RF}$$

(13)

Here, WSatm-CO2 (t C km−2 yr−1) is the weathering sink for atmospheric CO2. RF (mm) is the RF, [HCO3] (mmol L−1) is the bicarbonate concentration (simulated by an NPP–pCO2 and MPD) driven by carbonates and silicates weathering, and mr (g mol−1) is the molar mass of carbon. \({f}_{{{{\rm{c}}}}}\) and \({f}_{{{{\rm{s}}}}}\) are the fractions of HCO3 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 [HCO3] is from atmospheric CO2 and the other part is from carbonate rock. In this study, we only calculate the WSatm-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 model54 by using the PyMC3 package (version 3.10) in a Python 3 environment to evaluate the sensitivity of input parameters on [HCO3] and WSatm-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 [HCO3] results (or WSatm-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 [HCO3] (or WSatm-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 [HCO3] (or WSatm-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 [HCO3] 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 silicates55,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 pCO2 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 [HCO3] 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 [HCO3], 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).

Fig. 7: The structure and the carbonate weathering simulations of the Puding experimental test site.

a Five constructed artificial carbonate catchments (each 100 m2) dominated by bare rock (S1), bare soil (S2), cropland (S3), grassland (S4), and restored land (S5). The groundwater [HCO3-] was continuously measured from 2015 to 2021. The bare rock-dominated catchment (S1) had no vegetation cover during 2015–2021, which simulates the scenario of carbonate weathering controlled only by warming. The restored land–dominated catchment (S5) experienced nearly 7 years’ vegetation restoration during 2015–2021, which simulates the scenario of carbonate weathering controlled by greening and warming. b and c are drone photos for the Puding Comprehensive Karst Research and Experimental Station and experimental test site, respectively.

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 conditions37,57. The water flow of each simulated catchment was controlled and the water temperature and [HCO3] 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 [HCO3] could thus be further verified and quantified for this site.