Introduction

Several lines of evidence suggest that climate warming is inducing a decrease in global soil organic carbon (SOC) stocks, mainly by accelerating soil microbial activity1,2. Despite huge uncertainties associated with the magnitude of SOC loss3,4,5, it is acknowledged as a powerful climate-carbon cycle feedback6. For this reason, modelling efforts on various spatial scales have been undertaken to estimate potential losses under different future warming scenarios6,7,8. It has been shown that SOC sequestration efforts for climate change mitigation, focusing on agricultural soils in particular, can be severely hampered by climate change per se6,7,9.

Several agricultural soil monitoring networks on a national to continental scale have reported SOC losses that are partly interpreted as being caused by climate change10,11. In one prominent case, the observed loss of SOC in agricultural soils in England and Wales has been attributed to climate change10. This interpretation was challenged in a subsequent study12, which found changes in management to be a more obvious cause of SOC depletion. This is likely to be the case in many parts of the world. Between 2009 and 2015, average losses in cropland SOC were also detected by the LUCAS soil inventory across the European Union and similarly in different member states of the European Union11,13,14,15. Despite a strong increase in yields and thus in biomass production when compared with pre-industrial agriculture, intensification may have caused declines in SOC due to agricultural management. For instance, high appropriation of net primary production, breeding to maximise C allocation to the harvested product, the application of pesticides, drainage, erosion, soil compaction and very frequent soil disturbances have most likely contributed to SOC depletion to varying extents12,14,16,17,18. Furthermore, land-use change can have long-lasting legacy effects on SOC trends19. It is thus likely that SOC stocks in many agricultural soils are not in equilibrium under current management practices14. Simultaneously, the global average temperature has risen since the beginning of the last century by approximately 1.07 °C20. Globally, the spatial variability of measured climate change is huge, with a tendency towards greater warming in higher latitudes. This reflects the difficulty of decoupling the effects of agricultural management from climate change-driven background losses. While efforts have been undertaken to quantify the total global soil carbon debt of land use21, no such assessment exists for past climate change.

In 2019, the Food and Agriculture Organization of the United Nations (FAO) published a map of global SOC stocks that was prepared with the participation of numerous countries22. Based on this map, the FAO subsequently created a framework to estimate the global carbon sequestration potential of agricultural soils23, which is high on the political agenda as a nature-based solution to combat global changes. This participatory framework (GSOCseq) used the well-known and widely applied SOC model RothC24 in combination with the simple and frequently used net primary production (NPP) model MIAMI25 to estimate how much carbon can additionally be stored in agricultural soils under certain simple scenarios of carbon input alterations. This pragmatic framework was established by leading experts in SOC modelling and can be readily applied in various contexts.

Here, the GSOCseq framework was used to derive the first spatially explicit estimate of past climate change effects on SOC stocks in global agricultural soils. Starting from today’s SOC baseline as derived from the FAO GSOC map, one century of climate data was used to backwards predict SOC trends for a total of 932 k points, equally distributed over the agricultural area in all climatic zones. Modelling was restricted to the area of applicability of the model framework, thus excluding organic, forest and urban soils. Two scenarios were modelled: in the first scenario, the steady state C input of 2018 was kept constant throughout the time series, so that only the climate change effect on SOC decomposition would be simulated (NPPconst). In the second scenario (NPPvar), the C input was scaled by the modelled NPP response to climate change (temperature and moisture only) to estimate how much of the SOC changes in NPPconst could have been compensated by increased biomass production. Those may serve as two different, context-specific baseline SOC trends of the past 100 years with a high spatial resolution.

Results

On average, global air temperature at the assessed data points increased by 1.03 °C between 1919 and 2018. Depending on the scenario, this has led to an average decline in total SOC stocks by 2.5 ± 2.3 Mg ha−1 (NPPconst) and 1.6 ± 3.4 Mg ha−1 (NPPvar) (Fig. 1a,b). The comparison of both scenarios shows that climate-driven changes in NPP (excluding CO2 fertilisation) were able to compensate for 36% of predicted SOC losses on average. For both scenarios, the majority of sites had losses in the range of 0–3 Mg ha−1 (Fig. 1c, 57% in NPPconst and 52% in NPPvar). In line with the temperature trend, SOC losses were more pronounced in the latter half of the simulated timespan. A marked decline in average SOC stocks has been particularly evident since the late 1960s (Fig. 1d). Two different change rates were therefore calculated, i.e. for the first half and latter half of the timespan. This led to estimates of average climate change-driven SOC losses of 0.019 ± 0.045 (NPPconst) and 0.010 ± 0.037 (NPPvar) Mg ha−1 year−1 for 1919–1968 and 0.032 ± 0.059 and 0.019 ± 0.045 Mg ha−1 year−1 for 1969–2018 respectively.

Figure 1
figure 1

(a) Average climate change-driven SOC stock change since 1919 with standard deviations assuming constant C input (NPPconst). (b) Average climate change-driven SOC stock change since 1919 with standard deviations assuming variable C input (NPPvar). (c) Density distribution of SOC stock changes at all 932 k modelled points for both NPP scenarios. (d) Average global SOC stock of the past century at ten-year intervals as modelled from 2018 to 1919 with two C-input scenarios.

Full size image

SOC changes in response to warming also varied greatly between climatic zones (Table 1). Cold regions with high initial SOC stocks and the most pronounced increases in air temperature were predicted to have the highest SOC losses per area (Fig. 2a,b). However, also relative losses tended to be highest in colder regions (Table 2). Temperate ecoregions showed intermediate losses, and both tropical and arid regions tended to be the least affected (Tables 1, 2, Fig. S1). When NPP and thus biomass input was kept constant, on average SOC was lost from all climatic zones, ranging from − 10.8 ± 7.2 Mg C ha−1 in the temperate Cfc climate to − 0.6 ± 1.2 Mg C ha−1 in the hot arid BSh climate (Table 1). However, SOC losses were not only driven by the initial SOC stocks and mean annual temperature of the respective climate zone; relative SOC changes were also correlated to changes in water balance and temperature, i.e. higher losses were observed when the climate became wetter and warmer (Fig. 2c,d). Accordingly, the largest SOC losses in the NPPconst scenario could be expected where high initial SOC stocks, high increases in mean annual temperature (MAT) and mean annual precipitation (MAP) coincide (ΔSOC = 4.501–0.057 × SOCinitial − 4.318 × ΔMAT − 0.019 × ΔMAP, R2 = 0.74). When climate change-driven NPP alterations were considered, on average SOC was also lost from all climatic zones. In the NPPvar scenario, SOC changes in the past century ranged from − 6.4 ± 7.0 Mg C ha−1 in a cold humid Dfa climate to − 01 ± 2.3 Mg C ha−1 in a hot arid Bsh climate. However, variability across climatic zones was less well explained by the considered variables, most likely due to similar effects of temperature and moisture change on C input and C mineralisation. Initial SOC was the sole explanatory variable of the best model (R2 = 0.43).

Table 1 Average absolute soil organic carbon (SOC) stock changes for 1919–2018 (100 years) and 1968–2018 (50 years) for the scenarios NPPconst (constant net primary production) and NPPvar (variable net primary production) for all climatic zones based on the Köppen–Geiger classification with standard deviation (sd) and number of grid cells (n).
Full size table
Figure 2
figure 2

Average centennial (1919–2018) changes in soil organic carbon (SOC) stocks per climatic zone in the scenario with constant net primary production (NPPconst) as a function of (a) initial SOC stock, (b) mean annual temperature (MAT), (c) changes in water balance and (d) degree of warming. Absolute (a,b) and relative (c,d) changes are shown.

Full size image
Table 2 Average relative soil organic carbon (SOC) stock changes for 1919–2018 (100 years) and 1968–2018 (50 years) for the scenarios NPPconst (constant net primary production) and NPPvar (variable net primary production) for all climatic zones based on the Köppen–Geiger classification with standard deviation (sd) and number of grid cells (n).
Full size table

When mapped globally, spatial patterns of SOC stock change revealed increases in areas with negative temperature trends (cooling). The most pronounced SOC increases could be seen in south-west Canada, south-east USA and parts of South America, Africa, China and Australia (Fig. 3). However, in some extreme cases, such as in parts of Australia, even the negative trend in precipitation caused the model to predict an increase in SOC, with the decay of unchanged organic matter input becoming water-limited (Fig. S2). This is in line with climatic zones in which the water balance became more negative tending to lose less SOC (Fig. 2). This was not the case for the NPPvar scenario, since the soil moisture responses of plants and soil microorganisms have opposing effects on SOC stocks. Therefore, the difference between SOC losses in the NPPconst and NPPvar scenarios was negatively correlated to changes in water balance (Fig. 4). This indicates that in regions that became drier over the past century, NPP was negatively affected and thus more SOC was lost in NPPvar compared with NPPconst. However, in the majority of climate zones (20 out of 27), increased plant growth partly compensated for SOC losses (Table 1).

Figure 3
figure 3

Global distribution of climate change-driven soil organic carbon (SOC) stock changes in agricultural topsoils between 1919 and 2018 for scenarios with (a) constant net primary production (NPPconst) and (b) variable net primary production (NPPvar). White areas are those with non-agricultural land cover. 0 changes (grey) indicate changes of 0±0.3 Mg C ha−1 and all other values indicate changes ranging up to the given number. Maps were created in R, version 4.1.1 (https://cran.r-project.org/) using the package terra26.

Full size image
Figure 4
figure 4

Difference in soil organic carbon (SOC) stock changes between the two model scenarios (constant net primary production, NPPconst, and variable NPP, NPPvar) for each climatic zone as a function of the change in water balance (a positive change implies more available water). Linear regression with standard deviation is shown.

Full size image

A geothermal warming experiment was used to validate the scenarios used in this study. Figure 5 clearly depicts, that RothC, as well as the combination of RothC and MIAMI were able to predict the general trend in SOC stocks with relatively high accuracy. Relative root mean square errors (RMSE) were 26% for NPPconst and 19% for NPPvar. The NPPvar scenario was slightly closer to the observed results, which can be explained by the fact that the positive impact of warming on NPP partly compensated for enhanced mineralisation. This plausibility check increases the confidence in the global modelling results.

Figure 5
figure 5

Measured and modelled steady-state soil organic carbon (SOC) stocks along a geothermal warming gradient on an Icelandic grassland soil. Boxplots represent measured data (n = 5) with black dots indicating outliers, while blue and organge dots represent the two modelling scenarios assuming constant C input (NPPconst) and variable C input (NPPvar).

Full size image

Discussion

This is the first time that a comprehensive and spatially explicit estimate is provided on the contribution of climate change to global agricultural SOC stock dynamics over the past century. As expected, the average impact was negative and more pronounced in the latter half of the century, in accordance with the global average temperature trend (IPCC). Over the course of 100 years in which the global air temperature was warmed by on average 1.03 °C, agricultural soils lost on average 1.6–2.5 Mg C ha−1 (depending on the scenario) or were relatively depleted by 2.5–3.9%. This is in line with the order of magnitude observed in the very few long-term warming experiments that have been undertaken. In a century-long geothermal soil warming gradient in subarctic Canada, whole-soil SOC stocks have been found to decline by 3% per 1 °C27 in a deciduous forest. This was also observed for a grassland topsoil in a geothermal gradient in Iceland (2.8 Mg C ha−1 °C)28. Finally, the Harvard Forest soil, the longest-running soil warming experiment, lost about 3 Mg C ha−1 °C−1 in 26 years of warming, which is also similar to the estimated absolute values in the present study29. Long-term warming experiments in agricultural systems are scarce, so we used the mentioned Icelandic, geothermal warming experiment on a semi-natural grassland soil to validate the applied modelling approaches. Despite the facts that i) the studied soil was an Andosol which can have distinct properties regarding SOC stabilisation, (ii) warming occurred abruptly, from below and up to extreme temperatures, both modelling approaches were predicting the observations surprisingly well. The NPPconst scenario tended to slightly overestimate SOC losses, indicating that NPP and C inputs were also affected by soil warming. The validation exercise increased our confidence in the model results and proves that also relatively simple first-order kinetics models can be applied to estimate climate-change driven SOC stock changes.

With an average absolute change rate of − 0.019 (NPPvar) and − 0.032 Mg C ha−1 year−1 (NPPconst) since 1968, the results of this study suggest that the contribution of climate change to observed trends in SOC stocks is relatively small, yet not negligible. In relative terms, SOC was lost at rates of 0.03–0.05% per year, which is one order of magnitude below the SOC losses (0.6% year−1) observed in England and Wales in 1978–2003, for example10. However, the values are well in line with a model estimate of the climate change-driven annual SOC loss for the same region (− 0.08, − 0.04 and − 0.05% for cropland, grassland and forest mineral soils respectively without taking NPP changes into account)12. Both comparisons, i.e. with the long-term warming experiments (including the validation site) and the regional model exercise, indicate that the results of the present study are realistic. Thus, climate change is having an evidently persistent and extensive negative effect on global SOC stocks, but is about one order of magnitude lower than SOC losses observed in inventories of past decades10,13,30.

SOC dynamics are strongly driven by soil temperature and moisture. Both are highly variable on a global scale, on average and in their temporal trends over the last century. As expected, this has resulted in a diverse pattern of estimated changes in SOC stocks. Gains in SOC in the NPPconst scenario can be explained by either decreasing temperature or, in some cases, a more negative water balance (declining water availability), both of which can limit SOC turnover31. In some areas, the stimulation of NPP slightly overcompensated for increased SOC mineralisation, causing an accumulation of SOC with past climate change. This fits well with the fact that biomass production is expected to increase most in colder climates32. The combination of the two scenarios highlights the complexity of spatially explicit prediction of this particular climate-carbon cycle feedback33. However, in the modelling framework used here, no average SOC gains per climate zone were predicted for either of the two scenarios, pointing to the fact that increased NPP and thus C input is unlikely to compensate for SOC losses due to stimulated heterotrophic respiration34.

This study focused on agricultural mineral soils since they are the main target area for SOC sequestration measures and the area of applicability for the modelling framework23. Forest soils, organic soils and soils in other ecosystems such as wetlands and urban areas were not modelled. Although the directions of change might be similar in these systems, it was not possible to derive estimates for the whole land surface. Therefore, no attempt was made to estimate total climate change-driven losses of SOC, which would be a rather incomplete and imprecise estimate, also because the agricultural area has been far from constant in the past century35. Furthermore, the RothC model works in a monthly timestep, which flattens out climate extremes that might have more severely affected NPP and microbial activity36. However, it is believed that such short-term effects would not have changed the magnitude of change over the course of one century. Increased atmospheric CO2 is also a feature of climate change and has a fertilising effect on NPP37. Thus, the NPP response was probably slightly underestimated. However, several studies showed that the effect might be rather insignificant38,39. Moreover, backwards modelling was undertaken based on a steady-state assumption of SOC stocks to obtain the annual C inputs for the year 2018. This assumption was a rough approximation for the majority of soils, which may have led to an incorrect initial pool distribution. However, it is believed that this is not problematic as century-scale model projections have been shown to be very robust to variations in initial pool distributions40. Moreover, the steady-state was a necessary precondition in this study to exclude any legacy effect and clearly focus on climate change as the only driver of SOC dynamics. It should also be highlighted again, that the modelling approach deliberately ignored all changes in agricultural land use and management in order to isolate the climate change effect on SOC stocks from all other effects. Consequently, the modelled values are of a theoretical nature, particularly for the NPPvar scenario. The technological achievements in the last century had huge impacts on agricultural NPP41. This was not all accounted for in the present study, while the SOC stocks used to initialise the model runs (2018) were most likely more depleted than they would have been if only climate change affected them in the past century. In consequence, since SOC stock changes are a function of initial SOC stocks42, the absolute estimates in this study for the past century are slightly lower than they might have been in reality. However, (i) this is not true for reported relative changes and (ii) the advantage of using current SOC stocks is that the absolute estimates of climate change-driven alterations in the past 50 years can be expected to be close to the rates of today and the near future, since global warming has been a linear process since that time. Finally, Bradford et al. mentioned that many biogeochemical models, which are most often based on first-order kinetics, partly represent outdated knowledge of C turnover5. They argue that only the realisation of new ideas about C cycling in soils can increase confidence in modelling results regarding the effects of global warming. One specific example is the physiological response of microbes to warming, which may or may not decrease their growth efficiency2,43,44,45,46 and thus may or may not affect the size of the microbial biomass pool and also SOC stabilisation and mineralisation47. Indeed, experimental results on such specific mechanisms tend to scatter in all possible directions, as do ensembles of more complex biogeochemical models when applied to data of experimental warming48. Bradford et al. estimated, that a timespan of about 25 years would be needed to reduce current uncertainties and build more confidence on the overall effect of global warming on SOC stocks5. Currently, we see the use of a globally applied first-order kinetics model as a conservative and relatively robust way for such a global modelling effort with results being comparable to the multiple contexts in which RothC or similar models are still state of the art (such as national greenhouse gas inventories). We consider it unlikely, that the steady increase in temperature of on average 1 °C strongly altered microbial physiology and SOC stabilisation mechanisms2,44 and also tested the approach on an experimental dataset that was predicted surprisingly well.

Despite several shortcomings and uncertainties, this study provides a robust and spatially highly resolved estimate of climate change-driven SOC dynamics over the past century. This is an important product that may help interpreting observed local SOC changes and provide a certain reference scenario in various contexts. One example could be carbon farming, where currently the discussion is much about the quantification and verification of SOC sequestration on a certain land area49. In fact, not all measures implemented will lead to an absolute increase in SOC stocks, but rather avoid SOC losses. Without a seriously estimated reference scenario, it cannot be quantified how much CO2 loss due to climate change was effectively avoided. This is taken up by carbon farming start-ups, seeing a chance to sell certificates even after detecting no changes in SOC over time (https://co2-land.org/?page_id=60). The danger of a too negative reference scenario is high. A similar problem exists in forests, in which the baseline scenario (how much of the area would potentially be deforested) determines how much CO2 certificates can be sold when the forest is protected. The presented comprehensive estimation of global SOC stock changes attributable to climate change as an omnipresent global phenomenon can serve as a calibration baseline, while it is recommended to use more exact local data and modelling approaches to estimate climate change and possibly other impacts on SOC trends in carbon certification schemes.

Conclusions

The last one hundred years have seen rapid transformation. Technological developments, such as synthetic N fixation, as well as a rapidly growing human population have paved the way for the onset of large-scale industrial agriculture as well as for ever-increasing CO2 emissions responsible for the greenhouse effect50. The latter has led to the most severe environmental crisis humanity has ever faced: the climate crisis. Both agricultural intensification and climate change have simultaneously affected agricultural SOC stocks to varying extents, leading to a significant feedback to atmospheric CO2. Despite the tremendous efforts being made to quantify changes in SOC stocks in almost all regions of the world, drivers for such changes often remain elusive and, until now, were a blind spot in national greenhouse gas inventories51. The contribution of climate change is often subject to speculation and either ignored or inflated. It is ignored when emission factors rather than dynamic models are used to estimate SOC stock change due to land-use or management change52, and potentially inflated when changes in observed SOC stocks are explained10. The global estimates presented here provide a reliable foundation for a variety of applications into which the persistent, climate change-driven baseline trend of SOC stocks can be integrated.

Materials and methods

The RothC24 model was applied to predict climate change-driven SOC changes on global agricultural soils using the RothC version implementation in the R package soilassessment53. It is an established first-order kinetics model with five C pools that has been successfully applied across multiple climatic zones and also in climate change projections54,55,56. To estimate only the climate change impact on SOC stocks of today’s global agricultural soils, we had to start with a steady state situation excluding any potential legacy effects on SOC trends. Global SOC stocks are only available in sufficient quality and coverage for recent years, and also the agricultural area of today is not comparable to the agricultural area 100 years ago. Therefore, we initialised the model with todays’ SOC stocks and used temperature and precipitation data of the past 100 years to model SOC development, as affected by climate change only, backwards. In this way, all major agricultural developments in the past century, as well as land use changes were explicitly excluded from the modelling exercise and the model only predicted what a “cooling” effect (warming backwards) might have on SOC stocks. The GSOCmap (FAO-ITPS, 2019, 1 km resolution) was used to specify the initial value of the soil carbon stocks for the 932 k modelled points. In addition to carbon stocks, model driving data were used that are similar to the data compilation suggested by the FAO initiative “Global soil organic carbon sequestration potential”23:

  • clay contents in the 0–30 cm topsoil layer taken from soil texture maps of the OpenLandMap with a 1 × 1 km resolution (https://doi.org/10.5281/zenodo.1476854)

  • monthly precipitation, mean temperature and potential evapotranspiration from the CRU TS v 4.03 raster dataset with a resolution of 0.5°40

  • land use of 2018 from the European Space Agency (ESA) Climate Change Initiative (CCI) – Copernicus Climate Change Service in a resolution of 300 m

  • monthly land cover derived by NVDI from MODIS–MOD13A2 datasets (https://lpdaac.usgs.gov/products/mod13a2v006/), after the method suggested in23

Prior to model runs, the model was initialised by spin-up runs to derive carbon input at equilibrium (Cinequi) and related pool distributions in 2018. These spin-up runs were done with an analytical solution of RothC42 to minimise computational time.

Two model scenarios were run, and both explicitly ignored any changes in agricultural practices on both SOC decomposition and C inputs. Instead, only potential climate change-driven alterations in SOC were modelled (Fig. 6). In scenario 1 (NPPconst), a constant annual carbon input similar to the input at equilibrium was assumed. In scenario 2 (NPPvar), the annual carbon input from 2018 to 1919 was derived by scaling the Cinequi using the ratio of the recent NPP (NPP(t)) and NPP in the reference period 1919–2018 (NPPref):

Figure 6
figure 6

Illustration of the modelling framework with exemplary average soil organic carbon (SOC) stock trends of the two scenarios NPPconst and NPPvar as well as the actual, yet unknown trend in SOC stocks of today’s agricultural soils.

Full size image
$$Cin \left(t\right)=Ci{n}_{equi}\times \frac{NPP \left(t\right)}{NP{P}_{ref}}$$

This is not exactly the same approach applied in23, however both approaches produce identical C-input estimates. NPP needed to scale C inputs for the RothC model was estimated by the simple and well established MIAMI model57 based on annual precipitation (P) and annual mean temperature (T):

$$NP{P}_{T}=\frac{3000}{1+exp\left(1.315-0.0119T\right)}$$
$${NPP}_{P}=3000\left(1-exp\left(-0.0000664P\right)\right)$$
$$NPP=min\left(NP{P}_{T}, NP{P}_{p}\right)$$

Monthly soil water deficit, required to derive the rate modifying factor b24 in the RothC model, was quantified in forward mode starting in January 1919. There are several reasons, why those two scenarios were modelled. First of all, there might be different definitions of what is the baseline of climate change-driven SOC dynamics. One could be interested on the accelerated decomposition alone, which might then lead to considerations of how much more C input would be required to (over)compensate such a loss6. However, the question could also be, how climate change might have potentially affected the total balance of C input to the soil and SOC mineralisation. This is particularly interesting, when considering the spatial resolution of such alterations in both fluxes. There might be regions in which climate change has actually led to increased SOC stocks due to an overcompensation of enhanced SOC mineralisation by higher photosynthetic activity, while in other regions e.g. increasing drought led to a decline in potential C input. Figure 6 illustrates the modelling approach.

A unique geothermal warming experiment on Iceland was used to validate the model scenarios applied in this study. For lack of a long-term warming experiment in an agricultural system, we used this semi-natural grassland site, on which five transects along strong soil warming gradients were sampled and analysed for SOC contents58. The site is located close to the village of Hveragerdi (64° 00′ 0ʺ N, 21° 11′ 09ʺ W) on an Andosol with the topsoil (0–10 cm) having a loamy texture (8% clay, 61% silt and 31% sand), a pH in water of 5.7 and an average bulk density of 0.6 g cm358,59. At the time of sampling, soils were warmed for only six years, by on average 0.6, 1.8, 3.9, 9.9, 16.3 and 40 °C, which was due to an abrupt shift in the location of geothermal channels59. However, a later resampling on several of these warming levels plus a sampling of a neighbouring older warming experiment showed that six years of strong and abrupt soil warming was enough to reach a new steady state28. Therefore, the RothC model was initialised with the SOC stock of the unwarmed reference plots and run into a new steady state for each of the warming intensities adding the average temperature increase to the sampled weather data from 1998 to 2018, while the SOC stock of the unwarmed reference soil and mentioned climate data was used to estimate steady-state C input. This was done for both scenarios (NPPconst and NPPvar) in the exact same way as described earlier.

Owing to computational limitations, scenario runs were done for a sample of 1% of all raster grids selected by random sampling, which resulted in about 932,389 runs given the 1 km × 1 km resolution of the underlying SOC map. Modelled SOC, carbon input, NPP, weather data and soil moisture deficit were stored for 10-year time intervals. These variables were aggregated according to the Köppen–Geiger climate classification map60.

Raster maps of modelled soil carbon changes were spatially aggregated to a 0.1° resolution using the terra package in R61.

Linear regression models were fitted to explain the variability of SOC stock changes (absolute and relative) across climatic zones. Initial SOC stocks, mean annual temperature (MAT), mean annual precipitation (MAP), water balance as well as changes in temperature, precipitation and water balance were used as explanatory variables. The best model was chosen based on the Akaike Information Criterion (AIC) and model residuals were visually checked for normal distribution using quantile–quantile plots.