Climate change critically affects the status of the land-system change planetary boundary

The planetary boundaries framework defines a safe operating space for humanity. To date, these boundaries have mostly been investigated separately, and it is unclear whether breaching one boundary can lead to the transgression of another. By employing a dynamic global vegetation model, we systematically simulate the strength and direction of the effects of different transgression levels of the climate change boundary (using climate output from ten phase 6 of the Coupled Model Intercomparison Project models for CO2 levels ranging from 350 ppm to 1000 ppm). We focus on climate change-induced shifts of Earth’s major forest biomes, the control variable for the land-system change boundary, both by the end of this century and, to account for the long-term legacy effect, by the end of the millennium. Our simulations show that while staying within the 350 ppm climate change boundary co-stabilizes the land-system change boundary, breaching it (>450 ppm) leads to critical transgression of the latter, with greater severity the higher the ppm level rises and the more time passes. Specifically, this involves a poleward treeline shift, boreal forest dieback (nearly completely within its current area under extreme climate scenarios), competitive expansion of temperate forest into today’s boreal zone, and a slight tropical forest extension. These interacting changes also affect other planetary boundaries (freshwater change and biosphere integrity) and provide feedback to the climate change boundary itself. Our quantitative process-based study highlights the need for interactions to be studied for a systemic operationalization of the planetary boundaries framework.


Introduction
The planetary boundaries framework specifies the boundaries of relatively stable planetary conditions for nine human-perturbed Earth system processes, each with one or more control variables.A severe or persistent boundary transgression signifies entering a zone of increasing risk and subsequently highrisk where Earth's biophysical functionality and selfregulatory capacity are threatened (Rockström et al 2009a, 2009b, Steffen et al 2015, Richardson et al 2023).Human activities have already led to the transgression of six of the nine boundaries: climate change, land-system change, biosphere integrity loss, biogeochemical flows (Steffen et al 2015), release of novel entities (Persson et al 2022) and freshwater change (Wang-Erlandsson et al 2022).Despite the framework's intrinsic systemic perspective and the potentially synergistic nature of boundary transgressions, the transgression level for a given planetary boundary is usually quantified without factoring in the altered status of other planetary boundaries.There is a dearth of research on the current and potential future implications of planetary boundary transgressions and, notably, a lack of process-based understanding of interactions among them.Previous studies have either simply acknowledged the theoretical existence of planetary boundary interactions (Rockström et al 2009b, Steffen et al 2015), or explored interactions semi-quantitatively either based on globally aggregated values (Lade et al 2020) or expert elicitation (Chrysafi et al 2022).In an initial effort to develop a process-based understanding of interactions, Lade et al (2021) studied changes in runoff, vegetation cover and climate (as proxies for planetary boundaries) using the dynamic global vegetation model (DGVM) LPJmL to assess cross-scale interaction strengths as a basis for a prototype Earth system impact metric for use in corporate sustainability reporting.This incomplete understanding of boundary interactions contributes to misperceptions of the framework itself, as reflected in the selective operationalization of single boundaries in policies (Keppner et al 2020).Without a substantive and communicable understanding of their interactions, it remains challenging for stakeholders to understand the importance of considering the boundaries together, let alone identify and prioritize actions with positive impacts across boundaries (Arup 2021).Given the growing uptake of the planetary boundaries framework (Downing et al 2019), there is an urgent need to develop robust methods to systematically investigate a wide range of boundary interactions that quantitatively underpin the framework and its operationalization.
The underlying systemic perspective of the planetary boundaries framework must be made explicit by analyzing the feedbacks among Earth system processes that define individual planetary boundaries.Here, we begin with a focus on climate change and land-system change.The planetary boundary for climate change is crucial for maintaining Earth system stability, because of the strong direct influence of climate on other boundaries (Steffen et al 2015, Lade et al 2020).Anthropogenic greenhouse gas emissions have already raised CO 2 concentrations to an unprecedented level in the last 23 m.y.(Cui et al 2020), far exceeding 350 ppm CO 2 -the planetary boundary for climate change (Rockström et al 2009b).Consequently, global warming may have already set several tipping points in motion (Lenton et al 2019, McKay et al 2022), increasing the risk of the emergence of an undesirable state of the Earth system (Steffen et al 2018).
The land-system change boundary quantifies the remaining natural forest area essential for sustaining critical Earth system functions (Steffen et al 2015).Forests constitute a dominant terrestrial ecosystem (Pan et al 2013) and play a crucial role in Earth system functions (Bonan 2008) operating across local, regional and continental scales (Ellison et al 2017).These processes include sustaining the hydrological cycle through evapotranspiration provision (van der Ent et al 2010), surface cooling (Zhang et al 2020), carbon sequestration (Sohngen and Mendelsohn 2003), biogeochemical cycling (Lang et al 2016), albedo changes (Lee et al 2011), biogenic volatile organic compound emissions (Unger 2014) and habitat provision (Brockerhoff et al 2017).The land-system boundary is quantified specifically for each of the three forest biomes-temperate, tropical and boreal-due to varying land-surface climate coupling importance (Snyder et al 2004, West et al 2011, Steffen et al 2015) and potential tipping points, such as a boreal forest dieback and Amazon forest decline (Lenton et al 2008(Lenton et al , 2019)).A description of the biome-specific boundary values and a map of the current status of each biome on each continent are shown in figure 2 The main objective of this study is to investigate how different degrees of transgression of the climate change boundary (i.e.different atmospheric CO 2 levels) induce shifts in the land-system change boundary (the remaining extent of natural forest biomes-tropical, temperate and boreal), thereby quantifying the one-way interaction between these two boundaries.Furthermore, we explore cascading effects on other planetary boundaries within the affected regions, considering the boundaries for freshwater change (root-zone soil moisture), biosphere integrity (habitat intactness) and feedbacks to the climate change boundary itself (through changes in albedo and net biome productivity).For this, we employ the process-based LPJmL DGVM forced by an ensemble of climate change projections of phase 6 of the Coupled Model Intercomparison Project (CMIP6), ordered to represent different levels of atmospheric CO 2 concentrations until the end of the century and the end of the millennium.

Methods
Owing to their capacity to simulate the key processes underlying the impacts of climate change on vegetation distribution and composition (growth, mortality and resource competition) and disturbances such as droughts and wildfires, DGVMs are a prime tool for studying the effects of different levels of climate change boundary transgressions on the land-system change boundary and on other terrestrial planetary boundaries.Therefore, to systematically assess these interactions, we employ the state-of-the-art DGVM LPJmL (von Bloh et al 2018) which has a long record of model improvements and validations (Sitch et al 2003, Gerten et al 2004, Bondeau et al 2007, von Bloh et al 2018, Schaphoff et al 2018a, 2018b) and has proven to be suitable for studying the planetary boundaries framework (Gerten et al 2013, 2020, Heck et al 2018a, 2018b).A model description is provided in the appendix.LPJmL is externally forced by climate input, here using high-emission scenario forcings following the representative concentration pathway 8.5 combined with the shared socioeconomic pathway 5 (hereinafter SSP5-RCP8.5) of a multi-model ensemble from ten different general circulation models (GCMs) from CMIP6 (Eyring et al 2016), to account for uncertainties between climate change projections.Anthropogenic land-use is held constant at the level of 2015 to isolate the climate effect.We analyze five climate scenarios for stabilization and transgression (figure 1).A CO 2 concentration of 350 ppm sets the planetary boundary for climate change, with 450 ppm demarcating the upper end of the increasing risk zone (Steffen et al 2015).The three other scenarios depict different magnitudes of boundary transgression within the high risk zone (modest, 550 ppm; strong, 750 ppm; extreme, 1000 ppm).To realize these scenarios, individual runs are branched off from the SSP5-RCP8.5trajectory once the scenario's CO 2 concentration has been reached.Following the 'warming slice' approach by Schleussner et al (2016), the 20 surrounding years are repeated and shuffled until the year 3000 (the simulation protocol is explained in the appendix).Outputs were subsequently analyzed for both the end of this century (year 2100) and-to investigate possible long-term effects of equilibrated CO 2 levels and corresponding climate forcing-the end of the millennium (year 3000) for each of the five climate scenarios.
Following the definition of the land-system change boundary (Steffen et al 2015, Richardson et al 2023), our analysis emphasizes climate-changeinduced area changes within the pristine forest extent under potential natural vegetation (PNV), which acts as the reference status.To model the pristine biome distribution for the three forest types under PNV, LPJmL was run with pre-industrial climate (period 1850-1859) and without human land-use (subsequent to the spin up period).For each GCM forcing, the expansion of forest biomes under each scenario as well as their extent under current and PNV conditions was calculated following a cell-based biome classification scheme (see Supplementary Information SI for a detailed description).Subsequently, the respective status of the landsystem change boundary for each ppm scenario and the two time periods was calculated by subtracting the total biome area under PNV by its remaining extent under a given scenario.
Here, we follow the biophysically mediated planetary boundary interaction scheme by Lade et al (2020) where alterations in pressures on one planetary boundary lead to an impact on another planetary boundary through a biophysical mechanism, causing changes in the respective transgression level.Knock-on effects of land system changes were subsequently analyzed by assessing shifts in the planetary boundaries of biosphere integrity, freshwater change and climate change, using simulatable proxies broadly aligned with the definitions of the planetary boundaries that also capture key mechanisms of systemic interaction: habitat intactness index (described in the SI), mean annual changes in plant available rootzone soil moisture based on Wang-Erlandsson et al (2022), and mean annual changes in albedo and net biome productivity (a measure of carbon accumulation related to the climate change mitigation ability of ecosystems, described in the SI).Shifts were assessed by comparing the current status (mean 2005-2014) of each proxy to the simulation period around year 3000 (mean 2996-3005) under the 750 ppm scenario in the affected regions.Our analysis spans the nexus of climate, land, ecosystems and water to assess boundary interactions, but a fully integrated Earth system model is required to perform a comprehensive analysis of individual interactions; however such models currently run at too coarse scales with a high computational burden.

Results
The reference status for the land-system change boundary (pristine PNV extent) of the three forest biomes is mapped in figure 2(a).The current boundary status is described as the remainder of the forested area.The color-coded boundary status of each biome on each continent is shown in figure 2(b), following Steffen et al (2015).We found the safe boundaries for boreal and tropical forests to be transgressed across all continents.Temperate forests, on the other hand, are still regionally below its boundary in our simulation, which is partially attributable to the lower boundary value for this biome (see legend figure 2(b)).Landuse and climate change are primarily responsible for the deterioration of forest biomes over the historical period, but future changes (as shown in figure 3) are driven by climate change alone (since land-use has been held constant at 2015 levels in our simulations).Note that figure 3 features global aggregates of the land-system change boundary transgression levels.

Boreal forest
The boreal forest biome shows the strongest simulated response to climate change in terms of boundary transgression levels (red zone, figure 3).Staying within the climate change boundary of 350 ppm maintains the biome's stability by keeping it at its current position, even in the long run (by 3000, as shown in figure 3).Further transgression of the climate change boundary above the current (mean 2005-2014) value is matched by further transgression of the boreal boundary (figure 3), and is likely to initiate a significant and non-linear loss of the southern fringe of the pre-industrial boreal forest (figure 5(a)), especially in scenarios of 550 ppm and higher.With increasing ppm levels and warming, the GCM inter-model projection uncertainty increases and a delayed climate-vegetation quasi-equilibrium is notable as shown by the expanding interquartile range and the stronger boundary transgression by year 3000 than by 2100, respectively (figure 3).Very high ppm scenarios even result in a near-complete dieback of pristine boreal forest, especially by the end of the millennium.The remaining forest cover extent is reduced to a median of 37.8% (min = 8.2%) and 14.5% (min = 0.5%) under GCM forcings of 750 ppm and 1000 ppm, respectively, by the end of the millennium (see below for parallel boreal forest expansion outside the PNV area).

Temperate forest
In contrast to the boreal forest, temperate forests demonstrate a higher resilience toward climate shifts, as our simulations do not indicate a major forest loss and, on the contrary, even a short-term stabilization under any climate scenario.Respecting the safe climate operating space (350 ppm) would preserve the temperate forest biome extent in the near and long term (figure 3), but long-term transgressions of the climate change boundary beyond 550 ppm lead to a reduction of pristine temperate forest cover extent, pushing the land-system change boundary for this biome into increasing risk territory (yellow zone, figure 3).However, climate change produces a significant expansion of temperate forest beyond the PNV area (figure 5(b)).

Tropical forest
Of the three major forest biomes, tropical forests are simulated to experience the smallest climate-caused change in area extent.That is, across all climate scenarios and GCM forcings considered, the status of the land-system change boundary of the tropical forest remains close to its current value (solidly in the zone of increasing risk but not entering the high-risk zone; figure 3).

Spatial patterns of biome shifts
Changes in forest cover within the PNV extent determine the position of the land-system change boundary.However, under many scenarios, new climatic conditions push forest ecosystems beyond the periphery of their PNV extent, mostly through poleward shifts of the three major forest biomes (figure 4(a)).While the boreal forest biome reaches its largest expansion under 550 ppm (including the area reaching beyond the PNV area), the overall area consistently declines under higher ppm levels; in particular the pristine area (formerly) covered by PNV faces a near-complete dieback (inner circle, figure 4(b)).Boreal forests dislocate to adequate ecoclimatic zones, allowing temperate forests to establish in formerly boreal forest territory as seen in figures 5(a) and (b) in the midlatitudes and in the opposing latitude sums in figure 4(a) (additional maps in the SI).In contrast to boreal and temperate forest biomes, tropical forests are found to experience only small shifts and modest net gains in biome cover.

Biosphere integrity boundary
Forest habitat intactness changes are illustrated for the 750 ppm climate change scenario in figure 5(c).Different levels of change (10%, 30% and 50%) are displayed for forested areas whose biosphere integrity index (BII) is currently within the planetary boundary for terrestrial biodiversity (Newbold et al 2016).We find that most forested regions are subject to strong vegetation structural change.In particular, the boreal dieback and subsequent temperate forest succession as well as the poleward-moving treeline constitute large BII changes in the boreal and adjacent tundra zone.Additionally, the tropical forests witness a shift from evergreen to deciduous tree types (see part J of the SI).

Freshwater change boundary
Our results imply an overall increase in annual green water availability at high latitudes, where boreal forest is simulated to expand northward (figure 5(d) for the 750 ppm scenario).In contrast, decreases are projected along the southern boreal fringe (zone of boreal forest loss), Europe (temperate forest recession) and the Amazon rainforest.

Mutual feedback between climate change and land-system change boundaries
Our results attest a significant reduction in albedo resulting from the poleward boreal extension (figure 5(e)).While the area of temperate succession experienced a minor reduction in albedo, the tropical forest biome does not witness larger shifts.Net biome productivity (ecosystem carbon fluxes, including wildfires), in turn, experiences the most pronounced weakening in the tropical zone and minor decreases in the boreal and temperate zones (figure 5(f)).

Discussion
Our simulations reveal the direction and strength of change for the investigated planetary boundary interactions.That is, more or less sizable expansions or contractions of boreal, temperate and tropical forest biomes both within and beyond their respective PNV  (Steffen et al 2015).They delineate the safe operating space for the forest cover extent (green areas to the left) from the increasing risk and high risk zones (yellow and red areas to the right).The x-axis indicates the scenario forest extent relative to PNV, theoretically ranging from 100% (biome extent as under PNV) to 0% (no PNV forest cover remaining).To improve the readability, the x-axis was adjusted for each biome.The interquartile ranges capture the variations of the biome-specific land-system change boundary status for the ten GCM forcings for both the year 2100 and 3000, with the circle marking the median.
areas, are likely to occur in response to different levels of climate change boundary transgressions.The effects are long term, and the degree of disequilibrium (delta in biome shift between the years 2100 and 3000 as can be seen in figure 3) is largest for the strongest climate transgression scenarios in temperate and especially boreal forest biomes.This legacy effect of climate change calls for further long term environmental equilibrium studies as presented here.
In the following, we discuss these effects in more detail and put them into the larger context of overall planetary boundary interactions.

Boreal forest
High latitudes host one of the largest and most intact biomes, the boreal forests.The planetary boundary for the boreal forest biome, set to a precautionary 85% remainder of PNV area, is here found to be transgressed in both North America and Eurasia, which is in line with the assessments of Steffen et al (2015) and Richardson et al (2023).Globally, approximately 70% of the original extent remains in our simulations (figure 3).This biome is disproportionately affected by further transgression of the climate change boundary, perhaps because of polar warming amplification (Taylor et al 2013, Thomas et al 2020).Climate change-induced warmer and longer growing seasons benefit net primary production and atmospheric CO 2 sequestration, although the magnitude of this negative feedback is disputed (Zohner et al 2020).Higher temperatures also accelerate soil decomposition, permafrost thawing, insect dynamics and wildfires, leading to complex response patterns of boreal forests to global warming (Zhang et al 2013).
Accordingly, we find that the impacts of climate change are far reaching.On the one hand, global warming initiates a strong poleward shift and woody invasion of the adjacent tundra, extending the tree line and thus the boreal forest further north (i.e.forcing the boreal forest beyond the borders of its PNV area), thereby reinforcing dramatic changes to the Arctic and subarctic environment that are already observed today.On the other hand, the boreal forest biome is subject to a large-scale and non-linear forest loss in the southern fringe of the ecosystem in our simulations (>450 ppm).This is in line with earlier  et al 2022).This uncertainty is further exacerbated by the choice of mortality processes included in a DGVM (Pugh et al (2020) and particularly table 3 in their paper providing an overview on mortality processes).In addition to background mortality, LPJmL incorporates stress-related mortality arising from competition, fire disturbances and heat stress related mortality (exclusively affecting boreal trees).This heat stress-induced mortality is determined by quantifying the number of days surpassing a temperature threshold (23 • C) and associating it with a damage heat function (see table S4 in Schaphoff et al (2018b)).We found the emergence of the boreal to temperforest transition in our simulations to be primarily attributable to this heat stress-induced mortality, causing a dieback of southern boreal forest (section K in the SI).The inclusion of boreal heat stress mortality in LPJmL is potentially elucidating the lack of similar findings in comparable studies (Gonzalez et al 2010, Steinkamp et al 2015).Ito et al (2016) examined simulations from seven biome models under RCP8.5 and found that among them, LPJmL is the only DGVM that exhibits a decrease in vegetation carbon at lower high altitudes, likely due to the aforementioned heat stress function.In our simulations, less extreme climate change scenarios favor a poleward expansion over a southern dieback (figure 4(a)), but both the 750 ppm and 1000 ppm scenarios reverse this pattern in that southern boreal forest loss outweighs gains in the boreal north.This detrimental development is even more pronounced when only considering the locations covered by boreal forest under PNV (in accordance with the land-system change boundary).

Temperate forest
Among all forest biomes, temperate forests are the most severely affected (in our simulations, only 51% of the total PNV extent remains) due to agricultural expansions, especially in North America and Europe where the temperate forest biome is at increasing risk (figure 2(b)).Regarding future climate change effects, however, our simulations indicate quite high resilience with relatively little loss within their PNV areas, even under the most extreme climate scenarios.While this appears to contradict Sitch et al (2008), who found a loss of temperate forests based on simulations with LPJmL's predecessor model (Sitch et al 2003), this divergence is likely to originate from different climate forcing (and possibly from differences in processes in the models).This stability in terms of spatial extent (in PNV areas) needs to be evaluated in light of the altered disturbance regimes already observed under current climate change in Europe (Seidl et al 2014) and projected to increasingly threaten forest health (Seidl et al 2017)-with disturbances, such as severe droughts and wildfires or pathogens and insect outbreaks, becoming ever more likely under future climate change (Millar and Stephenson 2015).However, extended growing seasons (for example in Europe, Menzel et al (2020)) and the accompanying prolonged canopy duration result in an elongated carbon uptake period of forests, thereby increasing their net ecosystem production (Banbury Morgan et al 2021).This trend is predicted to continue into the future, supporting tree growth and productivity of some species (Vitasse et al 2011).
Moreover, the retreat of the southern boreal forest appears to benefit the northward-shift in temperate forest biome in outcompeting current species, i.e. there is a significant expansion of temperate forest outside the PNV area especially in the long term (for maps refer to section H in the SI).Besides the heat-stress related boreal dieback discussed above, the changing bioclimatic conditions are exceeding the chilling requirement for the establishment of boreal tree-type plant functional types (PFTs) (section K in the SI) benefiting the replacement by temperate forest.The complex role of wildfire in this biome shift is beyond the scope of this study.

Tropical forest
Tropical forests are hotspots of biodiversity (Gardner et al 2009), cool the planet due to their high evaporation rates, hold the largest amount of aboveground biomass among all forests (Pan et al 2011) and figure about one third of global terrestrial primary production (Beer et al 2010).At the same time, they constitute the nexus of global land-use change, experiencing the highest rates of ongoing deforestation among all forest biomes (Hansen et al 2013).According to our cropland data and PNV simulation, only 69.5% of the PNV extent remains today.Forest clearing has thus already transgressed the safe boundary and, if deforestation is not brought to a halt, is moving it close to high risk territory (>60% remaining; red area in figure 3).At the same time, projected future warming (and concomitant changes in precipitation, wildfire regimes and atmospheric evaporation demand) question the stability of the remaining tropical forests (Mitchard 2018).
The Amazonian tropical rain forest is of particular interest as its potential dieback under future climate change has been identified as a tipping element central to Earth system stability (Lenton et al 2008, McKay et al 2022).This assessment, however, refers to studies by Cox et al (2000), Cox et al (2004) who projected major forest loss under highly elevated CO 2 concentrations in the HadCM3 climate-carbon cycle model.This dieback is contested by Malhi et al (2009), who, when running a GCM ensemble approach of 19 members, report that the climate-induced intensified dry-season water stress is transforming parts of the Amazonian rainforest to seasonal forests, instead of causing a forest dieback and subsequent succession of savanna.Cox et al (2000), Cox et al (2004) were unable to simulate this distinction in structural forest features using HadCM3.In a recent study, Boulton et al (2017) reapplied HadCM3 with 57 model configurations (varying in parameters) and found that the Amazon rainforest is stable by the end of the century for most model configurations; however uncertainties increased for the highest emission scenarios.
In our simulations, the remainder of the tropical forest biome is remarkably resilient even to the strongest future warming levels (and associated precipitation declines).While a dieback of the tropical forest is not emergent in our scenarios, they still show a shift from evergreen to deciduous rainforest (for maps, see part J of the SI) which is in accordance with the transformation to seasonal forest reported by Malhi et al (2009).This important distinction of forest structure is lacking in the definition of the landsystem change boundary.Overall, the lack of a projected tropical forest dieback is consistent with earlier assessments, indicating that tropical forests may continue to act as a large terrestrial carbon sink even The possibility of moisture recycling feedbackinduced hysteresis behavior of the Amazon rainforest has been suggested (Staal et al 2020).This vital risk is not accounted for in our simulations since dynamic moisture recycling is not featured in LPJmL.Moreover, a simulation of the planned free-air CO 2 enrichment experiment in the tropics (AmazonFACE, https://amazonface.inpa.gov.br/;Fleischer et al (2019)) emphasize that models lacking phosphorus availability (such as LPJmL) highly overestimate the additional biomass growth through direct CO 2 effects, thereby distorting the resilience of tropical forests to climate change.Thus, the simulated stability and climate change resilience should be interpreted with caution.
Our results further show that global warming causes a moderate latitudinal expansion of tropical forest beyond the PNV extent which is consistent with the tropical belt expansion theory (Staten et al 2018) as well as modeling and reanalysis studies showing that changes in radiative forcing lead to shifts in the edge latitude of Hadley cells that define the tropics (Davis and Birner 2017).This result is also consistent with the widely observed phenomenon of shrub encroachment (Higgins et al 2000) which is linked to a savanna-forest transition over a wide range of land management practices as a response to elevated atmospheric CO 2 concentrations (Midgley and Bond 2015).Notwithstanding these expansions occurring beyond the PNV territory, the land-system change boundary for tropical forest remains transgressed due to historic and ongoing deforestation.

Further planetary boundary interactions
The above main analysis focuses on area changes (mostly reductions) within the PNV domains of the three major forest biomes, as the land-system change boundary defines limits to the loss of these areas.In addition to the extensions of forest biomes into regions outside their PNV areas (also shown above), we here provide a further analysis of knockon effects on other planetary boundaries associated with climate change-induced shifts in the landsystem change boundary (see figures 5(c)-(f)), using the control variables (or proxies thereof) as defined in the SI.
The habitat intactness proxy metric for the biosphere integrity boundary indicates extensive structural changes in ecosystems that are currently characterized by a high score in the BII.Climate-change induced relocation of biomes will strongly change ecosystem statuses worldwide (figure 5(c)), likely affecting how habitats can provide refugia to wildlife (Eigenbrod et al 2015).Future research is needed to incorporate a more complex metric of change such as gamma, an aggregated metric of combined structural and biogeochemical shifts in ecosystems (Ostberg et al 2013), to study shifts and disruptions and thus interactions between the climate change and biosphere integrity boundaries.
Our analysis of green water shows that LPJmL projects spatially heterogeneous changes in plant available soil moisture in the multiannual mean (figure 5(d)).In a recent study, Wang-Erlandsson et al (2022) highlighted the importance of green water in sustaining crucial Earth system functions like carbon sequestration, moisture recycling and evaporative cooling.Future studies will have to disentangle seasonal patterns and identify more systematically whether and where wet and dry departures lead to a (further) crossing of the planetary boundary for green water.Additional analysis is needed to link green water changes to either climate change directly or to the land cover changes resulting from it.Whether shifts in green water are responsible for process chains, such as boreal dieback and changes in net biome productivity, both emergent in our studies, is of particular interest.
Climate change-driven forest redistributions and tree line advancements may play an important role in climate regulation by changing the surface albedo (figure 5(e)) and the net carbon flux (figure 5(f)) to the atmosphere since the surface albedo feedback is characterized e.g. by dark forests reflecting less and absorbing more incoming solar radiation than tundra vegetation, thereby amplifying polar warming (Zhang et al 2013).Comprehensive Earth system models capable of dynamically simulating how changes in carbon flux and albedo modulate the Earth's radiative balance need to be employed in future studies to assess the strengths of the feedbacks back to the climate change boundary found here.

Conclusion
Our simulations of the effects of the transgression of one boundary on the transgression level of another in a systematic and process-based manner constitute a quantitative advancement of the planetary boundary framework, improving its usability by stakeholders concerned about the outlook for an Earth system in which multiple boundaries are being transgressed.We show that transgressed boundaries can be simulated dynamically and mapped, taking key knock-on effects on the climate, land, water and ecosystem nexus into consideration.We account for inter-model uncertainty in the latest climate projections and multi-centennial delays in the Earth system, thereby overcoming weaknesses in earlier planetary boundary interaction assessments (Lade et al 2020(Lade et al , 2021)).
Moreover, we demonstrate that the existing pressures on the land-system change boundary are exacerbated when the climate change boundary is further breached, although the effect is not equally felt across biomes.The spatial extent of tropical forests is characterized by a high degree of resilience toward future warming in our simulations, whereas temperate and especially boreal forest biomes are subject to ever stronger and increasingly non-linear change in their PNV extent, once climate change and the accompanying precipitation shifts go beyond the vegetation type-specific thresholds and inter-species competition is causing their replacement (figure 3).In particular, when reaching the high risk zone of the climate change boundary (>450 ppm), one-and possibly two-way interactions between these two boundaries exacerbate the changes in biome extent.The possibility that these new dynamics could result in the establishment of alternative stable states is a question that requires further research.Respecting the planetary boundary for climate change (350 ppm) would ceteris paribus stabilize all forest biomes, upholding them close to their current position, even by the end of the millennium.This, together with the non-linear biome shifts occurring beyond 450 ppm (the 'high risk' zone) indicates that the climate change boundary is well placed in relation to the land-system change boundary.
Our analysis also reveals that the current control variable of the land-system change boundary is not well-suited for resolving heterogeneous spatial patterns of climate change-induced shifts in forest biomes.For a more operational revision of this boundary, it would be crucial to account for biome shifts not only within but also outside the PNV area, and to consider further ecological, biophysical and biogeochemical changes that capture the stability and functioning of biomes (to the extent they are not covered by other planetary boundaries such as the one for biosphere integrity in particular).The brief analysis on additional pressures on a larger set of planetary boundaries highlighted here indicates the value of a comprehensive analysis of boundary interactions in future research, which in turn demands new directions in process-based Earth system modeling.In particular, model assumptions about the controls determining the distributions of PFTs have a substantial bearing on the spatial and temporal uncertainty in the analyzed interactions.Future studies of planetary boundary interactions could account for biosphere model-related uncertainties using an ensemble of multiple DGVMs that include a wider set of dynamic processes (e.g.insect disturbances, elaborated mortality and nutrient cycling) and that are parameterized and constrained using additional field measurements (Mevenkamp et al 2023).
The breadth of our findings flag that planetary boundary transgression levels change when interactions are considered, highlighting the need for further rigorous investigations of more combinations of boundary interactions.The resulting changes in transgression levels could imply that safe levels for some boundaries (e.g.biosphere integrity, land use change, freshwater change) may have to be tightened during periods of severe transgression of other boundaries (e.g.climate change).We stress that the power of the framework does not lie in the individual consideration of a single boundary, but in understanding the planetary boundaries as an interconnected and codependent stability landscape where political actions taken to mitigate pressures on one boundary affect Earth system stability in other planetary boundary dimensions (Gerten and Kummu 2021).To illustrate, large-scale plant-based carbon dioxide removal applications as a climate mitigation strategy have been found to worsen the status of other planetary boundaries (Heck et al 2018a).A richer understanding of interactions could reveal actions with the highest potential to mitigate pressure on various boundaries and inform policymakers about potential trade-offs.Improved systemic understanding and quantification of the planetary boundary framework is critical for stakeholders striving to identify what actions to take (or avoid) for positive impacts across boundaries, a missing piece for systemically integrated and coherent planetary boundary operationalization.
LPJmL.LPJmL is a state-of-the-art dynamic global vegetation model operating on a 0.5 • × 0.5 • grid spatial resolution and at daily time steps.The terrestrial biosphere is represented by 11 individually parameterized plant functional types (PFTs) of natural vegetation (for an overview see appendix S1) and 16 crop functional types (CFTs) to account for human land use.While the distribution of CFTs, through agriculture, is prescribed, the establishment and distribution of PFTs (excluding Antarctica) is simulated as the result of succession created by environmental conditions and inter-PFT competition for resources, disturbances such as wildfires, background mortality, and heat-stress mortality (for boreal trees).Spatiotemporal dynamics are the result of calculations of these parameters for each vegetation type in each grid cell at each time step.The fundamental underlying processes of the model are comprehensively described (Schaphoff et al 2018b) and successfully validated against local in situ, satellite derived and agricultural yield statistics data.This includes a comparison of simulated biome distribution to vegetation cover from remote sensing data (Schaphoff et al 2018a).Natural and human-induced wildfires are simulated via the process-based fire regime model SPITFIRE (Thonicke et al 2010).LPJmL5.1, the model version employed here, features the nitrogen cycle (von Bloh et al 2018), accounting for nitrogen shortage limiting plant production and constraining enhanced plant growth in response to elevated atmospheric CO 2 concentrations (Wang et al 2020).
Simulation protocol.LPJmL5.1 is forced by ten different GCMs from phase 6 of the Coupled Model Intercomparison Project, CMIP6 (Eyring et al 2016); for an overview of the GCMs see section G in the SI, and for analysis protocols and code see (Tobian et al 2024).This ensemble approach accounts for uncertainties between climate change projections.Initiated from zero, LPJmL was forced over an 8000 year spin-up period to establish a global PFT distribution and water and carbon reservoirs that are in quasi-equilibrium for each GCM.Subsequently LPJmL was run over both the historical period (1850-2014), a potential natural vegetation simulation (1850-1859) and the future projections for the high-emission scenario SSP5-RCP8.5 (2015-2100).A total of five scenarios was branched off from the future projection as described below.We chose SSP5-RCP8.5, the most extreme climate forcing of CMIP6 since it covers the largest range of potential future atmospheric CO 2 levels (i.e.different degrees of transgression of the climate change boundary).The daily climate data are bias-corrected and downscaled to 30 arc-min resolution to meet the requirements of LPJmL (Lange 2019).Human land-use (i.e.ecosystems converted to cropland and pasture) is responsible for the historic transgression of the land-system change boundary.In our simulation, historical land-use patterns (up until the year 2015) were prescribed (Frieler et al 2017), based on HYDE 3.2 (Goldewijk et al 2017).But land-use was held constant at the 2015 level in scenario simulations to isolate the climate change effect.Nitrogen limitation, a key nutrient constraining vegetation growth, was enabled.Deposition scenarios for both the historical period and the SSP5-RCP8.5scenario were obtained from the ISIMIP3 database (Jägermeyr et al 2021).
We analyze climate-stabilizing as well as overshoot scenarios ranging from a return to 350 ppm, a position at the 'high risk' zone (450 ppm), as well as a moderate (550 ppm), a strong (750 ppm) and an extreme transgression scenario of 1000 ppm (figure 1, main text).To realize these climate scenarios, individual runs are taken from the SSP5-RCP8.5trajectory once the target value of the desired atmospheric CO 2 concentration has been met.Following the 'warming slice' approach by (Schleussner et al 2016), each scenario run is forced with shuffled data from the warming slice of the 20 years surrounding the year when the scenario-specific atmospheric CO 2 concentration has been reached (see SI for warming slice years).This approach accounts for oscillations and thus allows to minimize interannual climate variability (Arguez and Vose 2011).The climate forcing of each scenario-and GCM-specific warming slice is repeated until the year 3000 (thereby passing the year 2100, figure 1).Simulations of long time horizons are necessary to reach near-equilibrium states in carbon pools (Schaphoff et al 2013) and to account for inertia in the response of the terrestrial biosphere to new climatic regimes (Jones et al 2009, Wu et al 2015).Vegetation adaptations were reported to take decades to manifest as Jones et al were able to show in a coupled climate-vegetation model (2009).Here, we follow an 'equilibrium vegetation' approach, where the climate is held constant after reaching the scenario-specific ppm level (figure 1) and the climate-induced biome shifts are studied both for the year 2100 and 3000.
(b).Agricultural practices such as burning, grazing and cultivation are the main drivers of landscape transformation and declining forest cover extent (Bhagwat et al 2014, Campbell et al 2017).Moreover, observations and modeling studies indicate first changes in vegetation distribution and composition, such as poleward shifts of the treeline in response to global warming (Parmesan and Yohe 2003, Lucht et al 2006, Boit et al 2016, Ostberg et al 2018, Boonman et al 2022).

Figure 1 .
Figure 1.Planetary boundary interaction simulation set-up.LPJmL5.1 is forced by SSP5-RCP8.5 scenario climate forcings provided by ten different GCMs from the CMIP6 ensemble (solid black line) succeeding the historical CO2 forcing (solid gray line).The blue dot represents the 'current' value for the year 2015 (398 ppm).Five ppm levels representing different statuses of the planetary boundary for climate change (ranging from 'safe' green to 'unsafe' red) were taken and held constant by recycling and shuffling a 20 yr window of surrounding climate and CO2 concentrations (colored rectangles).Changes resulting by the end of the century and millennium were assessed in the status of the land-system change boundary, and in the freshwater, climate change and biosphere integrity boundaries.In parallel to the historical forcing, a potential natural vegetation simulation was conducted to assess the undisturbed status of the land-system change boundary as the reference case.

Figure 2 .
Figure 2. Forest biome extent under potential natural vegetation, and the current status of the land-system change boundary.(a) Potential forest biome extent without land use and pre-industrial climate (mean 1850-1859), representing the reference status for the land-system change boundary.Biome classification based on Ostberg et al (2013) and described in SI.Only cells with GCM-forcing agreement of >30% are shown.(b) Current status (mean 2005-2014) of the land-system change boundary, based on the share of remaining forest biome in each continent of figure 2(a).Green areas are below the boundary; other areas have breached the planetary boundary and are subject to increasing risk (yellow) or high risk (red).Planetary boundary threshold values taken from (Steffen et al 2015, Richardson et al 2023).

Figure 3 .
Figure 3. Effects of transgressions of the climate change boundary on the status of the land-system change boundary.The bold vertical lines indicate biome-specific planetary boundaries(Steffen et al 2015).They delineate the safe operating space for the forest cover extent (green areas to the left) from the increasing risk and high risk zones (yellow and red areas to the right).The x-axis indicates the scenario forest extent relative to PNV, theoretically ranging from 100% (biome extent as under PNV) to 0% (no PNV forest cover remaining).To improve the readability, the x-axis was adjusted for each biome.The interquartile ranges capture the variations of the biome-specific land-system change boundary status for the ten GCM forcings for both the year 2100 and 3000, with the circle marking the median.

Figure 4 .
Figure 4. Simulated changes in forest biome area, ensemble median.(a) Latitude shifts of forest ecosystems by the end of the millennium, based on percentage change in biome area.While negative values indicate a loss, positive values describe a gain of the specific forest biome at the latitude under a particular scenario (irrespective of the PNV extent).(b) Biome areas under PNV, currently, as well as under different levels of the climate change boundary (by 3000).The size of the circles corresponds with the area of the biome (in Mha).The inner lighter circles indicate the remainder of the PNV biome and the darker color the biome area shifted beyond this extent.

Figure 5 .
Figure 5. Changes in the status of other planetary boundaries under the 750 ppm scenario between their current status (mean 2005-2014) and the simulation period around year 3000 (mean 2996-3005).Only changes in forested areas (under 750 ppm) are where >30% of the model ensemble is in agreement and where anthropogenic land-use constitutes less than 40% of a cell's area (in 2015).(a) & (b) Biome shift with red colors indicating a high GCM agreement of forest loss and blue colors a gain in forest for boreal forest (a) and temperate forest (b).(c) Habitat intactness change, depicting changes in plant composition in areas characterized by still integer ecosystems with a high biosphere integrity index (BII, Newbold et al 2016).(d) Changes in root-zone soil moisture, where red colors indicate a drying, blue a wetting trend (annual average).(e) Absolute surface albedo changes, where red colors indicate a reduction in albedo (darker surfaces).(f) Changes in net biome productivity, where red colors indicate a weakening of the terrestrial carbon sink.(c) Is chosen as a proxy for the biosphere integrity boundary, (d) represents the freshwater change boundary) while both (e) and (f) indicate biophysical feedbacks to the climate change boundary (see SI for details).
propositions (Lenton et al 2019) yet not seen in all vegetation models depending on which processes are dynamically represented (Friend et al 2014).Already, climate change causes more severe, frequent and intense boreal wildfires (de Groot et al 2013) and, since the dynamics of insect outbreaks and related increases in boreal tree mortality (Kurz et al 2008) are not incorporated in the LPJmL model used here, the emergent boreal dieback in our results may even be underestimated.Anderegg et al (2015) present a critique regarding the inadequate representation of drought-insect outbreak interactions in DGVMs, which are exacerbated by drought-induced stress in host trees under climate change.Instead, plant death is mostly simulated as 'background' mortality only, averaged over space and time (McDowell et al 2022), while short-term tree mortality mechanisms due e.g. to more intense drought conditions are not captured well (Allen et al 2015), underscoring uncertainties about future forest conditions (Hartmann under future climate change (Huntingford et al 2013, Schimel et al 2015, Fleischer et al 2019).However, this is in contrast to an earlier DGVM study conducted by Sitch et al (2008) and observation-based statistical model estimates which show a saturation and future drop in the carbon sink behavior in African and Amazonian tropical forests (Hubau et al 2020, Koch et al 2021).