biospheremetrics v1.0.2: an R package to calculate two complementary terrestrial biosphere integrity indicators – human colonization of the biosphere (BioCol) and risk of ecosystem destabilization (EcoRisk)

Ecosystems are under multiple stressors, and impacts can be measured with multiple variables. Humans have altered mass and energy flows of basically all ecosystems on Earth towards dangerous levels. However, integrating the data and synthesizing conclusions is becoming more and more complicated. Here we present an automated and easy-to-apply R package to assess terrestrial biosphere integrity that combines two complementary metrics.

• (i)

  The BioCol metric that quantifies the human colonization pressure exerted on the biosphere through alteration and extraction (appropriation) of net primary productivity.

• (ii)

  The EcoRisk metric that quantifies biogeochemical and vegetation structural changes as a proxy for the risk of ecosystem destabilization.

Applied to simulations with the dynamic global vegetation model LPJmL5 for 1500–2016, we find that large regions presently (period 2007–2016) show modification and extraction of >20 % of the preindustrial potential net primary production. The modification (degradation) of net primary production (NPP) as a result of land use change and extraction in terms of biomass removal (e.g., from harvest) leads to drastic alterations in key ecosystem properties, which suggests a high risk of ecosystem destabilization. As a consequence of these dynamics, EcoRisk shows particularly high values in regions with intense land use and deforestation and in regions prone to impacts of climate change, such as the Arctic and boreal zone.

The metrics presented here enable spatially explicit global-scale evaluation of historical and future states of the biosphere and are designed for use by the wider scientific community, being applicable not only to assessing biosphere integrity but also to benchmarking model performance.

The package will be maintained on GitHub and through that we encourage its future application to other models and data sets.

1 Introduction

Earth system stability relies on functioning ecosystems (McKay et al., 2022), providing, e.g., carbon sequestration, moisture recycling, and resilience to disturbances and disruptions (Friedlingstein et al., 2022; Aragão, 2012; Oliver et al., 2015). The global status and functioning of ecosystems is being evaluated using a range of approaches, including assessments of human footprint (e.g., Venter et al., 2016), indicators based on empirically collected biodiversity data (e.g., Hudson et al., 2017; Newbold et al., 2016), and ecological computer models (e.g., Sakschewski et al., 2015).

Ecosystems depend on photosynthesis as major energy source, producing the primary biomass that is the foundation of almost all food webs. During the Neolithic Revolution, humanity decoupled its biomass demand from the natural cycle, which led to large-scale modification of Earth's surface (Weisdorf, 2005). Today more than 75 % of the ice-free land area of the Earth is affected by human use (Watson et al., 2018; Arneth et al., 2019). However, the level of management and appropriation differs widely from extensive (e.g., occasionally livestock-grazed steppes or extensively used forests) to intensive (e.g., machine-aided agriculture with high mineral fertilizer use and irrigation). The aggregate effect of land use on net primary production (NPP – for a list of all abbreviations see Appendix A), i.e., altered productivity and removed biomass due to agricultural and forestry harvest, is often referred to as the human appropriation of NPP (HANPP) (Haberl et al., 2004; Vitousek et al., 1986; Rojstaczer et al., 2001; Imhoff et al., 2004). Utilizing this concept, it was shown that humanity has doubled its impact during the 20th century and thereby substantially decreased the NPP remaining in ecosystems (Krausmann et al., 2013; Kastner et al., 2022).

As one consequence of land use, only 40 % of remaining forests are still characterized by high ecosystem integrity (Grantham et al., 2020). Globally, only 18.6 % of highly intact habitats are currently protected (Mokany et al., 2020a), while modeling pressure–impact relationships suggest an increasing mean species abundance loss until 2050, even under optimistic scenarios with decreasing land use (Schipper et al., 2020; Williams et al., 2021). Over the course of this century, the dilemma of negative effects from either climate change or climate change mitigation via large-scale biomass plantation highlights the need to also integrate biodiversity in planning for negative emission technologies (Hof et al., 2018).

To prevent further degradation of the biosphere and to reverse the current loss of nature, land use scenarios should take into account the regional risk for ecosystem destabilization (Rockström et al., 2021; Obura et al., 2022). We define ecosystem destabilization as a severe change in ecosystem functioning, resulting in, e.g., a decline in carbon sequestration, species composition, or water provisioning. Shifts in biogeochemical conditions can act as a proxy for this risk based on the argument that substantial changes in either basic biogeochemical properties or vegetation composition are likely to imply far-reaching, potentially self-amplifying transformations in the underlying system characteristics, food chains, and species composition (Heyder et al., 2011). The Γ metric (Gamma) proposed by Heyder et al. (2011) represents such a metric, has been used to separate the historical drivers of ecosystem change (Ostberg et al., 2015), compares the effects of climate warming and land use under future climate scenarios (Ostberg et al., 2018), and finds temperature thresholds above which severe change is to be expected with high probability (Ostberg et al., 2013). Additionally, it has been applied component-wise to model outputs from the ISIMIP fast track ensemble, indicating an increasing area under risk of severe ecosystem change with a rise in global mean temperature (Warszawski et al., 2013).

The original definition of Γ, however, did not include nitrogen variables, and the code to compute it was hardly accessible to the scientific public. HANPP, on the other hand, has so far been based on census statistics and inventory data and has not been calculated purely from vegetation model outputs.

We therefore propose an easy-to-apply R package with two complementary biosphere metrics, BioCol and EcoRisk, building on the existing indicators HANPP and Γ (Haberl et al., 2004; Heyder et al., 2011), which are here exemplarily calculated and evaluated based on simulations with the global vegetation model LPJmL5 (von Bloh et al., 2018).

BioCol quantifies the human colonization pressure on the biosphere through extraction of biomass and prevention of natural NPP by photosynthesis. The metric basically follows the HANPP approach by Haberl et al. (2007), who spatially explicitly sum extracted and inhibited biomass amounts based on biomass inventory data and compare them to potential NPP in a counterfactual world without human land use. We replace biomass inventory data with the corresponding LPJmL5 model outputs and thereby add the possibility of also computing BioCol from deep historical and future simulations.

EcoRisk illustrates state shifts in ecosystems as a result of land, water, and fertilizer use, as well as climate change based on the Γ metric (Heyder et al., 2011; Ostberg et al., 2015). It quantifies, on a scale from 0 (no change) to 1 (very strong change), the dissimilarity of an ecosystem state from a reference condition and comprises four subcomponents (vegetation structure, local change, global importance, ecosystem balance) that are aggregated as a multidimensional proxy for the risk of biosphere destabilization. We follow the original publications but (in addition to water and carbon) now include nitrogen fluxes and pools. The addition of nitrogen variables fills a major gap because nitrogen limitation or surplus is a key determinant for plant growth and the overall ecosystem status.

2 Material and methods

In this section, we detail the calculation of both BioCol and EcoRisk, our biome classification for spatial aggregation of results, and the relevant specifics of the vegetation model LPJmL5.

2.1 BioCol

We define BioCol (see Fig. 1) as the flow of biomass (in terms of NPP) that is intentionally extracted in the form of crop, residue, and other biomass harvests (NPP_harv), plus the inhibited natural biomass production as a result of land use changes, management, and differences in fires (NPP_luc). NPP_luc is calculated as the difference between potential natural NPP (i.e., the NPP that would prevail without land use but with current climate, Haberl et al., 2014, – NPP_pot) and the actual biomass production (NPP_act), which are both calculated under the same climate:

$$\begin{array}{}\text{(1)}& {\text{NPP}}_{\text{luc}}={\text{NPP}}_{\text{pot}}-{\text{NPP}}_{\text{act}}.\end{array}$$

Using time series of NPP_pot and NPP_act allows the removal of the climate change (e.g., CO₂ fertilization) effect in the component NPP_luc. NPP_luc can become negative if the actual NPP is higher than potential NPP (e.g., through management or land use legacy effects, especially on managed grasslands). Absolute HANPP and relative BioCol are computed as

$$\begin{array}{}\text{(2)}& \text{HANPP}={\text{NPP}}_{\text{luc}}+{\text{NPP}}_{\text{harv}},\end{array}$$

where NPP_harv denotes the NPP withdrawn from ecosystems via harvest.

$$\begin{array}{}\text{(3)}& \text{BioCol}=\text{HANPP}/{\text{NPP}}_{\text{ref}}\end{array}$$

Here NPP_ref represents a reference NPP. In earlier applications, the respective NPP_pot of each year has been used as reference (Kastner et al., 2022; Krausmann et al., 2013). In contrast, here we use NPP_ref, which refers to the NPP_pot of the preindustrial period (mean of 1550–1579), which is the same time frame as for EcoRisk. NPP_harv sums the corresponding LPJmL5 model outputs for harvested and extracted carbon from crop areas (including residues), second-generation biomass plantations (not existing in historical period land use input), grassland, timber harvests from land use conversion combined with external timber extraction from managed forests, and human-induced fire carbon emissions. Residue harvest was assumed to contain 70 % of the remaining above-ground biomass after harvest. Carbon extraction of grasslands includes carbon contained in dairy products and methane emissions plus respiration and is based on prescribed livestock densities (Heinke et al., 2023), which are calibrated to match grazing amounts for the year 2000 originally provided by Herrero et al. (2013) and modified by Heinke et al. (2020) (about 1.1 Gt C yr⁻¹ globally). Currently, LPJmL5 is not able to model managed forests and does not separate human-induced and lightning-induced fire carbon emissions. Therefore, timber harvest from managed forests is included as an external dataset from LUH2-v2h (Hurtt et al., 2020). Human-induced fire carbon emissions can be included through external data on the human fire ignition fraction of total fire carbon emissions based on fire models such as LPJmL-Spitfire (Thonicke et al., 2010). For this study, the fire carbon emissions are not included since the updated LPJmL-Spitfire version is still being tested.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\text{human}& \text{ fire ignition fraction}=\\ & {\displaystyle \frac{\text{human ignition}}{\text{human ignition}+\text{lightning ignition}}}\end{array}\end{array}$$

Figure 1 Calculation scheme for BioCol. The basis for our analysis is the preindustrial potential NPP (NPP_ref) assessed from 1550–1579. The effects of CO₂ fertilization of plants resulting from historical anthropogenic CO₂ emissions, changes in atmospheric N deposition, and climate change lead to a net increase in NPP today (labeled as the “CC effect” in the biosphere) both for hypothetical “potential vegetation” without human land use and the “actual vegetation” including land use. HANPP is calculated as the sum of direct human biomass extraction (NPP_harv) and inhibited natural productivity through replacing natural vegetation with land use (${\text{NPP}}_{\text{luc}}={\text{NPP}}_{\text{pot}}-{\text{NPP}}_{\text{act}}$). BioCol is subsequently computed as the fraction of HANPP compared to NPP_ref.

NPP_pot, representative of the NPP under potential natural vegetation (under transient climate), is obtained from a model simulation without human land use but with otherwise identical settings to the run providing NPP_act. BioCol or HANPP and sub-components are available as spatially explicit values per grid cell for every time step but also as global sums over time, or they are aggregated per biome or world region. When aggregating grid cell values of BioCol, negative values can be treated as such, reducing the overall pressure, or the absolute values of all grid cells can be summed up (which is mainly used in this analysis). For further details on the LPJmL5 simulation setup, see Sect. 2.5.

2.2 EcoRisk

EcoRisk is computed as the average of four subcomponents: vegetation structure (vs), local change (lc), global importance (gi), and ecosystem balance (eb), each on the same scale of 0 to 1 and internally scaled with the respective change-to-variability ratio S(x,σ_x), following Ostberg et al. (2018).

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\text{EcoRisk}& ={\displaystyle \frac{\text{vs}+\text{lc}+\text{gi}+\text{eb}}{\mathrm{4}}}\\ & ={\displaystyle \frac{\left\{\begin{array}{l}V\cdot S(V,{\mathit{\sigma }}_V)+l\cdot S(l,{\mathit{\sigma }}_l)\\ +g\cdot S(g,{\mathit{\sigma }}_g)+e\cdot S(e,{\mathit{\sigma }}_e)\end{array}\right\}}{\mathrm{4}}}\end{array}\end{array}$$

The unscaled “vegetation structure” component V is evaluated based on the dissimilarity of the vegetation composition of the whole grid cell for both natural and managed areas based on Sykes et al. (1999); Heyder et al. (2011); Ostberg et al. (2018). For this, LPJmL5 provides outputs of the ground area covered by natural and cultivated plant types. Changes in vegetation structure are computed between ecosystem state i and j with respect to the total area (G) of each ground cover type (GCT) k (tree, grass, or barren), further detailed by the differences between the area (A) specific to each plant functional type (PFT, p) regarding attribute l (termed evergreenness, needleleavedness, tropicalness, borealness, and naturalness). For PFT-specific attributes (a), see Table 1. Barren is defined without subcategories. The attribute-specific weighting factor ω_kl is per a default value set to 0.2 (= “equal” weighting). Alternatively, attribute-specific weights as in Ostberg et al. (2018) can be applied.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& V(i,j)=\mathrm{1}-\sum _k\mathit{\{}\underset{\text{total area of GCT}k}{\underbrace{min\left(G_{ik},G_{jk}\right)}}\cdot \\ & [\mathrm{1}-\sum _l(\underset{\text{specific area difference regarding attribute}l}{\underbrace{{\mathit{\omega }}_{kl}\left|\sum _p\right(A_{iklp}\cdot a_{klp})-\sum _p(A_{jklp}\cdot a_{klp}\left)\right|}}\left)\right]\mathit{\}}\end{array}\end{array}$$

Table 1 PFT-specific attributes a_klp for GCTs tree and grass.

^∗ Dynamic share due to climate specific grass mix.

The remaining components l, g, and e are computed from two ecosystem state vectors s₁ (reference state) and s₂ (changed state), composed of 30-year averages of biogeochemical variables v_i,1 and v_i,2, with $i=[\mathrm{1},\mathrm{\dots },n]$.

$$\begin{array}{}\text{(7)}& {\mathit{s}}_{\mathrm{1}}=\left(\begin{array}{c}{v}_{\mathrm{1},\mathrm{1}}\\ \mathrm{⋮}\\ v_{n,\mathrm{1}}\end{array}\right),{\mathit{s}}_{\mathrm{2}}=\left(\begin{array}{c}{v}_{\mathrm{1},\mathrm{2}}\\ \mathrm{⋮}\\ v_{n,\mathrm{2}}\end{array}\right)\end{array}$$

They represent the major process cycles on the grid cell level determining plant growth: carbon, nitrogen, and water. The system described by the metric is the terrestrial land system, which has pools, influxes, and outfluxes across the system boundary. For carbon and nitrogen, stocks are aggregated into a vegetation and soil pool, while for water there is only a soil pool. In total this gives 11 dimensions, plus a single dimension reserved for other relevant processes (for LPJmL filled with the fire frequency) (Table 2). The user can define which model output variables constitute each process dimension. Hereby, it is important to know the model-specific details to keep consistency. If possible, output variables should use per area units because of the per area weighting for the “global importance” component.

Table 2 Processes and associated variables describing the land system and their aggregation from LPJmL5 outputs.

“Local change” describes changes compared to the local reference state (e.g., the preindustrial time) as the magnitude change in the difference vector of the biogeochemical properties (how strongly the local conditions have changed). For this, the state variables are normalized with the local values of the reference state:

$$\begin{array}{}\text{(8)}& {\mathit{s}}_{l_{\mathrm{1}}}=\left(\begin{array}{c}\mathrm{1}\\ \mathrm{⋮}\\ \mathrm{1}\end{array}\right),{\mathit{s}}_{l_{\mathrm{2}}}=\left(\begin{array}{c}{v}_{\mathrm{1},l}\\ \mathrm{⋮}\\ v_{n,l}\end{array}\right),\end{array}$$

with

$$\begin{array}{}\text{(9)}& v_{i,l}={\displaystyle \frac{v_{i,\mathrm{2}}}{v_{i,\mathrm{1}}}},\text{ for }{v}_{i,\mathrm{1}}\ne \mathrm{0}.\end{array}$$

For variables that are 0 in both vectors, both are set to 1, resulting in no change. If only the reference case is 0, the unscaled values are used for both vectors.

“Global importance”, in contrast, puts these local changes in the perspective of the global mean reference condition, taking into account that even moderate changes on the local scale may feed back to larger scales if they are large enough in absolute terms. Therefore, the state vectors are normalized with the global spatially averaged reference mean value $\stackrel{\mathrm{‾}}{v_{i,\text{refg}}}$:

$$\begin{array}{}\text{(10)}& {\mathit{s}}_{g_{\mathrm{1}}}=\left(\begin{array}{c}{v}_{\mathrm{1},g,\mathrm{1}}\\ \mathrm{⋮}\\ v_{n,g,\mathrm{1}}\end{array}\right),{\mathit{s}}_{g_{\mathrm{2}}}=\left(\begin{array}{c}{v}_{\mathrm{1},g,\mathrm{2}}\\ \mathrm{⋮},\\ v_{n,g,\mathrm{2}}\end{array}\right)\end{array}$$

with

$$\begin{array}{}\text{(11)}& v_{i,g,t}={\displaystyle \frac{v_{i,t}}{\stackrel{\mathrm{‾}}{v_{i,\text{refg}}}}},\text{ for }\stackrel{\mathrm{‾}}{v_{i,\text{refg}}}={\displaystyle \frac{\mathrm{1}}{\sum a_p}}\sum _{p=\mathrm{1}}^z{a}_p{v}_{i,p}\ne \mathrm{0},\end{array}$$

for cells $p=\mathrm{1},\mathrm{\dots },z$ with cell area a_p. If the global mean reference state is 0, the mean scenario state is used for scaling instead. If both are 0, both vectors are set to 1, as was done for local change. Afterwards, state vectors for both global importance and local change are multiplied by $\frac{\mathrm{1}}{\sqrt{n}}$ to scale them down according to the number of variables (EcoRisk can also be computed for simulations without nitrogen, and thus it is missing the corresponding variables). The difference between the two states is now characterized by the length of the difference vector between them, which for local change and global importance are defined as follows:

$$\begin{array}{}\text{(12)}& d_l=|{\mathit{s}}_{l_{\mathrm{2}}}-{\mathit{s}}_{l_{\mathrm{1}}}|,d_g=|{\mathit{s}}_{g_{\mathrm{2}}}-{\mathit{s}}_{g_{\mathrm{1}}}|.\end{array}$$

“Ecosystem balance” quantifies shifts in the relative magnitude of biogeochemical properties with respect to each other as an indicator for qualitative changes in the balance of dynamic processes, which may signal a breakdown of ecological functioning (Ostberg et al., 2018). It is calculated from the angle between the two state vectors with local normalization (as for local change):

$$\begin{array}{}\text{(13)}& b^{\prime }=\mathrm{1}-\cos \left(\mathit{\alpha }\right)=\mathrm{1}-{\displaystyle \frac{{\mathit{s}}_{l_{\mathrm{1}}}\cdot {\mathit{s}}_{l_{\mathrm{2}}}}{\left|{\mathit{s}}_{l_{\mathrm{1}}}\right|\left|{\mathit{s}}_{l_{\mathrm{2}}}\right|}}.\end{array}$$

b^′ is scaled to a range between 0 and 1, assigning a value of 1 if the angle between state vectors is larger than 60°:

$$\begin{array}{}\text{(14)}& e=\left\{\begin{array}{ll}{b}^{\prime }\cdot \mathrm{2}& \text{if }\mathit{\alpha }\le \mathrm{60}\mathit{°},\\ \mathrm{1},& \text{otherwise}\end{array}\right..\end{array}$$

Values for metric components l and g are derived by scaling d_l and d_g to a range between 0 and 1 using the sigmoid transformation function T:

$$\begin{array}{}\text{(15)}& l=T\left(d_l\right),g=T\left(d_g\right),T\left(x\right)=A+{\displaystyle \frac{\mathrm{1}-A}{\mathrm{1}+e^{-\mathrm{6}(x-\mathrm{0.5})}}},\end{array}$$

with $A=-\frac{\mathrm{1}}{e^{\mathrm{3}}}$.

The year-to-year variability is accounted for by the “change-to-variability ratio”. This variability factor is based on the standard deviation (σ_x) of the component in the 30-year reference period, which is in turn based on the assumption that ecosystems are adapted to the variability they are regularly exposed to but may be vulnerable if it is exceeded. The change-to-variability ratio S(x,σ_x) for components $x\in (V,l,g,e)$ is calculated as

$$\begin{array}{}\text{(16)}& S(x,{\mathit{\sigma }}_x)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-\mathrm{4}(\frac{x}{{\mathit{\sigma }}_x}-\mathrm{2})}}},\end{array}$$

with σ_x the interannual standard deviation of x under reference conditions.

For the specific status of variables describing the same process (e.g., carbon pools, or water fluxes), the change metric is also evaluated only for those variables, by default as $(\text{lc}+\text{gi}+\text{eb})/\mathrm{3}$, but for compatibility with earlier applications it is possible to only use lc.

2.3 Comparison to other indicators

To contextualize EcoRisk and BioCol with other biosphere integrity indicators, we transformed a set of seven widely used biosphere integrity indicators (Table 3) to the interval [0,1] (with 0 meaning high integrity, no pressure, and low risk) and show them alongside EcoRisk and BioCol, together with scatterplots of EcoRisk and BioCol against the average value across all indicators. We additionally extracted thresholds for each indicator from the literature that demarcate the transition between the low- and high-risk zones (sources listed in Table 3). With these we calculated the number of indicators per grid cell that show up as transgressed.

Oldeman et al. (1990)Oldeman et al. (1990)Venter et al. (2018)Venter et al. (2016)Newbold et al. (2016)Soergel et al. (2021)Schipper et al. (2020)Grantham et al. (2020)Grantham et al. (2020)Mokany et al. (2020a)Mokany et al. (2020b)Cherlet et al. (2018)Cherlet et al. (2018)

Table 3 Sources for other biosphere integrity indicators and associated thresholds.

2.4 Biome classification

To aggregate BioCol and EcoRisk from grid cells to the biome level (to also enable analysis for these larger-scale ecological units), we updated the biome classification for LPJmL5, originally published by Ostberg et al. (2013). It is based on the vegetation structure and its fractional coverage per cell, the tree-specific leaf area index, the temperature, and the carbon stored in vegetation. Cells are classified primarily based on the total tree cover. Thresholds of 60 %, 30 %, and 10 % are chosen to demarcate the boundaries between forests, woody savanna/woodland, savanna, and grassland based on the IGBP land cover classification system (Loveland and Belward, 1997). The dominant tree or grass species then determines the type of biome. Compared to Ostberg et al. (2013), the thresholds can be adjusted and an additional differentiation between Tropical Rainforest and Warm Woody Savanna/Woodland can be either based on the tree leaf area index or the vegetation carbon (in this paper, we use a tree leaf area index of 6 as a threshold). Additionally, Montane Grassland can be differentiated from Arctic Tundra by elevation or latitude (here, we use an elevation threshold of 1000 m).

2.5 LPJmL5

We employ the dynamic global vegetation model LPJmL 5.7 (Schaphoff et al., 2018; von Bloh et al., 2018), which can be run standalone (forced by climate inputs) or coupled into the Potsdam Earth Model POEM (Drüke et al., 2021).

LPJmL5 simulates key ecological and physiological processes such as photosynthesis, respiration, carbon allocation, and turnover for natural and managed vegetation based on historical data or future projections of land use, soil, nitrogen inputs, and climate conditions. The composition of natural vegetation, represented by 11 PFTs, dynamically develops based on climatic parameters, growth efficiency, competition, and fire disturbance (Sitch et al., 2003), while agricultural crop composition is prescribed considering 12 crop functional types (CFTs) (Bondeau et al., 2007), grassland/pastures (Heinke et al., 2023), and three second-generation bioenergy crop functional types (BFTs) (Beringer et al., 2011). The remaining group of “other” crops is simulated as temperate wheat or tropical maize, depending on the latitude. CFT-specific sowing dates are dynamically calculated based on optimal season length and fixed after the year 2000 (Waha et al., 2012). Local runoff is routed through a global network of waterbodies and river channels as accumulated discharge (Gerten et al., 2004). Irrigated crop area is prescribed, but irrigation requirements and water withdrawals are dynamically simulated based on the soil water deficit and available renewable water in lakes, rivers, reservoirs, and neighboring cells with available water (Jägermeyr et al., 2015). Different agricultural management practices and their impacts on soil processes and yields are simulated, including tillage, manure, fertilizer application, and optional growth of off-season cover crops (Lutz et al., 2019; Porwollik et al., 2022). Fluxes and stocks of carbon, nitrogen, and water are by default resolved on 0.5°×0.5° spatial and daily timescale. For this study, all outputs are aggregated per year.

Output from LPJmL5 to compute BioCol and EcoRisk is read and processed using the R package lpjmlkit (Breier et al., 2023).

We ran the model with a fused climate input from the ISIMIP project (GSWP3-W5E5). It combines daily GSWP3 data from 1901 to 1978 with a bias-adjusted version of ERA5 (W5E5, adjusted to better match CRU and GPCC) from 1979 to 2016 (Kim, 2017; Lange, 2019). Manure, fertilizer, and crop-specific cultivation areas from 1500 to 2018 are taken from a new hybrid dataset (Ostberg et al., 2023). Fallow land was included in the CFT class “other”, while the distribution of irrigated area into irrigation systems is based on Jägermeyr et al. (2015). Livestock densities on grasslands were determined based on estimates of grazing levels at both regional and production system levels reported by Herrero et al. (2013) using a set of simulations presented by Heinke et al. (2023). Tillage input is based on Lutz et al. (2019). For simulation years before 1901, the first 30 years of the climate input are randomly recycled.

3 Results

The results of BioCol are (except for forest harvest) strictly model based and are in general well in line with the existing data-based assessments (Kastner et al., 2022; Krausmann et al., 2013). The global sum of HANPP, based on LPJmL5 simulated outputs, increases from below 2 [2.5] Gt C yr⁻¹ before 1700 to more than 13 [15] Gt C yr⁻¹ today, depending on whether the absolute values are summed up for cells with negative values or not (Fig. 2a – left y axis). Over the course of the 20th century, relative values of BioCol thus approximately doubled from about 0.1 to over 0.2 (10 % to over 20 %) compared to the potential preindustrial NPP_ref of the 16th century (Fig. 2a – right y axis). Since 2000 the values increased further. When taking the absolute for negative cells, BioCol is approaching 0.25 (25 %) in 2016 – meaning that almost a quarter of the total terrestrial preindustrial biomass production on Earth is rerouted to human use or inhibited compared to a world without humans.

Figure 2 (a) Global BioCol and components over time. BioCol values use the orange axis on the right and are calculated as a sum of absolute values (“BioCol abs sum” – cells with negative value increase the global sum) or simple sum (“BioCol sum” – cells with negative value reduce the global sum). (b) Map of the relative values for the year 2000 (average 1995–2005). Relative values for BioCol are expressed in comparison to the average NPP from 1550–1579 from a run without human land use.

The global relative BioCol pattern for the year 2000 (1995–2005 average) shows high values for areas of intense agricultural use (Fig. 2b). The abundance of cells with negative BioCol values (higher productivity in the simulation with land use than in the similar run with only potential natural vegetation) might be explained by the explicit simulation of agricultural management, especially in regions with low natural NPP (e.g., irrigation in the Middle East) and potential legacy effects from earlier land use on now natural areas included in our land use dataset (e.g., in the Amazon, NPP upon regrowth of natural vegetation on abandoned agricultural areas is simulated to be higher than potential NPP in the equilibrium state; see Fig. B1).

Isolated cells with very high absolute BioCol values (e.g., in the boreal zone) are locations of reservoirs established before the year 2005, with associated decline in NPP and thus high values of NPP_luc.

EcoRisk calculated between 1550–1579 and 1985–2016 shows large areas with values above 0.5 (Fig. 3), mostly in regions with high land use intensity today. Land use change is reflected in the vegetation structural change component (vs), which quite closely resembles the current land use extent (see inset in Fig. 3). However, the components building on changes in biochemical variables (lc, gi, eb) indicate a much larger extent for regions with values >0.9. Figures 4b and B2 show that in many regions changes in nitrogen fluxes are responsible for these strong changes.

Figure 3 (a) Change in biochemical compositions computed by EcoRisk between 1550–1579 and 1985–2016. (b) Current land use extent for reference. (c–f) EcoRisk components are as follows: vegetation structure change, local change, global importance, ecosystem balance.

Figure 4 (a) Present (1987–2016) biomes classified from vegetation structure, plant-specific leaf area index, temperature and elevation in LPJmL5. (b) Change in biochemical compositions computed by EcoRisk between 1550–1579 and 1985–2016 as the median (Q10 and Q90 for whiskers) across the 16 most relevant biomes (“Temperate Broadleaved Evergreen Forest” is effectively non-existent as only two cells are classified as such, while “Rocks and Ice” and “Water” are skipped for lack of vegetation). See Table B1 for biome names and abbreviations.

Climatic changes are also reflected in EcoRisk. They are most visible in the components ecosystem balance (eb) and vs as bands of higher values along the frontier of boreal vegetation onset to the North in Canada. In Eurasia (and other land-use-free regions of the Arctic) this vegetation shift appears to be more patchy. The local change component (lc) indicates high values among especially vulnerable ecosystems with low absolute values (Arctic, deserts) that do not show up in the global importance (gi). Values of gi are high in regions with strong absolute changes, e.g., with intense land use or loss of tropical rainforests. Ostberg et al. (2018) provides a more detailed discussion of separate and combined effects of land use and climate change (albeit for a somewhat differently defined metric using an earlier model version and different input database).

For the original Γ-metric, a threshold of 0.3 had been established, which indicates the transition from moderate to high risk of ecosystem destabilization. To determine whether this threshold is still valid for EcoRisk (given the adaptations to the computation, particularly the inclusion of nitrogen flows and pools that show more variability), we performed two sets of synthetic simulations. We assume that a meaningful threshold should meet the following two criteria: (i) it should be higher than internal variability within biomes but (ii) lower than the variability between distinct biomes, such that a simulated EcoRisk above this threshold is indicative of changes equivalent to a shift in biome. For (i) we checked the homogeneity within forest biomes by computing EcoRisk with values from 1550–1579 between each cell of a biome and the average cell of this biome. In all forest biomes, internal biome variability for at least half the cells is higher than 0.3 compared to the average biome cell (Fig. B3). Thus, values <0.5 could better describe most of the internal variation.

For (ii), we compared the average biome cells against each other by computing EcoRisk between the current states of each of them. Most combinations show EcoRisk values >0.3 when compared against each other (Fig. B4), with some exceptions: Boreal Needleleaved Evergreen Forest is classified as being relatively similar to Temperate Broadleaved Deciduous Forest (0.26), a fact that we cannot explain and that is also complicated by the reverse direction showing a high value of 0.62. Tropical Rainforest is classified as being similar to Tropical Deciduous Forest. This can be partially explained through their similar locations and the overlap in species. Arctic Tundra and Montane Grassland are classified as being similar. They are in fact the same biome but at different elevations. Comparably, Temperate Savanna, Temperate Woodland, and Warm Grassland are classified as similar biomes that are only differentiated by total tree cover fraction and $\mathrm{C}\mathrm{3}/\mathrm{C}\mathrm{4}$ grass shares. There are further biomes that show intermediate values of $\mathrm{0.3}<\text{EcoRisk}<\mathrm{0.5}$ when compared against each other, all of which are partially explainable through compositional similarity or similar average conditions.

However, Fig. B4 also highlights that EcoRisk is not symmetric, mainly because of the normalization to different reference conditions for local change and ecosystem balance. Strong directional discrepancies can, for example, be observed for the difference between the three biomes Temperate Needleleaved Evergreen Forest, Temperate Broadleaved Deciduous Forest, and Boreal Needleleaved Evergreen Forest.

Aggregating EcoRisk grid cell values to the biome level according to the biome classification yields three classes (Fig. 4): (i) those with a median EcoRisk<0.3, i.e., Tropical Rainforest, Tropical Deciduous Forest, Boreal Broadleaved Deciduous Forest, Boreal Needleleaved Deciduous Forest (only few cells), and Arctic Tundra; (ii) those with $\mathrm{0.3}<\text{median EcoRisk}<\mathrm{0.65}$, i.e., Boreal Needleleaved Evergreen Forest, Warm Woodland, Warm Savanna, Warm Grassland, Temperate Savanna, Temperate Grassland, Montane Grassland, and Desert; and (iii) those with a median EcoRisk>0.65, i.e., Temperate Broadleaved Deciduous Forest, Temperate Needleleaved Evergreen Forest, and Temperate Woodland. The subcomponent local change is generally the one with highest values (except for TrRF, TrDF, and BoND), and nitrogen fluxes show stronger changes than those for carbon and water in all cases.

Comparison of EcoRisk and BioCol with other biosphere integrity indicators (Table 3) shows similar trends despite high scattering (Fig. 5k and m). The maps highlight that the pattern of BioCol is more similar to that of other transgressed indicators (Fig. 5a–j); however, this might stem from most of them being directly affected by land use. The additional benefit of EcoRisk is that it also captures effects due to climate and deposition changes.

Figure 5 Contextualization of several indicators of biosphere integrity, transformed to the interval [0,1], with 0 meaning high integrity, no pressure, and low risk: (a) GLASOD representing human-induced soil degradation (Oldeman et al., 1990), (b) HF representing human footprint (Venter et al., 2016), (c) BII representing the biodiversity intactness index (Newbold et al., 2016), (d) intactness representing GLOBIOM 2015 MSA (Schipper et al., 2020), (e) FLII representing the Forest Landscape Intactness Index (Grantham et al., 2020), (f) CI representing contextual intactness (Mokany et al., 2020b), (g) CoE representing the Convergence of Evidence from World Atlas of Desertification (Cherlet et al., 2018). (h) The number of the previous seven indicators that show up as transgressed per grid cell (see Table 3 for thresholds indicating transition between low- and high-risk zones), (i) EcoRisk and (j) BioCol (l) average of metrics shown in (a)–(g), (k) scatterplot of EcoRisk versus average, and (m) scatterplot of BioCol versus average.

4 Discussion

We present a model-based indicator set that allows for the assessment of the state of the biosphere. We show that large regions presently (period 2007–2016) show modification and extraction of >20 % of the preindustrial potential net primary production according to the indicator BioCol, along with climatic changes leading to drastic alterations in key ecosystem properties and suggesting a high risk of ecosystem destabilization according to the indicator EcoRisk.

Generally, the indicators presented in this package can serve both as an analytic tool to assess biosphere integrity from model simulations or as a means of benchmarking model performance (e.g., after new development). Therefore, depending on the context, the performance of the vegetation model is important. This paper primarily describes the methodology behind the biospheremetrics package, with the application to model results only being secondary. LPJmL5 is currently in a recalibration and validation phase, with major changes to code and key processes, following the implementation of tillage and the nitrogen cycle. Further model development may particularly focus on an improved distribution of PFTs and thus biomes (see Harper et al., 2023 for comparison) and the explicit simulation of multi-cropping, which has been neglected here. A better representation of human-induced fire emissions and including process-based forest management will also be important, since the effects on NPP are currently not considered.

The next step for us is thus to extend the biospheremetrics package to be compatible with other vegetation models (e.g., utilizing outputs from the ISIMIP3 ensemble) and thus also allow for intermodel comparison.

As presented in the previous section, the addition of nitrogen fluxes leads to a strong increase in values for EcoRisk, when compared to earlier results of the Γ-Metric (Figs. 3 and B2 and Ostberg et al., 2015, 2018). Relative changes in nitrogen fluxes (Fig. B5) as a result of nitrogen fertilizer and manure application are much stronger than those for carbon and water fluxes and pools (Figs. B6 and 4). This is intuitively plausible; however, the question remains, whether the associated high values for EcoRisk mainly resulting from the relative changes in nitrogen fluxes really correspond to a strongly increased risk for ecosphere destabilization. Generally, a theory of how changes in different components (vs, lc, gi, eb) or classes of state variables (e.g., water pools or nitrogen fluxes) can be translated to risk of ecosphere destabilization is lacking (and could be different among components and classes). Theoretically, a weighted downscaling of any component would be possible; however, the literature base for such changes is currently lacking. We thus refrained from any weighting until further research is done on this end.

The results for EcoRisk show a strong increase in the overall values. A new threshold between moderate and high risk (replacing 0.3 in Ostberg et al., 2018, 2015; Heyder et al., 2011) would be as high as 0.48, when picking the “optimal” EcoRisk threshold, balancing between forest biome internal homogeneity on the one hand and inter-biome dissimilarity on the other hand (Fig. B7).

Our way to calculate BioCol differs substantially from previous approaches for HANPP. First, our calculation is exclusively based on model output, with the aim to also simulate future scenarios and historical periods outside of available biomass inventory data. However, we also use a different baseline when comparing it to the static preindustrial NPP and taking the absolute of cells with negative HANPP values. Both processes are done deliberately in order to (i) not let ourselves get distracted from the fact that global NPP is mainly rising because the biosphere is in a resilience response phase to increasing atmospheric CO₂ levels, meaning that this additional NPP is not ours to use, and (ii) acknowledge that management-driven NPP increases beyond the potential natural values (negative NPP_luc) usually go together with modification of the water, carbon, or nitrogen cycles stressing local ecosystems.

The underlying principle of BioCol, NPP appropriation, is a function of human demand, i.e., land use change and biomass harvesting for food, fiber, fodder, bioenergy, and bioeconomy. It is therefore closely associated with issues not just of Earth system stability but also of justice, access, and sustainable management of resources (see Gupta et al., 2023). The challenge is to maintain the productivity of the biosphere, i.e., the sustainably available NPP, while ensuring the stability of the biosphere and also including climate and the ecosphere, i.e., the Earth system as used by humans.

Complementarily, EcoRisk indicates where and when critical transitions in ecosystems occur (as a result of NPP appropriation, but also from climate change, or other environmental pollution). It is thus a very useful indicator to assess the risk of ecosystem destabilization today and also going forward, depending on which pathway humanity will take in the future. We deliberately call it a risk metric because the mathematical property of measuring a non-directional change for us directly translates to a proxy for the risk of ecosystem destabilization. Of course there could be regions where ecosystems can benefit from changing biogeochemical properties; however, we would argue that when comparing to a long-term stable, land-use-free preindustrial reference situation (like the Holocene), most ecosystems would be at the risk of being thrown out of this equilibrium when presented with changing conditions. This choice of reference might need to be changed for application of EcoRisk to future climate stabilization scenarios with “Earth-System Stewardship” (Rockström et al., 2021; Steffen et al., 2018).

5 Conclusions

Humanity faces the challenge of stabilizing the Earth system in the Anthropocene. This requires stabilizing the climate and maintaining a resilient biosphere, which represents, according to Steffen et al. (2015), a key pillar of the Earth system. Here we present two computable biosphere integrity indicators that allow for the assessment of the historical and future risk of biosphere destabilization. These are designed to complement other biosphere integrity indicators.

Human appropriation of NPP is a direct result of land use changes, modulated by climatic changes, irrigation, fertilization, and other management, leading to biosphere integrity loss that can be measured by EcoRisk. On the other hand, natural NPP is also modified by changes in climate, water availability, biosphere integrity, and human land use. Thus, BioCol and EcoRisk are metrics that integrate multiple terrestrial planetary boundaries (Rockström et al., 2009; Steffen et al., 2015) and human drivers into two numbers, comparable to what global mean temperature does for the underlying more complex climate processes. Therefore, BioCol and EcoRisk have the potential to be included as indicators in an updated planetary boundaries framework. While well established as a concept, so far there has been no simulation model for assessing and projecting HANPP on the global scale for example to analyze future scenarios not accessible through inventory data.

Both BioCol and EcoRisk are spatially explicit, process-based, and computable metrics that can be aggregated over space and time. They can therefore serve as integrative meta-level proxies for biosphere integrity and Earth system stability and their dynamic changes. The code used here for their calculation is distributed as an open-source R package, and we invite external contributions. With the next update, we plan to equip it to deal with netcdf files, a common data exchange format for spatial data, which allows for future application to different models. Thereby, we hope to encourage others to join our quest to better understand the role of ecosystems for Earth system stability and find ways to preserve it.

Appendix A: Abbreviations

BFT	Bioenergy functional type
CFT	Crop functional type
GCT	Ground cover type
HANPP	Human appropriation of NPP
LPJmL	Lund–Potsdam–Jena managed Land
	(a dynamic global vegetation model)
BioCol	Metric for human biomass colonization pressure
EcoRisk	Metric for risk of ecosystem destabilization
NPP	Net primary production
PFT	Plant functional type

Appendix B

Table B1 Names and abbreviations for all biomes according to our biome classification.

Figure B1 (a) Average land use between 1988–2017 in percentage of grid cell area. Total historical sum of yearly land use fractions [0–1] with an (b) exponential and (c) linear legend.

Figure B2 (a) Change in biochemical compositions computed by EcoRisk (as in Fig. 3, but excluding nitrogen variables) between 1550–1579 and 1985–2016. (b–e) EcoRisk components, i.e., vegetation structure change, local change, global importance, and ecosystem balance.

Figure B3 Biome internal difference distribution for the forest biomes computed using EcoRisk between each cell of the biome with the average cell of the biome (for 1550–1579 states). Note that only very few cells are classified as BoND or TeBE; therefore, their distribution appears very homogeneous. For biome names and abbreviations, see Table B1.

Figure B4 Comparison of the state difference (EcoRisk) between biomes according to the average biome cell (y: reference state; x: scenario state) evaluated for the average state of 1550–1579. The color code is the same as in Fig. B3. Since only very few cells are classified as TeBE, this biome is left out. Biomes RoIce and Wat do not host vegetation and are also left out. For biome names and abbreviations, see Table B1.

Figure B5 Change in (a) leaching and (b) nitrogen volatilization between 1550–1579 and 1985–2016.

Figure B6 Change in biochemical compositions computed by EcoRisk between 1550–1579 and 1985–2016. (a) Total EcoRisk and an evaluation of $(\text{gi}+\text{lc}+\text{eb})/\mathrm{3}$ for only components related to (b) carbon, (c) water, and (d) nitrogen.

Figure B7 Cumulative area for all cells defined as not similar (EcoRisk > threshold) according to Fig. B4 (green) and defined as similar (EcoRisk < threshold) to the biome average according to Fig. B3 (blue). The value where the sum of both curves is maximal (an optimal threshold value fulfilling both criteria) is indicated by the dashed vertical line.