Quantifying the impacts of land cover change on gross primary productivity globally
GPP dataAs our primary productivity product we used the GOSIF GPP dataset21 which utilizes the linear relationship between GPP and remotely-sensed SIF34. GOSIF GPP is available globally at 0.05° spatial resolution for the period 2000–2021, with the period 2001–2015 selected here (for a short summary of all datasets used in this study see Supplementary Table 3). GOSIF GPP is based on a gridded SIF product (GOSIF)34 which uses MODIS enhanced vegetation index and meteorological data for spatial scaling and is trained with millions of SIF observations from the coarser-resolution Orbiting Carbon Observatory-235. The global coverage of GOSIF and the close relationship between SIF and GPP allow for an independent assessment of how land cover changes affect GPP in different regions around the world. For instance, SIF has been shown to capture the high GPP in the US Corn Belt derived from flux towers, while ecosystem models underestimated it36. While GPP can thus be empirically estimated from satellite SIF observations relatively reliably (even though some assumptions like the linear GPP–SIF relationship and its universality across biomes are still debated20,37,38,39), the calculation of NPP needs additional assumptions of autotrophic respiration. Therefore, we focused our study on GPP, but we included an NPP product in our uncertainty analysis. In addition to that, to account for the challenges and uncertainties in global GPP estimates we included four alternative GPP products in our sensitivity analysis (see below).Land cover mappingGridded land cover was derived from ESA-CCI22, a global land cover product designed for climate science. ESA-CCI is available at 300 m spatial resolution for the 1992–2020 period (https://cds.climate.copernicus.eu/). We first classified ESA-CCI land covers to forests, grasslands, and croplands according to IPCC classification: classes 50–100, 160, 170 forests (2,022,283 grid cells); classes 110 and 130 grasslands (509,297 grid cells); classes 10–40 croplands (950,025 grid cells). We focus on these three major land cover types to facilitate our analysis. We then converted the resulting map to 0.05° resolution by determining the prevalent (i.e., mode) land cover for each grid cell using the aggregate function from the raster package40 and only included grid cells in our training data in which the prevalent land cover was constant over the period 2001–2015. Other classes (e.g., cropland/natural vegetation mosaics) and grid cells where the land cover changed over the 2001–2015 period were not used for the RF training.Random forestsRF is a popular and efficient supervised machine learning technique which can be applied for classification and regression problems41. While complex, it is still easier to interpret compared to other machine learning methods such as Artificial Neural Networks. It has recently been applied to a wide range of ecological research questions, including the prediction of food42 and bioenergy43 crop yields, potential natural vegetation31, forest aboveground biomass44, soil respiration45, and soil carbon emissions from land-use change5 and is thus well suited for our approach. The “Forests” refer to a number of individual decision trees. For each tree, a random sample of the training data is selected and split multiple times based on a random subset of variables from which the one minimizing the weighted variance is selected by the algorithm. Model performance is computed directly on out-of-bag (OOB) data which is randomly omitted from the training data (36.8% of all grid cells). When RF is applied to new data, a weighted prediction of each individual decision tree contributes to the overall prediction. Variance in the individual trees, e.g., by selecting random subsets of the observations and random variables at each node improves the overall RF predictive skill. Model training and prediction were done using the R ranger package46. After initial testing (see Supplementary Fig. S11) we decided to set the number of individual decision trees to 800 and the number of variables to possibly split at in each node to 10. As the good evaluation measures of RF algorithms can be related to spatial autocorrelation24 we also tested a coordinate-only model and performed a leave-one-out cross validation including spatial buffers (Supplementary Discussion 2, Supplementary Fig. S3). Due to the large computational effort we reduced the number of decision trees to 100 for the buffered leave-one-out cross validation.Predictor variablesWe predicted forest, grassland, and cropland potential GPP using the following 20 predictor variables in our RF algorithm: mean annual surface temperature (Tmean), mean diurnal temperature range (Tdiurnal), temperature seasonality (Tseason; standard deviation), minimum temperature of the coldest month (Tmin), annual temperature range (Tannual), mean temperature of the warmest quarter (Twarmest), mean annual precipitation (Pmean), precipitation seasonality (Pseason; coefficient of variation), precipitation of the wettest quarter (Pwettest), precipitation of the driest quarter (Pdriest), precipitation of the warmest quarter (Pwarmest), mean annual solar radiation (SR), growing degree days (GDD), relative humidity (RH), soil clay content (Clay), elevation (EL), nitrogen deposition (Ndep), nitrogen fertilization (NF), pesticide application (Pest), and gross domestic product (GDP; a proxy for agricultural management input other than NF and Pest). Overall Tmean, Tannual, and Pmean were the most important predictor variables (see Supplementary Discussion 3 and Fig. S12). We also tested other predictors (including additional bioclimatic variables, soil pH, irrigation, or phosphate fertilization) but found only negligible improvements in RF evaluation metrics and hence decided to restrict our analysis to the 20 predictors mentioned above.Climate variables were taken from the CHELSA dataset47,48, remapped to 0.05° spatial resolution using the aggregate function from the raster package40. To only include years overlapping with our GPP data we used the CHELSA time-series data for the 2001–2013 period if available and 1979–2013 climatologies elsewise. Clay was derived from the Regridded Harmonized World Soil Database v1.249. Ndep was taken from ISIMIP2b50, bilinear remapped from 0.5° to 0.05° spatial resolution using Climate Data Operators33. Elevation was obtained from WorldClim51. NF and Pest were derived from country-specific FAO data (e.g., https://ourworldindata.org/grapher/pesticide-use-per-hectare-of-cropland), i.e., we used the same value for all grid cells in a country. GDP was obtained from ref.52.Suitable areaFor the comparison of potential forest, grassland, and cropland GPP in Fig. 1g–i we only included grid cells suitable for all three land cover types. For forests, we assumed forest cover possible if the grid cell is currently forested (e.g., all grid cells of our forest training data) or if the potential natural forest cover exceeds 36.3%. This threshold represents the 5th percentile of all currently forested grid cells. Potential natural forest cover was derived from a potential natural vegetation map, available for 20 biomes at 0.00833° spatial resolution31. To convert these biomes into potential natural forest cover we assumed 100% forest cover for the ten forest biomes and 30% forest cover for tropical savannah. Other biomes were not considered. We then aggregated the map to 0.05° spatial resolution by computing the mean of 36 grid cells using the aggregate function form the raster package40 (see Supplementary Fig. S5 for the resulting map). For grasslands and croplands, we computed the 5th percentile of Tmean and Pmean in the training data (− 9.9 °C and 165 mm for grasslands and 2.7 °C and 295 mm for croplands, respectively) and removed all grid cells below those thresholds, assuming these areas to be too cold or too dry for the respective land cover type. Finally, we calculated the land cover with the highest potential GPP for all overlapping grid cells.Sensitivity analysisTo explore the sensitivity and uncertainty of our RF approach we repeated our prediction using different input datasets, potential forest cover, and machine-learning approaches. The importance of the underlying potential forest map was estimated by replacing our potential forest map (Supplementary Fig. S5) by the LUH2 potential forest map (Supplementary Fig. S13)23. To explore the dependency on the land cover product we repeated our RF prediction using the spatially aggregated MODIS land cover map (MCD12C1; IGBP scheme), available at 0.05° spatial resolution53. We classified grid cells of classes 1, 2, 3, 4, 5, (all forests), 8 (woody savannahs) and 9 (savannahs) as forest. Classes 8 and 9 were included in forest because otherwise forest cover would be underestimated in the temperate and boreal zone. Class 10 was classified as grassland and class 12 as cropland. A comparison of ESA-CCI with MODIS reveals a substantially larger cropland area in ECA-CCI but a smaller grassland area (Supplementary Fig. S14).The sensitivity to the climate product was tested by repeating our analysis using predictor variables from the WorldClim climatologies (1970–2000)51, aggregated from 30 s to 0.05° spatial resolution using the aggregate function from the raster package40. In contrast to CHELSA, growing degree days and relative humidity were not available from WorldClim but we included water vapour pressure as additional predictor.We also tested four alternative global GPP products. The vegetation photosynthesis model (VPM) product, available for the period of interest at 0.05° spatial resolution, is based on improved light use efficiency theory and is driven by remotely sensed datasets and reanalysis climate data and land cover classification which also distinguishes C3 vs. C4 photosynthesis pathways54. The second product is derived from remote sensing considering radiation and canopy conductance limitations on GPP and is available at 0.05° resolution for the 2001–2012 period55. Land cover is not an input variable. The third product, FLUXCOM, uses machine learning to scale FLUXNET site GPP to the globe56,57. FLUXCOM is available at 0.0833° resolution and was conservative remapped to 0.05° using Climate Data Operators33 meaning that the GPP of different land cover types might be mixed in regions with heterogeneous land cover patterns. The forth product is the MODIS MOD17A3 GPP product58, available for the 2001–2013 period and aggregated to 0.05° resolution using the raster package40. It is derived from meteorological data, fraction of absorbed photosynthetic active radiation/leaf area index, and land cover. As there is also a MOD17A3 NPP product available we additionally conducted a prediction for potential NPP. The MOD17A3 NPP product is calculated as GPP minus maintenance and growth respiration estimated from allometric relationships linking daily biomass and annual growth of plant tissues to leaf area index58. This leads to additional uncertainty compared to the MOD17A3 GPP product.To test the effect of an alternative RF algorithm we repeated our prediction with the RF algorithm from the Python scikit-learn library59 using the same number of decision trees (800). Additionally, we tested another machine-learning technique, a deep neural network (DNN), using the PyTorch library60. We selected 10 linear layers with 5 times alternating 128 and 256 nodes and a sigmoid output function. All layers were connected using the rectified linear unit activation function. We used the adamW optimizer with 0.0003 learning rate and 2000 epochs of training. To prevent overfitting, we included a 10% dropout after the 7th layer. Lastly, we included a very simple estimate of the most productive land cover based on the nearest neighbour using scikit-learn’s BallTree implementation together with the Haversine formula. For each grid cell we searched for the nearest forest, grassland, and cropland grid cell and assigned the respective GPP also to this grid cell. We thus assumed that environmental conditions are more or less identical in these grid cells, which might be a reasonable assumption for many locations but less reliable in complex terrain or in large homogeneous regions like the central Amazon rainforest where the nearest cropland/grassland grid cell might be located far away.Land-use change scenariosTo estimate the effects of historical and potential future land cover changes on global GPP we applied LUH2 scenarios23 which also serve as input data for ESMs participating in CMIP6. Land-use changes over the historical period are based on the HYDE reconstruction3, while future projections were developed by different Integrated Assessment Models combining various assumptions of socio-economic behaviour (SSPs) with climate mitigation targets (RCPs). Annual fractions for the two land cover classes cropland (sum of 5 crop types) and managed grassland (sum of pasture and rangeland) were available for each scenario at 0.25° resolution (https://luh.umd.edu/). We converted to 0.05° resolution assuming the same land cover fractions for all 25 grid cells around the LUH2 grid cells. We considered the following land cover transitions: forest to managed grassland, forest to cropland, and natural grassland to cropland (and reverse transitions for future scenarios). Transitions in areas suitable for only two land cover types were also included. Conversions of natural grasslands to managed grasslands were assumed not to affect productivity. We assumed the original land cover of a grid cell to be either forest (i.e., potential forest cover > 36.3%) or natural grassland and accordingly multiplied the converted areas by the differences in potential GPP derived from our RF approach. Our broad forest definition including open tree cover (see above) and the fact that we assumed a change from 100 to 0% forest area in deforested grid cells results in a total historical deforestation area substantially larger than estimated in a recent study (2.4 Mkm2 vs. 1.6 Mkm2)61. These assumptions, however, do not impair our GPP estimate as our approach implicitly accounts for gradients in forest productivity (open forests tend to have lower GPP than closed forests). To test the sensitivity of the resulting GPP reduction we also applied the potential GPP maps from our uncertainty analysis to historical land-use changes (Supplementary Fig. S6). For future land cover changes we investigated changes over the 2015–2100 period for all available LUH2 scenarios: SSP1-1.9, SSP2-2.6, SSP4-3.4, SSP5-3.4, SSP2-4.5, SSP4-6.0, SSP3-7.0, and SSP5-8.5. Land-use activities in these scenarios range from large-scale deforestation (e.g., SSP3-7.0) to reforestation (e.g., SSP1-1.9) (Supplementary Fig. S7).Earth System ModelsWe compared the potential GPP estimated by our RF algorithm to simulations of eight ESMs participating in CMIP6 (CESM2-CLM562, CNRM-ESM2.1-Surfex 8.0c63, EC-Earth3-Veg-LPJ-GUESSv464, GFDL-ESM4-GFDL-LM4.165, IPSL-CM6A-LR-ORCHIDEEv2.066, MIROC-ES2L-MATSIRO6.0 + VISIT-e ver.1.067, MPI-ESM1-2-LR-JSBACH3.2068, UKESM1-0-LL-JULES-ES-1.069) with an explicit representation of natural vegetation and at least one agricultural land cover class (cropland or managed grassland) in their vegetation sub-model. We selected these ESMs so that all vegetation models implemented in more than one ESM were represented only once (e.g., the JSBACH vegetation model is a component of both MPI-ESM1-2-LR and AWI-ESM). For each ESM, the variable gppLut was downloaded from the CMIP6 archive (https://esgf-data.dkrz.de/search/cmip6-dkrz/) for the historical simulations. These files contain simulated GPP for natural vegetation, pasture, and cropland for which we calculated the 2001–2014 mean (2014 is the last year of the historical period). ESMs use fractional land covers for each grid cell, meaning that climatic drivers are inherently the same for all land cover types within a grid cell and simulated productivities can therefore be directly compared. As ESMs differ in their spatial resolution we bilinear remapped all output to 0.05° resolution using Climate Data Operators33 to allow for a fair comparison across models. To assess the sensitivity of our results to the interpolation method we also tested conservative remapping which results in slightly different maps (Supplementary Fig. S15) and usually larger model biases (Supplementary Table 2). In addition, ESMs differ in where they simulate forests in natural vegetation areas, and therefore we removed all grid cells from the comparison where at least one ESM simulated no tree productivity/cover/biomass in order to avoid comparing the GPP of natural grasslands to managed grasslands. We provide maps based on the original output for each ESM in Supplementary Fig. S10.FLUXNET dataWe compared our predictions of potential GPP to FLUXNET Tier 1 eddy covariance measurements (Supplementary Fig. S16)70. We included all forest, woody savannah (classified as forest), grassland and cropland sites21 which were located in suitable areas for the respective land cover. Mean GPP was calculated as the mean of the GPP estimates based on the night-time (GPP_NT_VUT_REF) and day-time (GPP_DT_VUT_REF) partitioning method. As some sites only had a few years of data, all available years were considered (i.e., site mean GPP was calculated for a different time period than 2001–2015). Comparisons were made with the potential GPP in the respective grid cell in which the site was located (i.e., not calibrated to site conditions). More
