Rise and fall of vegetation annual primary production resilience to climate variability projected by a large ensemble of Earth System Models’ simulations

Climate change is affecting natural ecosystems and society. Anticipating its impacts on vegetation resilience is critical to estimate the ecosystems’ response to global changes and the reliability of the related ecosystem services, to support mitigation actions, and to define proper adaptation plans. Here, we compute the Annual Production Resilience Indicator from gross primary production (GPP) data simulated by a large ensemble of state-of-the-art Earth System Models involved in the last Coupled Model Intercomparison Project (CMIP6) of the Intergovernmental Panel on Climate Change. In the Sustainability (Taking the Green Road) and Middle of the Road scenarios (ssp126 and ssp245), the areas where vegetation shows increasing GPP resilience are wider than the areas with decreasing resilience. The situation drastically reverses in the Fossil-fuel Development (Taking the Highway) scenario (ssp585). Among the larger countries, Brazil is exposed to the highest risk of experiencing years with anomalously low GPP, especially in the Taking the Highway scenario.


Introduction
The Sustainable Development Goals, formally embraced by the 2010 Conference of Parties, recognize the importance of ensuring conservation, restoration and sustainable use of terrestrial ecosystems and their services, and strengthening the resilience and adaptive capacity to climate-related hazards (SDG-15 and SDG-13, respectively;United Nations 2016). Stable ecosystems, characterized by small variations from their average state despite changes in environmental conditions, are indeed considered healthy and reliable in terms of the services they provide (MAE 2005, Costanza et al 2014. Ecosystems in good condition are necessary to secure the sustainability of human activities and human well-being (Maes et al 2020).
The concept of resilience is closely connected to ecosystem stability. Resilience has been defined either as the larger disturbance that a system can absorb without losing its structure, relationships and functionalities (Holling 1973, 1996, Walker et al 2004, Yi and Jackson 2021 or as the time required by an ecosystem to recover and return back to the equilibrium state after a disturbance (Pimm 1984, Yi andJackson 2021). These definitions-termed 'ecological resilience' and 'engineering resilience' , respectivelyare conceptually clear but do not directly provide a practical way to measure resilience (Morecroft et al 2012, Scheffer et al 2015. In fact, a quantitative estimation of resilience requires objective methods to identify and measure the external stresses and shocks (Meyer 2016). Also as a result of such indeterminacy, a large number of indicators was proposed to measure different aspects of resilience (de Keersmaecker et al 2014, Scheffer et al 2015, Meyer 2016, Yi and Jackson 2021. Up to date, none of these methods has been used to evaluate vegetation resilience at the global level and in future climate scenarios yet. Gross primary production (GPP)-the total carbon fixation by plants-is a primarily important terrestrial ecosystem function, at the point that it was also considered to be strongly related to resilience itself (Moore et al 1993, Stone et al 1996, Bryant et al 2019. Climate change is indeed expected to alter vegetation GPP resilience by potentially compromising the availability of water for vegetation in dry regions (Santini et al 2014, Zampieri et al 2019 and in general by increasing the frequency, amplitude and duration of extreme events that are detrimental for vegetation productivity (Dosio et al 2018, Naumann et al 2018. At the same time, the increase of atmospheric CO 2 concentration coming along with global warming is expected to bring positive effects in terms of vegetation photosynthetic rate (although acclimation should be also considered) i.e. the so called 'CO 2 fertilization effect' (Sage and Kubian 2007) and water use efficiency (Peters et al 2018).
Here, we used the Annual Production Resilience Indicator (R p ), defined as the squared mean annual GPP divided by its variance, which was recently proposed for a statistical evaluation of the production resilience of natural vegetation (Zampieri et al 2019) and agricultural systems from annual production time-series (Zampieri et al 2020b). R p is a simple but powerful indicator with several interesting properties. Being inspired by the ecological definition of resilience, R p is proportional to the amplitude of the largest disturbances that the system can absorb (measure by their rareness) and it is potentially consistent with the engineering definition (i.e. the return timing) as well (Zampieri 2021). It increases with diversity (number of species) and it accounts for memory effects, i.e. for perturbation recovery timings longer than a season (Zampieri et al 2020b).
We compute the Annual Production Resilience Indicator from an ensemble of 480 Earth System Models (ESMs) simulations included in the Sixth Coupled Model Intercomparison Project (CMIP6, www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6) and involved in the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, www.ipcc.ch/report/ar6/wg1/). ESMs are global climate models with an explicit representation of carbon processes and cycling over land, atmosphere and the oceans (Dahan 2010, Randall et al 2019, allowing to explore future climate variability of vegetation GPP according to different greenhouse gases emission scenarios. We quantify the relative changes in resilience of the GPP production with respect to period 1985-2015 for the near and the mid-long terms (2021-2050 and 2051-2100) under three scenarios of socio-economical global changes, corresponding to different levels of greenhouse gases emissions and land-use (i.e. the Sustainability, the Middle of the Road, and the Fossil-fuel Development scenarios) and we compute country level statistics of the larger projected changes.

Data and methods
The Annual Production Resilience Indicator (R p ), is based on two assumptions (Zampieri 2021). For annual producing systems such as agriculture or natural ecosystems that are sufficiently adapted to the environmental conditions and to the local climate, it is sensible to assume that the largest disturbances are rarer compared to the 'normal' conditions (assumption 1) and that the largest disturbances result in larger impacts of the annual production values (assumption 2). Under such conditions, the size of the disturbance can be univocally measured by its rareness e.g. the return period of production anomalies (T * ). Focusing on the production function only, the ecological resilience can be then simply measured by T * MAX , which is the return period of the largest adverse event that the system can cope with before losing completely the production ability. This approach is not sensitive to changes in composition and structure of the ecosystems, so it may allow for adaptation according to a more 'modern' interpretation of ecological resilience (Walker et al 2004).
For homogeneous production systems, it can be demonstrated that T * MAX is proportional to the annual production resilience indicator, defined as: where µ is the mean and σ is the standard deviation of the annual production (see appendix A available online at stacks.iop.org/ERL/16/105001/mmedia). In case the annual production resilience indicator is evaluated over a region including bare ground, the indicator is sensitive to the vegetated portion only (Zampieri et al 2019). It is interesting and potentially useful to disentangle the effects of changes in the mean and in the variability of GPP on the annual production resilience indicator. This can be accomplished by approximating the R P change with a first order 2D Taylor expansion of equation (1) as a function of the changes in the mean and in the standard deviation of GPP as follows: where the s stands for 'scenario' , the h stands for 'historical' , and ∆ represents the difference between two periods: where ∂ is the partial derivative. By computing the derivatives and dividing both members of equation (3) by R p one obtains:  Thus, an indication on the changes induced by the mean and the variability on the production resilience may be obtained by comparing the projected relative changes of the mean and of the variability, using the same weight. Equation (4) provides a normalized indicator of such comparison: which varies from −1 (variability dominates) to +1 (mean dominates), which is useful to assess and compare the dominant relative changes in different locations.
In this study, R P and its changes are computed on a large ensemble composed of all the climate change simulations for vegetation GPP available from all the Earth System Grid Federation (ESGF) portals up to 31 December 2019. The full list of simulations is provided in table 1, along with the detailed reference to the land surface component of the ESMs.
The ESMs land surface components include a prognostic representation of the biosphere with spatially distributed vegetation processes such as evapotranspiration, photosynthesis, carbon allocation and growth of leaves, stems and roots interacting with near surface meteorological variables such as temperature, radiation, wind and CO 2 concentration, and soil variables such as moisture, carbon and nitrogen (see references in table 1). Therefore, GPP variability in ESMs is the result of both bio-geophysical and biogeochemical processes such as soil moisture dynamics and energy budget, permafrost thawing, atmospheric CO 2 fertilization and nitrogen limitation as well as land use changes defined as a function of the different future scenarios. Some of the models also include a simplified representation of abiotic stresses on vegetation, such as fires (Lawrence et al 2019, Bastos et al 2020.
Memory effects linked to antecedent drought conditions are well reproduced since soil moisture dynamics and the related physical feedbacks were quite well developed already in the GCMs (Seneviratne et al 2010), which are the predecessors of the ESMs and from which ESMs inherits the representation of abiotic processes. In general, ESMs provide a reasonable representation of the GPP response to drought (see citations in table 1), which is, however, largely variable among models (Knauer et al . This motivates the use of a large ensemble for a robust assessment of GPP changes such as the one used here. The annual GPP is derived by summing up the monthly GPP outputs for each simulation listed in table 1 consistently with the spatial variability of vegetation seasonality (Peano et al 2019). The annual GPP data of each simulation is interpolated on a common 0.5 degrees regular grid with a second conservative remapping method (Jones 1999, Chen andKnutson 2008). The simulations' ensemble mean R p is computed for each ESMs and for each period and scenario. Finally, the overall median of the R P changes with respect to the historical period is computed for each future period and scenario. The robustness of the results is assessed by highlighting the areas where at least 75% of the models agree on the sign of changes.

Results
Future climate projections display significant changes of GPP variability resilience (figure 1) compared to period 1985-2014. The annual vegetation primary production resilience indicator is anticipated to generally increase in the lower emission scenarios (ssp126 and ssp245, figures 1(a)-(d), respectively). The larger positive changes are expected to occur especially in the snow dominated bioclimatic regions (see table S2 (available online at stacks.iop.org/ERL/16/105001/ mmedia)). The amplitude and the area covered by these changes are comparatively larger in the ssp245 scenario than in the ssp126 scenario and increase with time towards the end of the 21st century (table 2). Positive changes are also estimated for Central Africa and the Sahel regions, India and over the Himalayan Plateau. However, regions with loss of GPP resilience Table 2. Fraction of global land area where the relative resilience indicator (∆Rp/Rp) exceeds different thresholds (5%; 10%; 15%; 20%) in the simulation ensemble median. The first estimation (third and fourth columns) considers all areas displaying changes larger than the thresholds. A second estimate (fifth and sixth columns) is restricted to the areas where at least the 75% of the models agree on the sign of changes. appear as well, especially in Brazil, China and surrounding countries of equatorial America. Under the ssp245 socio-economic scenario, the CMIP6 ESMs project resilience losses also in Mexico and the southern part of the US, the Mediterranean region, Southern Africa and Australia. This occurs not only in the mid-long term (2051-2100, figure 1(d)), but also in the near term (2021-2050; figure 1(c)). Under moderate emission scenarios (ssp126 and ssp245), about one third of land area is going to experience an increase of vegetation annual GPP resilience over the period 2021-2050 (see table 2). This proportion is slightly lower, about one forth, when considering only the areas where 75% of the models agree on the sign of changes. Differently, the area with positive changes will cover almost half of the global land area (less than one third when the constraint on models' agreement is introduced) over the period 2051-2100. Regions losing resilience cover a smaller percentage of the global area, about 10% under ssp245 in the near term. The differences between the plain estimate and the one constrained on models' agreement become negligible for variations of resilience higher than 15% (table 2).
The results for the ssp585 scenario stand out significantly compared to the lower emission scenarios. Broad areas with negative change (i.e. loss of vegetation GPP resilience) appear already in the period 2021-2050 (figure 1(e)) in the Amazon region, the Unites States, South Canada, Western Europe, the Mediterranean basin, as well as in the Middle East, Central, Western and Southern Africa, Southeast Asia and China, and Oceania. Areas with at least 5% loss of vegetation GPP resilience are projected to cover approximately one fifth of the global land (12% considering models' agreement); while areas with more than 15% losses are projected to be limited to 3%. Positive gains of vegetation GPP resilience in boreal regions are simulated to be more limited with respect to the lower emission scenarios. Gains of at least 5% are expected to cover about one fourth of the global areas (15% considering models' agreement); while areas with changes larger than 15% are limited to the 6%, similarly to the ssp126 scenario.
The severity of the projected losses is expected to further exacerbate in the period 2051-2100. In the Taking the Highway scenario, less and less regions are expected to experience gains in vegetation GPP resilience. These regions are: La Plata basin in Argentina, part of the Sahel region, Eastern Africa, Western India, North-western China and some regions along the coast of the Arctic Sea. In general, areas gaining at least 5% resilience are simulated to be limited to 14% (8% considering model agreement) of global areas, while regions with more than 15% increase of vegetation primary production resilience are limited globally to 6% of the land area. The areas losing resilience are expected to outbalance those ones increasing resilience and cover 43% (25% considering models agreement) of global land area with more than 5% resilience losses. Globally, 13% of land areas are predicted to lose more than 15% vegetation primary production resilience.
The GPP resilience changes can be determined either by the change in the GPP mean and by the changes in the GPP variability due to climate change (see section 2). Positive resilience changes in the near term under moderate emission scenarios are often linked to positive changes in the mean annual GPP (figures 2(a)-(e) and S4) connected to overall higher levels of atmospheric CO 2 concentration and to higher mean growing temperature in Boreal Regions. Negative resilience changes are generally associated to an increase in the interannual variability of GPP (see figure S5). The areas affected by an increase of variability largely change across the scenarios and reach their maximum extent under the scenario ssp585 (figures 2(e), (f), S5(e) and (f)).
Gain and losses of resilience are quantified at the national level in order to provide country-specific information for adaptation options, and possibly to support ambitious mitigation policies. This analysis is displayed in figure 3 for the ten largest countries (and in table S3 for all World countries). Russia is characterized by the widest gains of resilience, which could cover almost 70% of the country area in period 2051-2100 under the ssp245 scenario. The spatial extent of areas expected to experience gains is reduced up to about 15% in the near term under the ssp585 scenario. This tendency continues towards the end of the century, under the ssp585 scenario, when also areas with GPP resilience losses start to appear. Canada shows a similar picture, but with less optimistic estimation of predicted losses largely outbalancing the gains in the 2051-2100 period under the ssp585 scenario. The US and China display similar figures, with gains predicted to reach 20% in the low emission scenarios (ssp126 and ssp245) and losses ranging from 10% to 15% in the ssp585 over period 2051-2100. Among the largest countries, Brazil is the one characterized by the most alarming projections, with the risk of losing resilience in 50% of its total territory under the ssp585 scenario at the end of the 21st Century. It is worth noting that these changes are likely to represent an underestimation as the current trend of land-use change (Freitas et al 2018) is only partially considered in the ESMs (Hurtt et al 2020). Australia is estimated to undergo negligible losses . Nevertheless, Australia will experience comparatively large losses of resilience towards the end of the century under the ssp585 scenario. The European Union is characterized by a more stable situation with significant positive changes only under the ssp245 scenario over the period 2051-2100. India shows projections similar to the EU, but with significant areas of vegetation that both gain and lose resilience under the high emission scenario at the end of the 21st Century. Large positive and negative changes in resilience are also estimated for Argentina under the high emission scenario (2051-2100). Whether or not these compensating changes in different area are beneficial for the countries' adaptive capacity could be subject of specific follow-up investigations.
Results for the remaining World countries (see table S3) allow identifying severe cases, such as losses of resilience higher than 15% over more than 50% of the area under the ssp585 scenario over the period 2051-2100 in Gabon, Bhutan, Venezuela, Equatorial Guinea, Malaysia, Peru, Guyana, Lebanon, Japan, Congo, Bolivia, Honduras, Zambia, South Korea, Papua Nuova Guinea and other 16 countries. Under the same scenario, the list of 'winner' countries is much shorter, with only Somalia gaining at least 15% resilience over more than 50% of its territory. Countries having the largest increase of vegetation production resilience under the ssp245 scenario in the mid-long term are Russia and the Northern European countries. Under the ssp126 scenario, the vegetation production resilience increases are geographically spread into more continents. In both ssp126 and ssp245 scenarios, there are almost no countries experiencing more than 15% losses of resilience over more than 10% of their lands.

Discussion and conclusions
This study implements a simple and robust statistical indicator, the Annual Production Resilience Indicator (R p ), able to provide a bulk estimation of ecological resilience from annual production time-series. Being focused on production, R p is not sensitive to changes in ecosystem structure and other state variables that play a fundamental part in Holling's definition of ecological resilience (Holling 1973). Conversely, it is formally consistent with more modern conception of resilience allowing adaptation and other changes in ecosystem structure in a resilience framework (Walker et al 2004). Other indicators based on higher temporal resolution data could be employed as well in order to assess more detailed aspects of resilience in a non-linear dynamics framework (Yi and Jackson 2021). However, R p already provides a quick and useful indication of the changes in statistical stability of the production time series in different periods, and it is especially suitable for being computed on large data volumes. Here, the vegetation GPP resilience changes are estimated with reference to the climate variability simulated by a large ensemble of state-of-the-art ESMs.
While the main properties of R p are already demonstrated (Zampieri et al 2019, 2020b, Zampieri 2021, as also summarised in the Appendix A of the Supplementary Material (available online at stacks.iop.org/ERL/16/105001/mmedia), here we developed an interesting decomposition of the resilience changes in terms of the mean GPP and its variability. Such decomposition provides a potentially useful assessment for adaptation planning that should presumably account differently for the changes in the mean with respect to the changes in variability of vegetation GPP in the future climate projections (Mearns et al 1997, Smit et al 2000, van der Wiel and Bintanja 2021. However, the implications of the estimated GPP resilience changes need to be interpreted and confirmed through dedicated local-scale assessments, also identifying the specific drivers of the changes. For instance, while gains in resilience could be generally considered a positive feature in dry environments, their consequences for wet environments and over high latitudes should be further investigated considering more complex approaches (Drews andGreatbatch 2017, Laamrani et al 2020).
Our results aim at providing a useful firstorder indication that surely needs to be corroborated with further dedicated studies also accounting for different socio-economic perspectives. However, the analysis presented here is already potentially very useful to identify the world regions where there might be serious losses of vegetation GPP resilience as well as the countries that are subject to the most urgent necessity of improving adaptive capacity and resilience to climaterelated hazards under different future emission scenarios.
Our results show large differences in the changes of GPP resilience across the globe, depending on greenhouse gases concentration of the projected scenarios. Under low emission scenarios, as found in previous studies (Hubau et al 2020), the CO 2 fertilization effect often prevail over the negative effect due global warming (e.g. an increase in metabolic costs which in turn would lead to a reduction in forest productivity and their efficiency, Collalti et al 2020) and to the increase of climate variability. We find a general increase of vegetation GPP resilience over boreals regions that is in general agreement with the already observed changes of vegetation distribution in those areas (Myneni et al 1997). This tendency is expected to be reinforced in the future climate scenarios, especially those with higher greenhouse gases emissions (Zampieri and Lionello 2010). However, radiation will always be a limiting factor for the vegetation development at very high latitudes (Seddon et al 2016).
The main findings point to areas in the midlatitudes where vegetation resilience is estimated to decrease in the higher emission scenarios, such as the Mediterranean, the mid-West in the US, Central America, part of China, Southern Africa and Australia. The stability of agricultural production and the reliability of ecosystem services provided by the natural vegetation might be compromised in these regions, unless sensible adaptation actions are taken. The importance of climate change mitigation is most evident under the higher emission scenarios, where vegetation resilience is affected in most land areas and especially in tropical regions, where society is highly dependent on ecosystem services and more vulnerable to climatic changes.
The results of our analysis strongly support the SDG-13 on taking action to combat climate change and its impacts. Over areas with a high level of anthropization, our results are relevant for agricultural production, which is a main source of employment, livelihood and income for a large portion of population especially in developing countries (SDG-1, no poverty) as well as a main food source (SDG-2, no hunger). Our results are also relevant for the SDG-15, on the sustainable management of ecosystems and halting land degradation and biodiversity loss.
Overall, in the scenario with lower climate change mitigation (i.e. the Fossil-fuel Development scenario), the areas losing vegetation resilience are larger than the ones gaining resilience, jeopardizing the stability of the ecosystems structure (and of the related services). Adapting to changes in variability more than to changes in the mean production of vegetation will be critical for society and natural ecosystems in areas experiencing large vegetation GPP resilience losses.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// esgf.llnl.gov/.