Annual Green Water Resources and Vegetation Resilience Indicators: Definitions, Mutual Relationships, and Future Climate Projections

Satellites offer a privileged view on terrestrial ecosystems and a unique possibility to evaluate their status, their resilience and the reliability of the services they provide. In this study, we introduce two indicators for estimating the resilience of terrestrial ecosystems from the local to the global levels. We use the Normalized Differential Vegetation Index (NDVI) time series to estimate annual vegetation primary production resilience. We use annual precipitation time series to estimate annual green water resource resilience. Resilience estimation is achieved through the annual production resilience indicator, originally developed in agricultural science, which is formally derived from the original ecological definition of resilience i.e., the largest stress that the system can absorb without losing its function. Interestingly, we find coherent relationships between annual green water resource resilience and vegetation primary production resilience over a wide range of world biomes, suggesting that green water resource resilience contributes to determining vegetation primary production resilience. Finally, we estimate the changes of green water resource resilience due to climate change using results from the sixth phase of the Coupled Model Inter-comparison Project (CMIP6) and discuss the potential consequences of global warming for ecosystem service reliability.


Introduction
The contribution of ecosystems to human well-being has received growing attention in the public and scientific domains. Following the multilateral treaty of the 1992 International Convention on Biological Diversity of the United Nations, the 2010 Conference of Parties (COP 10) formally engaged on a strategic plan to "take effective and urgent action to halt the loss of biodiversity in order to ensure that by 2020 ecosystems are resilient and continue to provide essential services", explicitly recognizing the need for resilient ecosystems, able to withstand the increasing pressure from climate change and anthropization. More recently, in 2015, world leaders adopted the Sustainable Development Goals (SDG), which among others include climate action (SDG13) and life on land (SDG 15) to ensure conservation and restoration of vegetation and sustainable use of the related ecosystems services [1].
The scientific community has provided evidence on the key role of ecosystems for human well-being [2], and has discussed the economic valuation of ecosystem services to mainstream their However, how to measure terrestrial ecosystems' ecological status and their ability to withstand the increasing pressure of climate change and anthropization-i.e., their resilience-is still an open question.
Data from satellite remote sensing have been increasingly adopted for assessing vegetation resilience and ecosystem services, as they can provide information on the land cover typology, extent and characteristics at different points in time, and from the local to the global scale [11,12]. Remotely sensed data can be used to define vegetation indexes such as the Normalized Difference Vegetation Index (NDVI) that is strongly correlated with the photosynthetic activity of the plant canopies [13]. The NDVI is frequently used to assess ecological responses of vegetation to climatic conditions and environmental changes [12,14,15] and has often been adopted as proxy for ecosystem services [16][17][18][19][20].
The NDVI was used as a key indicator in a wide range of ecosystem services including (according to the Millennium Ecosystem Assessment classification): provisioning services, such as energy, food, raw materials, water provisioning, and genetic resources; regulating services, such as climate regulation, disturbance regulation, erosion regulation, flood regulation, hazard regulation, pollination, water purification, and waste treatment; supporting services, such as biological refugia, habitat, primary production in sea and rivers; cultural ecosystem services, such as cultural heritage and recreation [21].
Until recently the focus has been mainly on the assessment of the current flow of ecosystem services and their relative changes compared to past conditions [6,22], which has been possible by the availability of satellite data time series. Given the importance of ecosystem resilience for the reliability of the provided services, several studies have addressed the resilience of ecosystems through indicators of ecosystem stability [23,24]. However, additional research is needed to quantify the resilience of ecosystems to environmental change in spatially explicit assessments [25].
Here, we focus on the estimation of terrestrial ecosystem resilience related to the stability of primary production, which is generally recognized as one of the main constituent of ecosystem resilience [26]. We estimate vegetation primary production resilience by means of the NDVI and its inter-annual variability. More specifically, we adopt the recently proposed indicator for annual crop production resilience [27]. This indicator is derived from the original ecological definition of resilience i.e., the measure of the largest stress that the system can absorb without losing its function [28]. We apply the annual resilience indicator to remotely sensed NDVI, as proxy of vegetation primary production. Furthermore, we apply the annual resilience indicator to ground-based precipitation data, as proxy of annual green water resources (i.e., the water available to plants [29]).
We show and discuss three main results.

1.
Validity of the fundamental principle of the resilience indicator [27] applied to precipitation, i.e., its proportionality with drought return times.

2.
Relationship between the water resource resilience indicator and the vegetation primary production resilience.

3.
Changes of green water resource resilience due to global warming.

Materials and Methods
The original meaning of the term resilience-first introduced in ecology-refers to the largest pressure that a system can cope with without changing its internal structure and losing its functioning capacity [28,30]. While the concept of resilience has evolved in complexity while spreading into other fields such as engineering and social sciences, a recent study [27] brought back the original ecological definition of resilience. The only assumption of that study is that severe shocks are rarer than moderate perturbances. This allows measuring the severity of the shock by its return period and to formally derive a consistent statistical indicator that can be computed on annual production time series [27]. This section summarizes the main properties of the resilience indicator (R C ) and introduces its application to NDVI and precipitation data.

Annual Crop Production Resilience Indicator (R C ): Definition and Properties
The resilience indicator for annual crop production [27] is defined as the squared mean divided by the squared standard deviation of annual production time series: It has been demonstrated [27] that the R C indicator provides several advantages:

1.
It is formally derived from the ecological definition of resilience, thus, theoretically more grounded than similar indicators based on different functions of the µ over σ ratio such as the coefficient of variance.

2.
It is inversely/directly proportional to the frequency/return period of the extreme events leading to large production losses.

3.
It takes into account spatial heterogeneities and diversity in a simple and intuitive manner i.e., R C computed on the sum of n uncorrelated time series with same µ and σ is exactly n-times R C of the individual time series.

4.
It is simple to compute, and it can take into account the effects of non-linear long-term trends easily e.g., by normalizing the time series by the running mean baseline values prior to the indicator computation.
Unfortunately, R C suffers from the disadvantage of requiring relatively long time series for an accurate estimation of resilience. For normally distributed variables, a sample size over 30 values is typically needed to compute R C with and accuracy of 30% [27].

Annual Green Water Resources Resilience Indicator (R P ): Definition and Data
The annual green water resources resilience indicator (R P ) follows the same formulation as the R C indicator (Equation (1)) but it is computed using annual precipitation data i.e., the main input variable that defines the so-called green water [29]. In this work, we use global precipitation data since 1901 at half-degree spatial resolution from the Climate Research Unit (CRU) archive [31].
We test the rationale behind the R C indicator applied on precipitation (i.e., point 2 of Section 2.1) by comparing R P values to the frequency of meteorological drought estimated independently with the Weighted Anomaly of Standardized Precipitation (WASP) index [32,33] (i.e., the same data used in Figure 3 of Carrão et al.'s paper [33]).
To provide a preliminary estimate of future change of water resource resilience due to climate change, we compute R P using data from the sixth phase of the Coupled Model Inter-comparison Project (CMIP6). At the time this analysis has been developed, climate change simulations performed by six climate models (CanESM5, CNRM-CM6, IPSL-CM6A-LR, MIROC6, MRI-ESM2-0 and UKESM1-0-LL) were available for the historical period (1850-2014) and Shared Socio-economic Pathways SSP585 [34] for period 2015-2100. The SSP585 corresponds to a high emissions scenario consistent with the Representative Concentration Pathway RCP 8.5 used in CMIP5.

Annual Vegetation Primary Production Resilience Indicator (R V ): Definition and Data
The annual vegetation primary production resilience indicator (R V ) is analog to the R C indicator (Equation (1)) but it is computed using annual NDVI data, which we use as proxy of gross primary production. NDVI time series spanning the period July 1981 to December 2015 were obtained from the Global Inventory Modeling and Mapping Studies(GIMMS) NDVI3g v1 dataset [35]. The long-term GIMMS dataset is derived from measurements made by different Advanced Very High Resolution Radiometer (AVHRR) sensors on board of the National Oceanic and Atmospheric Administration (NOAA) polar-orbiting satellite series [35]. The satellite observations are corrected for various effects, including sensor degradation, orbital drift, atmospheric effects, and volcanic eruptions. The time series have a temporal resolution of about 15 days (i.e., two observations per month, one referring to day 1-15 and one to day 16-end of the month) and a spatial resolution of 1/12 degree. The 15-day composites were obtained using the maximum value composition (MVC) technique [36], which minimizes the scattering and absorption effects of the atmosphere.
We temporally smoothed the time series to reduce remaining noise. We applied an iterative Savitzky-Golay filter [37] to fit the upper envelope of the NDVI time profiles as described by Chen et al. [38]. Finally, the annual mean NDVI values were calculated for the period 1982-2015.
We aggregate the NDVI data at 0.5-degree resolution to allow analyzing the relationship between R V and R P . Water bodies are not considered in the aggregated values, which are achieved through a conservative mapping procedure [39]. We discuss the relationships between R V and R P for each world biomes grouping the results by homogeneous climatic areas according to the Köppen-Geiger (KG) classification [40]. The KG classification codes are listed in Table 1.

Properties of R P and R V
The R P and R V indicators inherit all properties of R C (Section 2.1). In addition, the R P and R V indicators computed at the grid cell level have an important property, i.e., they neither depend on the period of the year nor on the portion of area with no precipitation or vegetation within the analyzed cell (P = 0, NDVI close to 0).
Assuming linear mixing, this claim can be proven by computing the mean and the variance of a time series of values x i , i = 1,.., N and of derived values defined as y i = a × x i , where a is the fraction of area or the fraction of period where precipitation or NDVI are not null. It follows that µ y = a × µ x and σ y = a × σ x , so the resilience indicator computed on these two time series does not change (q.e.d.). Please note that if a trend exists in the period of the year or portion of area with no precipitation or vegetation, this should be removed as mentioned in point 4 of Section 2.1.
Because of the 3rd property of the resilience indicator (Section 2.1), when aggregating the resilience indicator computed at the grid cell level over an area of interest, the resulting resilience will be larger if groups of temporally uncorrelated (i.e., diverse) cells exist in the area.

Results
The three main findings of this study are here shown and discussed: Remote Sens. 2019, 11, 2708 5 of 13 1.
The test of the annual resilience indicator applied to precipitation to be inversely proportional to drought frequency (Section 3.1-2nd property).

2.
The relationship between R P and R V (Section 3.2). 3.
The effects of climate change on R P (Section 3.3).
3.1. Annual Green Water Resources Resilience Indicator (R P ) and Drought Frequency Since 1901 Figure 1 shows the mean (Figure 1a), the standard deviation (Figure 1b) of observed annual precipitation data since 1901, as well as their combination defining R P (Figure 1c). R P tends to be larger in regions with high precipitation values. However, the spatial variability is larger than the one of the precipitation magnitude due to complex patterns of precipitation inter-annual variability. Inter-annual variability is probably underestimated in regions with scarce data such as central Greenland, Western Arabian Peninsula, and Central Amazon Forest, hindering the accurate computation of the annual green water resource resilience indicator.    (Figure 1c). RP tends to be larger in regions with high precipitation values. However, the spatial variability is larger than the one of the precipitation magnitude due to complex patterns of precipitation inter-annual variability. Interannual variability is probably underestimated in regions with scarce data such as central Greenland, Western Arabian Peninsula, and Central Amazon Forest, hindering the accurate computation of the annual green water resource resilience indicator. Figure 1d shows the drought hazard computed as the probability of exceeding the median of global severe precipitation deficits for the reference period using the WASP [32]). The data displayed in Figure 1d is exactly the same as in Figure 3 of [33] and shows the frequency of abnormal precipitation deficits defining drought prone regions. Comparison between Figure 1c and Figure 1d strongly suggests an inverse relationship between RP and drought frequency. Figure 2 shows the relationship between RP and drought return periods. Their joint distribution is consistent, despite the large scatter of the data, with a global tendency toward the increase of drought return periods with increasing values of the precipitation resilience indicator. This consistency was tested using the proportionality coefficient that we found to be statistically  [28,29]). Figure 1d shows the drought hazard computed as the probability of exceeding the median of global severe precipitation deficits for the reference period using the WASP [32]). The data displayed in Figure 1d is exactly the same as in Figure 3 of [33] and shows the frequency of abnormal precipitation deficits defining drought prone regions. Comparison between Figure 1c,d strongly suggests an inverse relationship between R P and drought frequency. Figure 2 shows the relationship between R P and drought return periods. Their joint distribution is consistent, despite the large scatter of the data, with a global tendency toward the increase of drought return periods with increasing values of the precipitation resilience indicator. This consistency was tested using the proportionality coefficient that we found to be statistically significant with a p-value < 0.001.

Precipitation and Vegetation Resilience Indicator
Since 1982. Figure 3 shows the mean (Figure 3a), the standard deviation (Figure 3b) of annual mean NDVI data since 1982, as well as their combination defining RV (Figure 3c). RV is larger in tropical regions, where the NDVI is large and the inter-annual variability is relatively low. RV sharply decreases in subtropical regions, characterized by large inter-annual NDVI variability. RV also decreases towards the high latitudes, where vegetation is generally less vigorous and subject to the largest inter-annual variability.   Figure 3 shows the mean (Figure 3a), the standard deviation (Figure 3b) of annual mean NDVI data since 1982, as well as their combination defining R V (Figure 3c). R V is larger in tropical regions, where the NDVI is large and the inter-annual variability is relatively low. R V sharply decreases in subtropical regions, characterized by large inter-annual NDVI variability. R V also decreases towards the high latitudes, where vegetation is generally less vigorous and subject to the largest inter-annual variability.     Table 2 shows the overall statistics organized per Köppen-Geiger (KG) classification. In general, R V tends to be larger than R P . We note a relatively large spread of the precipitation data, i.e., large standard deviation with respect to the mean precipitation averaged in each KG climate zone. The spread of the drought frequency estimate by the WASP and the vegetation data measured by the NDVI is smaller. The latter result is somehow expected because the KG climate classification was developed exactly to characterize homogeneous land ecosystems. The spread of the green water resources and vegetation primary production resilience indicators over homogeneous climate zones is also relatively large, probably reflecting the dependency of the resilience indicator accuracy computed on time series of limited length [23,27]. However, the mean R V is often significantly different among KG classes. R P is higher in equatorial (A class) and warm temperate climates (C class) than in arid climate (B class). Similarly, R V decreases from equatorial (A class) to warm temperate (C class), arid (B class) and snow (D class) climates. Within warm temperate biomes (C class), R P is greater in fully humid (f class) than in winter and summer dry regions (w and s classes, respectively), while R V takes higher values in winter dry and summer dry (with warm summer) climates (wb and sb classes, respectively) rather than in fully humid climates (f class). In arid climates (B class), R P is greater in the steppe than in the desert regions (S and W classes respectively), while lower values of R V are found in cold areas than in hot areas (class k and h respectively). In snow climates (class D), regions with cool and warm summer (class c and b) have higher R V than regions with hot summer or extremely continental weather (a and d classes, respectively). The size of the markers represents the land area actually covered by the homogeneous climate zones. The color scheme is the same as in Figure 4d. Figure 4 shows the relationships between R P and R V averaged over the climate zones defined according to the Köppen-Geiger classification ( The size of the markers represents the land area actually covered by the homogeneous climate zones. The color scheme is the same as in Figure 4d. Figure 4 shows the relationships between RP and RV averaged over the climate zones defined according to the Köppen-Geiger classification ( Table 2).

Figure 4.
Ratios of vegetation primary production and green water resource resilience indicators (RP and RV, respectively) vs RP averaged over the homogeneous KG climate zones displayed in Figure  3d.
Despite the large scatter of the data that emerged from the analysis in Table 2, the comparison between RP and RV for different climate zones displays clear relationships. In fact, the ratio between the two is roughly constant over two wide ranges of biomes. On the one hand, the RV / RP ratio tends to stay around the 10-15 band for ecosystems ini arid climate (i.e., the B class), warm temperate climate (C class) and equatorial climate (A class). On the other hand, the RV / RP ratio tends to stay below 5 for snow and polar climates (i.e., the "D" and E classes, respectively) as well as the subpolar oceanic climate (Cfc), representing an exception to the previous statement. Another exception is represented by the warm-summer Mediterranean climate (Csb), characterized by the largest RV / RP ratio. This peculiar climate is found in the north-western Iberian Peninsula and central Anatolia (see Figure 3d).

Effects of Climate Change on Annual Green Water Resources Resilience Indicator (RP).
Figure 5 shows the changes in mean annual precipitation (Figure 5a), standard deviation ( Figure  5b) and of RP (Figure 5c) computed from an ensemble of six climate models from CMIP6 between the periods 1981-2014 and 2015-2100. The mean annual precipitation and the standard deviation are increasing in many regions around the globe and mainly in the Northern Hemisphere, because of the larger humidity levels that the warmer atmosphere can hold, because of alterations in the frequencies of large-circulation patterns and because of the general increase of climate variability [41,42]. According to these simulations, the relative ratio of the increases in the mean and in the standard deviation of annual precipitation is such that the annual RP is decreasing over most regions (Figure 5c and d). Despite the generalized increase in precipitation, the significant increase in inter- Figure 4. Ratios of vegetation primary production and green water resource resilience indicators (R P and R V , respectively) vs. R P averaged over the homogeneous KG climate zones displayed in Figure 3d.
Despite the large scatter of the data that emerged from the analysis in Table 2, the comparison between R P and R V for different climate zones displays clear relationships. In fact, the ratio between the two is roughly constant over two wide ranges of biomes. On the one hand, the R V /R P ratio tends to stay around the 10-15 band for ecosystems ini arid climate (i.e., the B class), warm temperate climate (C class) and equatorial climate (A class). On the other hand, the R V /R P ratio tends to stay below 5 for snow and polar climates (i.e., the "D" and E classes, respectively) as well as the subpolar oceanic climate (Cfc), representing an exception to the previous statement. Another exception is represented by the warm-summer Mediterranean climate (Csb), characterized by the largest R V /R P ratio. This peculiar climate is found in the north-western Iberian Peninsula and central Anatolia (see Figure 3d). Figure 5 shows the changes in mean annual precipitation (Figure 5a), standard deviation ( Figure 5b) and of R P (Figure 5c) computed from an ensemble of six climate models from CMIP6 between the periods 1981-2014 and 2015-2100. The mean annual precipitation and the standard deviation are increasing in many regions around the globe and mainly in the Northern Hemisphere, because of the larger humidity levels that the warmer atmosphere can hold, because of alterations in the frequencies of large-circulation patterns and because of the general increase of climate variability [41,42]. According to these simulations, the relative ratio of the increases in the mean and in the standard deviation of annual precipitation is such that the annual R P is decreasing over most regions (Figure 5c,d). Despite the generalized increase in precipitation, the significant increase in inter-annual variability is reflected in a reduction of R P over most regions around the globe. This means that even though there will be more water available, this will be less reliable. The only exception is found in desert regions, where small increases in mean precipitation are dominating the R P signal.

Effects of Climate Change on Annual Green Water Resources Resilience Indicator (R P )
annual variability is reflected in a reduction of RP over most regions around the globe. This means that even though there will be more water available, this will be less reliable. The only exception is found in desert regions, where small increases in mean precipitation are dominating the RP signal.

Discussion
The functional quality of ecosystems is directly linked to the amount and reliability of services they can provide [43,44]. With this paper, we enrich the set of analysis tools for assessing the quality and the stability of ecosystems through the resilience indicators that we propose to apply to green water resources and to vegetation primary production.
These indicators, formally derived from the ecological definition of resilience [27], require relatively long time series (30 years and more) that, despite not available from the most frequently used satellite mission (e.g., MODIS, SPOT-VGT and Prob-V), can be extracted from AVHRR-based long-term archive as the one used in this study (i.e., GIMMS). This long-term approach may thus complement more complex resilience measures based on the engineering definition of resilience i.e., the ecosystem recovery time more suited for short time-scales analyses [45,46]. Comparing these two different views would probably provide an interesting way forward to the understanding of ecosystems resilience and stability assessments.
Our analysis linking annual precipitation and vegetation productivity suggests a coherent relationship between the resilience indicators computed on precipitation and NDVI. The magnitude of the resilience indicator for vegetation was found larger than the one of precipitation. This may indicate the capability of vegetated ecosystems to maintain their functioning also in the presence of

Discussion
The functional quality of ecosystems is directly linked to the amount and reliability of services they can provide [43,44]. With this paper, we enrich the set of analysis tools for assessing the quality and the stability of ecosystems through the resilience indicators that we propose to apply to green water resources and to vegetation primary production.
These indicators, formally derived from the ecological definition of resilience [27], require relatively long time series (30 years and more) that, despite not available from the most frequently used satellite mission (e.g., MODIS, SPOT-VGT and Prob-V), can be extracted from AVHRR-based long-term archive as the one used in this study (i.e., GIMMS). This long-term approach may thus complement more complex resilience measures based on the engineering definition of resilience i.e., the ecosystem recovery time more suited for short time-scales analyses [45,46]. Comparing these two different views would probably provide an interesting way forward to the understanding of ecosystems resilience and stability assessments.
Our analysis linking annual precipitation and vegetation productivity suggests a coherent relationship between the resilience indicators computed on precipitation and NDVI. The magnitude of the resilience indicator for vegetation was found larger than the one of precipitation. This may indicate the capability of vegetated ecosystems to maintain their functioning also in the presence of precipitation variability. Nevertheless, the spatial patterns of the two resilience indicators do coincide to a great extent. In regions with reduced agreement, as for instance in the high latitudes of the Northern Hemisphere, the reduced vegetation resilience as compared to the precipitation resilience may be due to temperature and/or radiation, not explicitly considered in this study. Further studies are needed to assess the robustness of these results including, for instance, the effects of temperature and radiation as well.
This analysis could also be further evolved by accounting for the delay between precipitation input and vegetation response, as well as different time-scales [45,47,48], including different seasonalities such as double seasons and different periodicities as those found in the Southern Hemisphere and disentangling the differences between natural and anthropic vegetation cover, possibly isolating the effects of management and irrigation [49].

Conclusions
This paper provides an alternative framework for the assessment of ecosystem resilience by proposing two indicators, one for green water resources and one for vegetation productivity, based on the ecological definition of resilience [27].
Our global scale analysis strongly suggests a coherent relationship between these indicators over two wide ranges of biomes in the present climate. However, further studies are needed to quantify the robustness of these relationships in a changing climate, considering other definitions of resilience as well [23,45]. Including temperature and radiation would be also needed to fully explain the observed spatial patterns in vegetation resilience.
Quantifying the resilience of ecosystems and ecosystem services is of upmost relevance when considering future climate change [50][51][52][53]. Indeed, the variation of temperature and precipitation patterns caused by climate change is expected to affect the availability and seasonality of water resources and the growing season of vegetation in various regions of the world [54], affecting biodiversity and the provisioning of ecosystem services [24,46,55].
Our preliminary analysis of climate change projections from CMIP6 data displays alarming results. The increase of precipitation magnitude projected for the next future is counterbalanced by the concomitant increase of its variability, resulting in a general decrease of green water resource resilience at the global level.
In the future, it would also be interesting to carry on this analysis with an indicator of Blue Water Resources Resilience, which could be easily defined by computing the resilience indicator on the surface water budget given by precipitation minus evapotranspiration. This could be put in relation with soil moisture and river drought through the standardized Precipitation Evapotranspiration Index (SPEI [56,57]) or the Standardized River Discharge Index (SRDI [58,59]) to measure the reliability of freshwater resources for society and to help assessing the status of river ecosystems [60]. This additional resilience indicator may be potentially suitable for river related ecosystem services.