Regional Greening as a ‘Positive’ Tipping Phenomenon

Yu Sun School of Systems Science/Institute of Nonequilibrium Systems, Beijing Normal University, Beijing 100875, China Teng Liu School of Systems Science/Institute of Nonequilibrium Systems, Beijing Normal University, Beijing 100875, China Shang Wang School of Systems Science/Institute of Nonequilibrium Systems, Beijing Normal University, Beijing 100875, China Jun Meng School of Science, Beijing University of Posts and Telecommunications, Beijing 100876, China Potsdam Institute for Climate Impact Research, Potsdam 14412, Germany Yongwen Zhang Data Science Research Center, Faculty of Science, Kunming University of Science and Technology, Kunming 650500, China Saini Yang State Key Laboratory of Earth Surface Processes and Resource Ecology, Beijing Normal University, Beijing 100875, China School of National Safety and Emergency Management, Beijing Normal University, Beijing 100875, China Xiaosong Chen School of Systems Science/Institute of Nonequilibrium Systems, Beijing Normal University, Beijing 100875, China Deliang Chen Department of Earth Sciences, University of Gothenburg, Gothenburg 40530, Sweden Jürgen Kurths Potsdam Institute for Climate Impact Research, Potsdam 14412, Germany Department of Physics, Humboldt University, Berlin 10099, Germany Shlomo Havlin Department of Physics, Bar Ilan University, Ramat Gan 52900, Israel Hans Joachim Schellnhuber Potsdam Institute for Climate Impact Research, Potsdam 14412, Germany Jingfang Fan [email protected] School of Systems Science/Institute of Nonequilibrium Systems, Beijing Normal University, Beijing 100875, China Potsdam Institute for Climate Impact Research, Potsdam 14412, Germany

(December 7, 2024)

Abstract

Earth system tipping elements have been predominantly investigated for their potential to trigger negative ecological, climatic, and societal shifts. Yet, an overlooked but seminal avenue exists in the form of positive tipping phenomena, whose underlying mechanisms and benefits remain largely underexplored. To bridge this gap, our research introduces a fundamental percolation-based framework to assess the criticality and resilience of planetary terrestrial vegetation systems. Leveraging high-resolution satellite data, we focus on greening-induced positive tipping dynamics driven by global warming. We feature the Qinghai-Tibetan Plateau (QTP) and the Sahel region as contrasting yet analogous case studies. Our analysis uncovers an intriguing phenomenon where vegetation fragmentation aligns with a percolation threshold, exhibiting a scale-invariant pattern characterized by nearly perfect power laws with three critical exponents. Remarkably, contrary to conventional destructive tipping elements, these regions act as favorable tipping elements, transitioning from fragmented to cohesive vegetation patterns due to anthropogenic climate change and afforestation efforts. Furthermore, we propose an optimal resilience enhancement model to reinforce vegetation robustness while minimizing socio-economic costs. This study provides valuable insights into the favorable aspects of tipping elements under climate change and offers effective strategies for enhancing ecological resilience against environmental threats.

I Main

Tipping elements, integral components of the Earth System, possess the capacity for sudden and irreversible state shifts when specific thresholds are exceeded [1, 2, 3]. These tipping points can precipitate catastrophic cascading failures across ecological, climatic, and societal domains, often in response to minor perturbations. Such tipping points are not restricted to Earth System but manifest in a variety of other complex systems, including financial markets [4, 5], ecosystems [6, 7, 8, 9, 10, 11], and social dynamics [12, 13]. While the majority of research has focused on the negative aspects of tipping points—namely, their potential to induce collapses or disrupt ecosystems—certain tipping elements can also drive favorable outcomes. These positive tipping elements offer a largely untapped avenue for mitigating the adverse impacts of climate change by fostering resilience and adaptive capacity across diverse domains. In this paper, we introduce a novel percolation-based framework to assess the criticality and resilience of planetary vegetation systems. We focus on two contrasting yet analogous regions—the Qinghai-Tibetan Plateau (QTP) and the Sahel—which serve as case studies for examining the complex interplay of tipping elements within the broader context of climate change and human activities.

The QTP, often referred to as the “Third Pole”, spans an impressive 2.5 million km² with an average elevation of 4.5 km. Its massive glacier storage, second only to the Antarctic and Arctic, feeds Asia’s ten largest rivers, serving nearly 800 million people downstream [14, 15, 16, 17]. The Plateau’s unique topography influences various atmospheric, cryospheric, hydrological, and biospheric processes [16, 17, 14, 18, 19, 20], leading to broad climate patterns and establishes large-scale teleconnections with various climatic phenomena and tipping elements, such as the Indian Summer Monsoon [21, 22, 23], East Asian Summer Monsoon [24, 23], westerly jet [25], Amazon rainforest [26], Arctic Oscillation [27], El Niño-Southern Oscillation [28, 29], dust storms [30], and others (Fig. 1a). Alarmingly, it has experienced a rate of warming twice the global average over the past five decades [31], revealing its significant sensitivity to climate change. Terrestrial vegetation on the QTP, primarily alpine steppe and meadow, plays a vital role in both local and global systemic stability. Climate change has a profound impact on QTP vegetation, affecting processes such as photosynthesis [32, 33, 34] and tree-line patterns [35, 36]. This leads to feedback mechanisms including evaporative cooling [19], carbon sequestration [37], and complex, self-reinforcing interconnections involving vegetation, air, ice, water, wildlife, and human activity (Fig. 1b).

In contrast, the Sahel region acts as Africa’s climatic and ecological pivot, lying between the arid Sahara and the humid forests of West and Central Africa [38]. Its northern fringe is defined by the Sahara desert, while its southern edge reveals a sharp gradient in vegetation, demarcating the shift towards the monsoonal forests characteristic of West and Central Africa. Its unique position makes it particularly vulnerable to the impacts of climate change and anthropogenic activities, a fact extensively supported by scholarly research [39, 40, 41]. Recent fluctuations in precipitation patterns have exerted discernible effects on both ecological systems and human societies [42, 43]. Notably, the Sahara has undergone periodic wet phases throughout the Quaternary period, with a proliferation of vegetation [44]. Future projections suggest that climate change could further amplify this trend, resulting in even greater rainfall in the region [45]. Despite these insights, pinpointing the specific causative factors underlying these changes remains a complex analytical challenge [46, 47, 48, 49].

From a theoretical standpoint, understanding tipping points is crucial for predictive and adaptive strategies. Early-warning signals like lag-1 autoregression, variance, skewness, and detrended fluctuation analysis often rely on the principle of critical slowing down (CSD), which refers to the system’s progressively slower recovery from perturbations as it approaches a tipping point [50, 51, 52]. Yet, the capacity of CSD-based analyses to capture the geometric complexity and emerging functional structures that stem from spatiotemporal variability within the system might be limited. Additionally, the accuracy and effectiveness of such indicators may be hindered by inadequate empirical observations over time.

To overcome the limitations of conventional early-warning signals, our percolation-based framework employs high-resolution satellite data to evaluate the self-organized criticality (SOC) [53] and resilience of the vegetation in both the QTP and the Sahel region. Drawing a parallel to occupancy probability in percolation theory, each grid point represents different vegetation types, characterized by an associated Enhanced Vegetation Index (EVI) value indicating the degree of greenness. Higher EVI values correspond to more abundant and thriving vegetation, similar to higher occupancy probability signifies a greater likelihood of a site being occupied in a percolation model. In our framework, when neighboring cells have EVI values surpassing a predefined threshold (EVI_c), they form a vegetation fragment (or cluster) [54]. This process is illustrated in Fig. 1c-e. We present photos taken during June-July of 2021 and 2022, displayed in Fig. 1c, showcasing the diversity of vegetation present in the QTP. The corresponding EVI values are depicted as bar heights in Fig. 1d and values in Fig. S1a. Additionally, after removing nodes with EVI values below EVI_c, the configuration of vegetation fragments can be observed in Fig. 1e and Fig. S1b.

This framework delivers a unique spatial perspective, enabling the identification of functional fragmentation structures – the disintegration of ecosystems into isolated units – and the detection of geometric tipping properties like scale invariance and self-similarity in vegetation patterns. Notably, our calculated percolation threshold EVI_c approaches an empirically established ecological tipping point of approximately 0.2, delineating the accepted threshold between sparse or nonexistent vegetation and regions of healthy vegetation over the past two decades. Concurrently, the accompanying power law distributions also support the approach to this tipping point, exhibiting three identical exponents for fragment size distribution ( $\tau$ ), mass fractal dimension ( $D_{f}$ ), and external diameter fractal dimension ( $D_{e}$ ). Through this method, critical areas requiring protection can be identified, facilitating the development of effective strategies to tackle the challenges posed by factors such as climate change. This insight allows policymakers and conservationists to prioritize their efforts and allocate resources more efficiently by pinpointing regions that are particularly vulnerable or essential for ecosystem health.

II Results

To conduct our percolation analysis [55], EVI [56], a well-established index of greenness and more effective than other vegetation indices, such as the Normalized Difference Vegetation Index, is used. Our analysis utilizes high spatial resolution QTP EVI maps (25.99°N-39.82°N, 73.49°E-104.63°E), which contain approximately 85.5 million pixels, and Sahel EVI maps (9.45°N-20.07°N, 17.40°W-38.29°E), which contain approximately 117.3 million pixels. Each pixel measures 250m $\times$ 250m. Our study focuses on the QTP and Sahel region during the summer months (July-September) of the Northern Hemisphere. This period coincides with the peak growing season when vegetation cover is at its highest within these area. Please refer to the Data Section (IV.1) for comprehensive details and references regarding the data used in this study.

II.1 Criticality of vegetation in the Qinghai-Tibetan Plateau and Sahel

To assess the criticality of vegetation in the QTP, we perform a site percolation analysis using the average EVI image from the recent five-year period (2017-2021). The spatial pattern of the EVI is depicted in Fig. 2a. In our analysis, we’re systematically attacking and removing nodes (i.e., pixels) with EVI values that fall below a given threshold. The resulting fragments are tracked at each step using the Newman-Ziff algorithm [57]. Fig. 2b illustrates the relative sizes of the largest fragment, $G_{1}$ , and the second largest fragment, $G_{2}$ , as a function of EVI. We observe a sudden and dramatic jump in the size of the largest fragment ( $G_{1}$ ) at a percolation threshold (EVI_c), which is approximately 0.2 (indicated by the vertical black dashed line). Concurrently, the maximum size of the second-largest fragment ( $G_{2}$ ) also peaks at this threshold. Fig. 2c displays the greenness fragmentation structures at EVI_c, highlighting only fragments with a size (node number) larger than 100. The structure reveals two major fragments, with a critical node at an EVI_c value located near (36.07 ${}^{\circ}N$ , 101.96 ${}^{\circ}E$ ) acting as a connection between these fragments. In percolation theory, the percolation threshold is determined by $G_{1}$ [58] and $G_{2}$ [59]. Ecologically, an EVI value of 0.2 typically represents the threshold between unhealthy (without or sparse vegetation cover) and healthy vegetation [60]. The strong agreement between our theoretical percolation threshold EVI_c and empirical values highlights the criticality of vegetation in the QTP.

Next, we analyse the fragment size distribution, $n_{s}$ , at EVI_c, as described in Section (IV.4). Our analysis reveals a total of over 87,000 fragments distributed throughout the QTP, with sizes spanning seven orders of magnitude. Intriguingly, the fragment size distribution follows a perfect power-law relationship: $n_{s}\sim s^{-\tau}$ , where the exponent $\tau$ is approximately 2.04, as depicted in Fig. 2d. This empirical exponent aligns closely with the theoretical value of $\tau=187/91\approx 2.05$ as predicted by the classical percolation theory for a 2D lattice [55, 61].

Furthermore, we investigate the fractal properties of the giant fragment at EVI_c. Using the box renormalization method (refer to Methods) and examining the mass and external perimeter of the giant fragment at various image resolutions, we confirm that both the inner structure and external boundary of the giant fragment at EVI_c exhibit self-similarity (see Section IV.6). Our analysis, shown in Fig. 2e, yields empirical fractal dimensions for the mass ( $D_{f}\approx 1.93$ ) and external perimeter ( $D_{e}\approx 1.40$ ). Remarkably, these empirical exponents, $D_{f}$ and $D_{e}$ are consistent with the theoretical values of $D_{f}=91/48\approx 1.90$ and $D_{e}=4/3\approx 1.33$ at the percolation threshold for a 2D lattice (refer to Table. 1).

To provide a comparison, we also introduce a null model generated by shuffling the original EVI data (100 independent realizations), which represents the classical uncorrelated site percolation. The results, depicted in Fig. 2d, Fig. S2a, and Fig. S2b, demonstrate consistent power-law characteristics in the fragment size distribution and fractal substructures between the empirical and null models, with similar exponents $\tau$ , $D_{f}$ , and $D_{e}$ (Table. 1) for the QTP. These findings strongly suggest that the greenness structures of the QTP at EVI_c are in a critical state, indicative of a system operating at the SOC stage.

Table 1: Scaling fragmentation of vegetation cover on the Qinghai-Tibetan Plateau (QTP) and Sahel. Critical exponents

\tau

D_{f}

, and

D_{e}

are calculated in the empirical and null models, as well as theoretical exponents corresponding to 2D classical percolation theory at the percolation threshold. At the critical point, the exponent of the fragment size distribution

\tau

is further related to the fractal dimension

D_{f}

via the hyperscaling relationship

\tau=1+d/{D_{f}}

[55, 61, 62].

	QTP	Sahel	Null model	Theory
Fragments size distribution exponent $\tau$	2.04 ( $\pm 0.02$ )	2.06 ( $\pm 0.02$ )	2.06 ( $\pm 0.0006$ )	187/91 [61]
Mass fractal dimension $D_{f}$	1.93 ( $\pm 0.07$ )	1.90 ( $\pm 0.06$ )	1.91 ( $\pm 0.02$ )	91/48 [63]
External diameter fractal dimension $D_{e}$	1.40 ( $\pm 0.04$ )	1.19 ( $\pm 0.03$ )	1.33 ( $\pm 0.01$ )	4/3 [64]

Building on our analytical percolation framework, we extend our assessment to evaluate the Sahel region’s ecological criticality, conducting site percolation analyses for the period 2013-2017. The spatial distribution of the EVI values across the Sahel is represented in Fig. 2f. In a manner analogous to our work on the QTP, Fig. 2g presents the relative sizes of the largest fragment ( $G_{1}$ ) and the second-largest fragment ( $G_{2}$ ), respectively. These are charted against varying EVI values, allowing us to identify the percolation threshold (EVI_c) at approximately 0.2, as indicated by the vertical black dashed line. Fig. 2h illustrates the Sahel’s greenness fragmentation structures at this critical EVI_c level. Specifically, it highlights two major vegetation fragments centered around the critical node near (14.92^∘N, 3.86^∘W). The distribution of fragment sizes at EVI_c also follows a power-law formula $n_{s}\sim s^{-\tau}$ , where the exponent $\tau$ is calculated to be approximately 2.06, as demonstrated in Fig. 2i. Empirical fractal dimensions for mass ( $D_{f}$ $\approx$ 1.90) and external perimeter ( $D_{e}$ $\approx$ 1.19) are determined in Fig. 2j. Notably, the three empirical exponents $\tau$ , $D_{f}$ , and $D_{e}$ for the Sahel are almost consistent with the theoretical values established for a percolation in a 2D lattice (refer to Table. 1).

II.2 Greening-induced tipping point

To investigate the spatio-temporal patterns of vegetation on the QTP and Sahel, we employ percolation procedures on a series (2001-2021) of temporally evolving EVI images, each representing the temporal average of pixel-by-pixel EVI values within a 5-year sliding window. The sliding window moves forward by 2 years, resulting in adjacent windows with a temporal overlap of 3 years. We analyse the percolation phase diagram and determine the EVI_c for each sliding window, as illustrated in Fig. S3. In analogy to percolation theory, the resilience state of QTP’s vegetation can be evaluated through EVI_c. This threshold indicates the critical point at which the cohesive connectivity within the vegetation starts to break down [65]. A higher EVI_c denotes a better state that more extensive vegetation damage is required to disrupt the ecosystem. As depicted in Fig. 3b, we observe an overall trend of EVI_c increasing from 0.17 to 0.23 over the past two decades, indicating a shift toward a positive tipping point and enhanced ecosystem resilience. To further assess the changes in vegetation patterns, we compare the coverage area of the largest fragment in the QTP between the periods of 2001-2005 and 2017-2021 at the same threshold EVI_c=0.23 (before the jump). As shown in Fig. S4, we observe an increase of approximately 16.2% in the coverage area percentage of the QTP’s largest fragment, accompanied by greenness in the northeast and the clustering of numerous small greenness fragments. The enhancement in greenness resilience suggests that the QTP has been transitioning from a fragmented to a cohesive phase, acting as a positive tipping element.

Additionally, we verify the power-law form of the fragment size distribution for each sliding window, as shown in Fig. S5. Remarkably, we observe a significant trend in the empirical fragment size distribution exponent $\tau$ over the period of 2001-2021 for the QTP. As depicted in Fig. 3d, $\tau$ progressively approaches the theoretical value of $\tau\approx 2.05$ . This trend suggests increasing alignment of the fragmentation of QTP’s vegetation cover with classical percolation theory. Furthermore, we analyse the goodness of power-law fits for the fractal dimensions $D_{f}$ and $D_{e}$ (shown in Fig. S6 and Fig. S7). The results reveal that both fractal dimensions also exhibit a significant trend toward the theoretical values, $D_{f}=1.90$ and $D_{e}=1.33$ , as shown in Fig. S8. These findings suggest that since approximately 2010, the greenness of the QTP has been converging towards the tipping point, supported by the proximity of the critical exponents to its theoretical value.

Both observational and modelling studies have demonstrated that recent climate change has impacted the structure and ecological functioning of QTP vegetation [19, 66, 67, 68]. In marked contrast we find here, in Fig. 3a, the spatial pattern of significant EVI change trends from 2001 to 2021 reveals a general positively greening trend, particularly in the Northeast, Northwest, and Southwest regions of the QTP. Fig. S9 reveals distinct temporal fluctuations in the EVI values among the critical nodes identified within the QTP.

To investigate the potential climatic drivers behind the observed criticality and enhanced resilience of the QTP, we examine the relationship between the evolution of EVI and temperature as well as precipitation patterns across the region. Our investigation reveals a similar trend in both, the EVI_c and the warming of the QTP (Pearson correlation coefficient $r$ = 0.88, statistically significant at $p<10^{-4}$ ), as demonstrated by comparing EVI with the average temperature from June to August (Fig. 3c). Additionally, increased precipitation was identified as a key factor contributing to the enhancement of greenness. The variation in vegetation intensity aligns with the variation in precipitation, exhibiting a delay of approximately two years (Fig. S10). Pixel-by-pixel cross-correlation analysis with time lags reveals the relationship between the EVI series and precipitation series, identifying precipitation-induced spatial patterns of vegetation greening with the corresponding time lag, with a maximum lag of up to four years, as illustrated in Fig. S11a during the 2001-2021 period (Fig. 3b). Approximately 34.3% of the QTP’s total area exhibits significant vegetation growth. Interestingly, about 62% of these areas experiencing significant vegetation growth are preceded by a precipitation increase at least one year prior. In particular, in the arid northern regions of the QTP, vegetation exhibits substantial sensitivity to fluctuations in precipitation. In contrast, the southeastern areas, distinguished by plentiful rainfall and verdant vegetation, witness a decrease in vegetation sensitivity in response to an upswing in precipitation (refer to Fig. S11b) [69]. This differential response can be ascribed to the marginal effect of precipitation. In arid regions, even a minor uptick in rainfall can spur significant vegetation growth. Conversely, in more humid locales, intense storms could inflict damage on the terrain, thereby diminishing vegetation coverage. Our analysis highlights that the temporal changes in greenness exhibit a significant positive response to precipitation, with a delay of over one year in areas where greenness has experienced significant enhancement ( $p<$ 0.05). In addition to climate change, human afforestation efforts and projects aimed at ecological conservation and restoration, such as the Grain-for-Green Program [70, 71, 72], have contributed to land improvement and increased vegetation greenness. Millions of hectares of farmland and degraded land have been converted back to forest and grassland, resulting in an increase of forest ecosystems by $0.6\times 10^{4}$ km² [73]. Consequently, the mechanism of the greening-induced tipping point in the QTP is likely driven by a combination of climate change and human factors.

Similar to our analysis of the QTP, we examine the Sahel region’s vegetation dynamics (Fig. 3g shows the spatial pattern of significant EVI change and Fig. 3h indicates the precipitation-induced greening pattern) over a two-decade span from 2001 to 2021. As for vegetation resilience, Fig. 3i reveals a generally positive trend: the critical EVI_c has increased from 0.17 to 0.19 over the last two decades. This is indicative of a system shifting towards a more resilient tipping state. However, it is noteworthy that a minor decline in EVI_c was observed during 2017-2021, likely due to a period of reduced rainfall, aligning with observed climatic patterns. Fig. 3g highlights a macroscopic trend of expanding vegetation cover—often referred to as ‘greening’—particularly in the Sahel’s northern regions from 2001 to 2021.

Moreover, in Fig. 3j, the empirical value for the fragment size distribution exponent ( $\tau$ ) closely matches theoretical expectations for most of this period, suggesting that the vegetation in the Sahel is operating in a critical state (refer to Figs. S12–S17 for details). Once again, the one noticeable deviation occurs in 2017-2021, which can be attributed to significantly low rainfall levels during that period (refer to Fig. S18). Further corroborating this, the fractal dimensions for mass and external perimeter also align with theoretical values, as shown in Fig. S19. These results support the notion that Sahel’s vegetation is near a state of SOC, evidenced by the critical exponents closely matching their theoretical counterparts.

II.3 Origin of criticality

To comprehend why the QTP and the Sahel region appear to be on the brink of significant ecological changes, or ‘tipping points’, we explore two factors: correlation length ( $\xi$ ) and susceptibility ( $\chi$ ), as detailed in Section (IV.8). The correlation length measures the extent to which spatial fluctuations in a system are correlated, while susceptibility indicates the system’s sensitivity to external influences, like less rainfall or hotter temperatures. We observe that both $\chi$ and $\xi$ peak at EVI_c, as depicted in Figs. 3e and k. This indicates that when the system approaches the critical threshold, it becomes more sensitive and spatially correlated, making it more vulnerable to disturbances.

Furthermore, we find that both $\chi$ and $\xi$ diverge at EVI_c when extrapolated to an infinitely large system size, as illustrated in Figs. 3f and 3l. The divergence of correlation length and susceptibility at EVI_c indicates that the presence of inherent global correlation of the spatio-temporal dynamical processes within the greenness system [74]. This global correlation plays a vital role in the emergence of tipping criticality in the QTP and Sahel. In other words, the increased spatial correlation and sensitivity of the system near the tipping point which is close to EVI_c=0.2, lead to stronger responses to external influences, potentially triggering abrupt transitions.

II.4 Optimal enhancing resilience model

In ecology, maintaining a robust and well-connected giant ecological component (EVI $>$ 0.2) is pivotal as it facilitates essential ecological exchanges such as energy flow, species movement, information transfer, and nutrient cycling. These processes are integral to the health, resilience, and functioning of ecosystems. Energy flow, which initiates from primary producers and traverses various trophic levels, is a fundamental ecosystem process [75]. Similarly, species exchange promotes biodiversity by allowing for migration, dispersal, and colonization [76]. Information transfer, vital for species adaptation, is enabled through a well-connected ecological network [77]. Nutrient cycling, crucial for maintaining soil fertility and productivity, also relies on a well-connected ecosystem [75].

Hence, protecting the giant ecological component and maintaining its connectivity is essential to preserve the integrity, functionality, and resilience of ecosystems. To this end, we propose an Optimal Enhancing Resilience Model (OERM) to strengthen the ecosystem resilience in the QTP and Sahel. This model utilizes the percolation method to identify critical nodes within the system that, when disrupted, can lead to collapse and fragmentation, as shown in Fig. 2. Proactively protecting these nodes is an effective way to reduce systemic vulnerability and risks. In the OERM, we focus on the protection of a single node to minimize socio-economic costs associated with protection efforts. The specific details and methodology for selecting and protecting this single node can be found in Section. IV.10. By focusing on these critical nodes, OERM aims to enhance the resilience and structural integrity of the QTP and Sahel ecosystems, thereby reducing the risk of fragmentation and system failure.

For the QTP during the period 2017-2021, we identify the critical node $i_{c}$ at the percolation threshold EVI_c that is approximately located at (36.07^∘N, 101.96^∘E), as shown in Fig. 4a and b. We proceed to select $O$ nodes in the critical region and neighboring areas and protect them. To quantify the resilience of the QTP, we define the metric $R(O)$ ,

R(O)=\frac{G_{1}^{c}(O)-G_{1}^{c}}{G_{1}^{c}},

(1)

where $G_{1}^{c}(O)$ represents the largest cluster size after protecting $O$ nodes, and $G_{1}^{c}$ denotes the largest cluster pre-protection. Specifically, $G_{1}^{c}(0)=G_{1}^{c}$ . Fig. 4c demonstrates a significant resilience boost from $R=0$ to $10.9\%$ for $O\geq 1$ , with $O=1$ emerging as the optimal point. For comparison, we introduced a Random Enhancing Resilience Model (RERM), where $O$ nodes are randomly chosen for protection. The results in Fig. 4c (orange curve) reveal that RERM is not effective at improving resilience.

Subsequently, we apply both the OERM and RERM frameworks to the Sahel using data from 2013-2017. The results, shown in Fig. 4d-f, confirm that targeted protection of a critical node, as suggested by the OERM, significantly amplifies the system’s resilience. Here, protecting just one node ( $O=1$ ), results in a $R=$ 23.0% increase in resilience.

In contrast, RERM does not deliver a comparable boost in resilience, emphasizing the utility of OERM as a pragmatic tool for environmental endeavors. This model could guide afforestation efforts by identifying areas that would most benefit from increased vegetation planting, in terms of bolstering overall system resilience. Additionally, it could aid in developing sustainable land use plans by highlighting regions where human activity could potentially push the system towards a negative tipping point. Furthermore, the OERM can assist in prioritizing areas for protection or restoration by singling out those crucial to the ecosystem’s overall well-being and robustness.

III Summary and Discussion

The present study explored the greening-induced positive tipping point of the Qinghai-Tibetan Plateau (QTP) and Sahel as well as their potential implications. Our findings unveiled that both the QTP and Sahel are on the verge of a tipping point, specifically, our calculated percolation threshold EVI_c approximates an empirically determined ecological critical state of roughly 0.2. At this EVI_c threshold, the vegetation architecture adheres to perfect power laws with consistent exponents, encompassing $\tau$ , $D_{f}$ , and $D_{e}$ . Importantly, our study discovered that their vegetation undergoes ongoing improvement in system resilience and sustainability, characterizing them as favorable tipping elements. We also propose an Optimal Resilience Enhancement Model. This model aids in boosting vegetation resilience with minimal socio-economic repercussions and presents a practical strategy to preserve the ‘giant’ ecological component. This emphasizes the need for adaptability and transformation in maintaining the resilience of ecosystems [78].

The greening-induced positive tipping points propose several potential consequences and implications for both the ecosystem and human society. Firstly, in terms of ecosystem functioning, it has the potential to enhance various aspects such as nutrient cycling, water, and energy balance, and habitat availability for wildlife. These alterations contribute to the overall resilience of the ecosystem, enabling it to adapt to future environmental changes [79]. Secondly, it can also influence the provision of essential ecosystem services. For instance, it can impact carbon sequestration [80], water regulation, and soil conservation. Changes in these services have far-reaching implications, particularly for local communities that rely on them for their livelihoods. Furthermore, it can have effects on climate change feedbacks. The increased vegetation cover resulting from the tipping points can bring about alterations in the regional climate by influencing surface albedo, evapotranspiration, and the overall energy balance [81]. These changes can have cascading effects on global climate patterns, potentially mitigating risks induced by climatic changes.

Despite these promising findings, concerns persist regarding the future prospects of the QTP and Sahel vegetation due to potential changes in climate conditions and continued warming [82]. For example, the ongoing loss of glacier water in QTP is particularly worrisome as it may result in water shortage in the future, which could disrupt the stability of vegetation growth [15]. Additionally, the increasing frequency of vegetation fires poses a significant threat to ecological stability, with young trees being particularly vulnerable to high-severity fires attributed to persistent warming [83]. Notably, fire regimes are not only influenced by climate change but could also lead to net carbon emissions from terrestrial carbon reservoirs.

Our findings shed light on the positive tipping phenomenon in vegetation resilience amidst climate change and human impact, emphasizing the urgent need for robust strategies to strengthen ecosystem resilience and promote sustainable development. Employing a percolation-based framework and resilience analytics, our study provides guidance for policymakers and environmental managers in preserving ecosystems and protecting the socioeconomic welfare of future generations.

Refer to caption — Figure 1: A schematic representation of the complex conditions within the Qinghai-Tibetan Plateau (QTP) and an illustration of the percolation-based framework. a, The QTP, as a key component of Earth’s climate system, dynamically interacts with various climatic phenomena. Arrows in different colours, along with descriptions, indicate interactions with the westerly jet, Indian Summer Monsoon (ISM), East Asian Summer Monsoon, Arctic Oscillation, El Niño-Southern Oscillation (ENSO), dust storms, and the Amazon rainforest area. b, A schematic representation of the QTP’s biospheric, cryospheric, and hydrological conditions, highlights the intricate interplay between these elements. c, Photos captured during June-July of 2021 and 2022 ( $\copyright$ 2023 Jingfang Fan. All rights reserved) illustrate the diverse vegetation conditions observed in the QTP, providing visual examples of the area’s vegetation. d, The corresponding Enhanced Vegetation Index (EVI) values to various vegetation conditions shown in c are represented by the heights of the bars. Taller bars indicate higher EVI values, indicating denser vegetation. e, The spatial distribution of fragments derived from the percolation-based framework, generated from the EVI configuration shown in Panel d, is displayed. The gray background represents the nodes that were removed (EVI $<$ EVI_c), whereas coloured bars represent fragments. Nodes belonging to the same fragment share the same colour, visually demarcating clusters of vegetation.

\contcaption

a, Spatial distribution of the average EVI across the QTP over 2017-2021. EVI is an indicator of an area’s greenness, with higher values indicating denser, healthier vegetation. b, Evolution of the relative sizes of the largest (denoted by $G_{1}$ , in green) and second largest (denoted by $G_{2}$ , in blue) fragments in relation to EVI. The vertical dashed line indicates the percolation threshold, EVI_c. c, Fragmentation patterns of greenness at EVI_c. Different colours depict distinct fragments, with fragments of over 100 nodes displayed. The largest fragment is illustrated in green and the second largest in blue. d, Observed (green dots) and null model (gray histogram) distribution of fragment sizes at EVI_c. The solid line represents the power-law fit with an exponent $\tau$ $\approx$ 2.04 ( $R^{2}$ =0.996). Statistics incorporate approximately 87,000 fragments. e, Mass ( $M$ , green dots) and external perimeter ( $L$ , blue dots) for the largest fragment at various rescaled box renormalization factors ( $b=2^{0},\ 2^{1},\ 2^{2},\ 2^{3},\ 2^{4},\ 2^{5}$ ). Power-law fits with exponents $D_{f}\approx 1.93$ ( $R^{2}$ = 0.904) and $D_{e}\approx 1.40$ ( $R^{2}$ = 0.971) are denoted by the green and blue dashed lines, respectively. f-j, The same as a-e, but for the Sahel region over 2013-2017.

\contcaption

a, The spatial pattern of significant EVI change trends ( $p<0.05$ ) from 2001 to 2021, demonstrating areas of remarkable vegetation change. b, Precipitation-induced greening patterns from 2001 to 2021. The cross-correlation coefficients (according to Eq. (12)) between precipitation and EVI series are computed, revealing regions with significant correlation ( $p<0.05$ ). The colours denote the corresponding cross-correlation coefficients. EVI series in the northern arid regions of the QTP display high sensitivity to precipitation variations. However, in the southeast regions characterized by high precipitation and lush vegetation, increased precipitation is linked to decreased vegetation sensitivity (refer to Fig. S11b). c, A comparison of EVI_c (green dots) and the average temperature $T$ in June-August (blue dots) as functions of time. d, The temporal evolution of the critical exponent $\tau$ over the past two decades (2001-2021). The green shade represents the error bars. The horizontal black dashed line indicates the theoretical value of $\tau\approx 2.05$ as predicted by classical 2D percolation theory. Each point is calculated within a sliding window of 5 years. e, The correlation length $\xi$ (blue dots) and susceptibility $\chi$ (green dots) with respect to EVI. The vertical black dashed line indicates EVI_c. f, The behavior of $\xi$ (blue dots) and $\chi$ (green squares) as functions of $1/b$ , where $b$ represents the rescaled box renormalization factor. The best-fit lines for $\xi$ and $\chi$ exhibit slopes of $k_{\xi}=1.02$ ( $R^{2}=0.998$ ) and $k_{\chi}=2.05$ ( $R^{2}=0.901$ ), respectively. Panels g-l, The same as a-f, specifically for the Sahel region.

IV Data and Methods

IV.1 Data

In this study, we utilize the Enhanced Vegetation Index (EVI) data derived from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite imagery. EVI is a remote sensing metric that quantifies the greenness or density of vegetation on the Earth’s surface. It considers the reflectance of visible and near-infrared light to optimize the vegetation signal while minimizing atmospheric and background noise. Compared to other vegetation indices like the Normalized Difference Vegetation Index (NDVI), EVI is more effective in areas with dense vegetation or under varying atmospheric conditions. EVI values range from -1 to 1, with higher values indicating healthier and more vigorous vegetation [56]. We obtained the spatial composites provided at 16-day temporal resolution and 250-meter spatial resolution (MOD13C1 Version 6; https://lpdaac.usgs.gov/products/mod13q1v006/).

Additionally, we incorporate monthly surface air temperature reanalysis data ( $0.1^{\circ}\times 0.1^{\circ}$ ) from ERA5 Land (https://doi.org/10.24381/cds.68d2bb30) and monthly precipitation accumulation reanalysis data ( $0.04^{\circ}\times 0.04^{\circ}$ ) from TerraClimate [84] (https://www.climatologylab.org/terraclimate.html).

These data were preprocessed and downloaded on the Google Earth Engine (GEE) cloud platform [85]. And, we focus on the Qinghai-Tibetan Plateau (QTP) area (https://doi.org/10.3974/geodb.2014.01.12.V1, accessed on 31 July 2022) and Sahel area. The Sahel region selects areas with EVI ranging between 0.1 and 0.4, corresponding to the shapefile provided in the attachment. We took summer (July-September) as the focused season. The summer months represent the peak growing season on the QTP/Sahel, making it crucial for studying vegetation patterns, and observing the highest proportion of vegetation cover in the region.

IV.2 Data filtering

To minimize the effects of annual fluctuations in our study, we employ a sliding-window average approach. We use a 5-year sliding window that progresses by 2 years at a time, resulting in adjacent windows having a temporal overlap of 3 years. This approach calculates the temporal average of EVI values pixel-by-pixel within each sliding window. By smoothing out short-term variations and noise, the sliding-window average provides a more stable and reliable representation of the underlying spatial patterns and trends in the data. The overlapping windows ensure a continuous and smooth analysis of the temporal changes, minimizing the potential for abrupt shifts or artificial breakpoints in the time series. Each window is defined by the start year of the corresponding period.

IV.3 Greenness Percolation Model

In our Greenness Percolation model, we represent the EVI image as a 2D lattice network, where each valid pixel in the image corresponds to a node in the network. Valid pixels are those located within the boundary of the QTP (Sahel) and have non-null EVI values. Each image contains approximately 48.5 (55.2) million valid pixels, providing a high enough spatial resolution suitable for percolation analysis.

The percolation procedure follows these steps: starting from a fully occupied lattice network, we rank all the nodes based on their EVI values in ascending order. We then remove nodes one by one, starting with the node having the lowest EVI value and progressing to higher-ranked nodes. During the removal process, we identify the fragments using classical percolation theory [86, 87]. A fragment is defined as a subset of nearest neighbour nodes, where each node in the subset is connected to at least one other node. We employ the efficient Newman–Ziff algorithm [57] to detect the different fragments at each step.

It’s worth mentioning that our model incorporates free boundary conditions. Each node interacts with four neighbors, excluding the nodes adjacent to invalid pixels or the boundaries of the QTP/Sahel. In percolation theory, the relative size of the largest fragment, often referred to as the ‘giant’ fragment, is typically identified as the order parameter [86]. Here, the ‘giant’ fragment $G_{1}$ is defined as [88],

G_{1}=\frac{\max\left[\left(\sum_{i\in\mathcal{H}_{1}}g_{i}\right),\left(\sum_{i\in\mathcal{H}_{2}}g_{i}\right),...,\left(\sum_{i\in\mathcal{H}_{m}}g_{i}\right),...,\right]}{\sum_{i=1}^{N}g_{i}},

(2)

where $\mathcal{H}_{m}$ denotes a series of disjoint fragments, $g_{i}$ stands for the size/area of node $i$ , and $N$ is the total number of nodes. Since most of all nodes in the QTP are located within $[30^{\circ}N,40^{\circ}N]$ , and most of all nodes in the Sahel region are located within $[10^{\circ}N,20^{\circ}N]$ , We assume their sizes are approximately equal, i.e., $g_{i}=1$ ; However, if the latitude spans larger scales, we use then $g_{i}=\cos{(\phi_{i}})$ , where $\phi_{i}$ is the latitude of node $i$ [88]. We have also performed the same analysis by using $\phi_{i}$ , and the results do not affect our conclusion. Throughout the percolation simulation, we measure the order parameter $G_{1}(P)$ at each time step $P$ and calculate the largest one-step gap $\Delta_{c}$ as defined by Equation (3),

\Delta_{c}\equiv\max_{P}\left[G_{1}(P-1)-G_{1}(P)\right].

(3)

Equation (3) is used to determine the percolation threshold [58]. We define the time step with the largest jump $\Delta_{c}$ as $P_{c}$ . The EVI value of the critical node at $P_{c}$ is defined as the percolation threshold EVI_c, at which the greenness state abruptly shifts from cohesiveness to fragmentation.

IV.4 Power-law fitting for the fragment-size distribution

In our analysis, we calculated the fragment-size distribution at the EVI_c. We used the maximum-likelihood fitting method [89] to fit the frequency of fragment sizes to a power-law distribution according to Eq. (4). The fitting was performed using a Python package [90] available online,

n_{s}\sim s^{-\tau}

(4)

where $n_{s}$ is the number of fragments of size $s$ . We assess the Goodness of fit using the Kolmogorov-Smirnov statistic and likelihood ratios. In the main body of our paper, we set the maximum size of the power-law fitted distributions $s_{\mathrm{max}}$ to 10,000. Additionally, we test alternative maximum sizes, 5,000, 6,000, 7,000, 8,000, 9,000, and 11,000 for the QTP and 9000, 11,000, 12,000, 13,000, 14000, and 15,000 for the Sahel. The results show robustness and consistency, as demonstrated in Figs. S20 and S21, supporting the findings and conclusions of the study. These consistent results provide confidence in the methodology used in the research.

IV.5 Box renormalization

Box renormalization approach was employed to quantify the self-similarity features of the largest fragment at the percolation threshold [61]. This procedure involves systematically reducing the resolution of a large lattice grid by eliminating fluctuations on scales smaller than a given length scale, $b$ . The goal is to re-scale the lattice grid into a sequence of smaller ones, enabling the study of properties such as fractal behaviour and finite-size effects. The box renormalization operation involves selecting a box size $b$ (e.g., $b=2^{1},\ 2^{2},\ 2^{3},\ 2^{4},\ 2^{5}$ , …), and re-grouping the lattice grid into $b\times b$ boxes. Each box is replaced with a super-node, where the super-node’s value (EVI) is the average of all the nodes in the corresponding box. By applying box renormalization with a sequence of re-scale factors $b$ , a sequence of re-scaled lattice grids is obtained. The percolation procedure can then be re-applied to these re-scaled lattice grids to examine properties such as fractal behaviour and finite-size effects. For consistency, the re-scale factor $b$ for the original lattice grid is defined as 1.

IV.6 Fractal dimensions

The fractal concept, first introduced by Mandelbrot [91], was applied to percolation by Stanley [92] to describe the complexity of fractal substructures within the giant fragment at the percolation threshold. At this threshold, the giant fragment exhibits self-similarity on all length scales and can be considered a fractal. Here, we focus on two properties: the mass and external perimeter of the giant fragment at the percolation threshold.

The mass of the largest fragment is defined as the number of nodes it contains, while the accessible external perimeter, introduced by Aharony [64], is defined as the number of nodes on the external boundary. The external boundary can be determined by probing (adsorbent) nodes that are adjacent to the largest fragment’s exterior.

To estimate the mass and external perimeter fractal dimensions ( $D_{f}$ and $D_{e}$ , respectively), we examine the power-law relationship between these quantities and the re-scale factors $b$ as given by Eqs. (5) and (6) where the mass and external perimeter of the giant fragment at the percolation threshold are represented by $M(b)$ and $L(b)$ respectively,

M(b)\sim b^{-D_{f}},

(5)

L(b)\sim b^{-D_{e}}.

(6)

The values $D_{f}$ and $D_{e}$ indicate the presence of self-similar fractal patterns and provide insight into the complexity of the fractal substructures within the fragment.

To obtain the mass and external perimeter of the giant fragment at EVI_c, we renormalize the EVI image for a sequence of re-scale factors ( $b=2^{1},\ 2^{2},\ 2^{3},\ 2^{4},\ 2^{5}$ ). The EVI_c value is defined as the EVI corresponding to the largest jump of the giant fragment within the selected EVI range of 0.15-0.3, which is mainly characterized by the alpine meadow. Finally, we calculate the mass and external perimeter fractal dimensions of the largest fragment at EVI_c using Eqs. (5) and (6) for $b=2^{0},\ 2^{1},\ 2^{2},\ 2^{3},\ 2^{4},\ 2^{5}$ .

IV.7 Null model

In this study, we create a null model to serve as a comparison standard. The null model is generated by randomly shuffling the EVI values to create 100 new samples. The percolation procedure was then applied to each of these new samples. In these new samples, the original spatial associations are destroyed, and percolation is performed on a random spatial profile, corresponding to the uncorrelated site percolation class. According to percolation theory, the null model is expected to exhibit critical phenomena at the percolation threshold. Therefore, it serves as a comparison standard for critical features. By analysing the behaviour of the null model, we can better understand the significance of spatial associations in the original data and assess their impact on the observed critical phenomena and tipping points in the QTP and Sahel ecosystems.

IV.8 Correlation length and susceptibility

The correlation length ( $\xi$ ) in percolation theory represents the average distance between two sites that belong to the same fragment. Equation (7) defines the correlation length as,

{}\xi^{2}=\frac{2\sum^{\prime}_{s}R_{s}^{2}s^{2}n_{s}}{\sum^{\prime}_{s}s^{2}n_{s}},

(7)

where $n_{s}$ denotes the number of fragments of size $s$ , $R_{s}$ is the average gyration radius of fragments of size $s$ , and the prime on the sums indicates the exclusion of the largest fragment in each measurement.

The gyration radius of a given fragment $a$ of size $s$ is defined by,

R^{2}_{a}=\frac{1}{s^{2}}\sum_{i=1}^{s}\sum_{j>i}^{s}|\mathbf{r}^{a}_{i}-\mathbf{r}^{a}_{j}|^{2},

(8)

The average gyration radius of all the fragments of size $s$ can be calculated using Equation (9),

{}R_{s}^{2}=\frac{1}{s^{2}n_{s}}\sum_{a=1}^{n_{s}}\sum_{i=1}^{s}\sum_{j>i}^{s}|\mathbf{r}_{i}^{a}-\mathbf{r}_{j}^{a}|^{2},

(9)

where $|\mathbf{r}_{i}^{a}-\mathbf{r}_{j}^{a}|$ denotes the Euclidean distance between the $i$ th and $j$ th node in fragment $a$ . The susceptibility $\chi$ is defined as:

\chi=\frac{\sum^{\prime}_{s}s^{2}n_{s}}{\sum^{\prime}_{s}sn_{s}},

(10)

where the prime on the sums indicates the exclusion of the largest fragment in each measurement. The susceptibility represents the system’s sensitivity to changes in the percolation threshold, and it is used to characterize the critical phenomena.

To demonstrate that the correlation length $\xi$ and susceptibility $\chi$ diverge at EVI_c as the size of the empirical system tends to approach infinity, we investigate the finite-size effects using box renormalization. We renormalized the EVI image for a sequence of re-scale factors $b=2^{1},\ 2^{2},\ 2^{3},\ 2^{4},\ 2^{5}$ and applied the percolation procedure to these re-scaled EVI images. Then we calculated the correlation length $\xi$ and susceptibility $\chi$ at EVI_c for $b=2^{0},\ 2^{1},\ 2^{2},\ 2^{3},\ 2^{4},\ 2^{5}$ , examining their behaviour as the system size is extrapolated to infinity. By analysing the behaviour of $\xi$ and $\chi$ as a function of the re-scale factor $b$ , we can determine whether these quantities diverge at the critical EVI value (EVI_c), indicating the existence of a critical point and a phase transition at the percolation threshold.

IV.9 Cross-correlation

The cross-correlation with time lag $\sigma$ between the EVI series $E(y)$ and the precipitation series $Q(y)$ are defined by:

C(\sigma)=\frac{\langle Q(y)E(y+\sigma)\rangle-\langle Q(y)\rangle\langle E(y+\sigma)\rangle}{\sqrt{\left\langle\left(Q(y)-\langle Q(y)\rangle\right)^{2}\right\rangle}\times\sqrt{\left\langle\left(E(y+\sigma)-\langle E(y+\sigma)\rangle\right)^{2}\right\rangle}}

(11)

where $\sigma\in[-\sigma_{\mathrm{max}},\sigma_{\mathrm{max}}]$ is the time lag, with $\sigma_{\mathrm{max}}=4$ year. Therefore, we can achieve $2\sigma_{\mathrm{max}}+1$ different cross-correlation values. The maximum absolute value of the cross-correlation function was identified, and the corresponding time lag is denoted as $\sigma_{0}$ . In this study, we redefine cross-correlation $C$ as,

C=C(\sigma_{0}).

(12)

IV.10 Optimal enhancing resilience model

We propose the Optimal Enhancing Resilience Model (OERM) as a solution to further improve the stability of the system. The OERM is based on the snapshot of fragments at EVI_c where the critical node $i_{c}$ , is identified as the node connecting the largest and second-largest fragments at EVI_c. The total number of protected nodes is denoted by $O$ (with a max value of $90$ in this study). When the protection steps $O=0$ , no nodes are specifically protected. When the protection steps $O=1$ , the critical node $i_{c}$ is protected, and its EVI value is changed to $1.0$ . Subsequently, unoccupied nodes at EVI_c are gradually protected based on their distance from the critical node $i_{c}$ , following a near-to-far approach. If multiple nodes have the same distance, the one with a larger EVI is protected first. The EVI of the protected nodes is also changed to $1.0$ .

The protection resilience $R$ , is defined as the ratio of the fragment size increase after protection to the fragment size before protection. It is measured using Eq. (1). To compare the effectiveness of OERM, we propose the Random Enhancing Resilience Model (RERM), where an equal number of protected nodes are randomly selected. By comparing the performance of OERM and RERM, we aim to evaluate the effectiveness of targeted protection strategies based on critical nodes and their surrounding areas. This comparison would help develop more efficient and cost-effective measures for enhancing ecosystem resilience and mitigating the impacts of climate change.

Data availability

The data used in Figs. 2–4 are available as Source Data, and other data supporting the plots within the paper and findings of the study can be obtained from the corresponding author upon request.

code availability

The Python codes used for the analysis are available on GitHub (https://github.com/fanjingfang/RGPTP).

Acknowledgments

We acknowledge the support of the National Natural Science Foundation of China (Grant No. 12275020, 12135003, 12205025) and the Ministry of Science and Technology of China (2019QZKK0906). D.C. is supported by the Swedish strategic research area MERGE. J.K. received support from the Germany BMBF grant 01LP1902A.

Author Contributions

J.F. designed the research. Y.S. and J.F. performed the analysis and prepared the manuscript, Y.S., T.L., S.W., J.M., Y.Z., S.Y., X.C., D.C., J.K., S.H., H.J.S. and J.F. discussed results and contributed to writing the manuscript. J.F. led the writing of the manuscript.

Additional information

Supplementary Information is available in the online version of the paper.

Competing interests

The authors declare no competing interests.

References

[1] Lenton, T. M. et al. Tipping elements in the Earth’s climate system. Proceedings of the National Academy of Sciences 105, 1786–1793 (2008).
[2] Lenton, T. M. et al. Climate tipping points—too risky to bet against. Nature 575, 592–595 (2019).
[3] Armstrong McKay, D. I. et al. Exceeding 1.5°C global warming could trigger multiple climate tipping points. Science 377, eabn7950 (2022).
[4] May, R. M., Levin, S. A. & Sugihara, G. Ecology for bankers. Nature 451, 893–894 (2008).
[5] Council, N. R. New Directions for Understanding Systemic Risk: A Report on a Conference Cosponsored by the Federal Reserve Bank of New York and the National Academy of Sciences (The National Academies Press, Washington, DC, 2007).
[6] Hirota, M., Holmgren, M., Van Nes, E. H. & Scheffer, M. Global resilience of tropical forest and savanna to critical transitions. Science 334, 232–235 (2011).
[7] Scheffer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. Catastrophic shifts in ecosystems. Nature 413, 591–596 (2001).
[8] Kéfi, S. et al. Spatial vegetation patterns and imminent desertification in mediterranean arid ecosystems. Nature 449, 213–217 (2007).
[9] Drake, J. M. & Griffen, B. D. Early warning signals of extinction in deteriorating environments. Nature 467, 456–459 (2010).
[10] Lever, J. J., van Nes, E. H., Scheffer, M. & Bascompte, J. The sudden collapse of pollinator communities. Ecology Letters 17, 350–359 (2014).
[11] Jiang, J. et al. Predicting tipping points in mutualistic networks through dimension reduction. Proceedings of the National Academy of Sciences 115, E639–E647 (2018).
[12] Centola, D., Becker, J., Brackbill, D. & Baronchelli, A. Experimental evidence for tipping points in social convention. Science 360, 1116–1119 (2018).
[13] Otto, I. M. et al. Social tipping dynamics for stabilizing Earth’s climate by 2050. Proceedings of the National Academy of Sciences 117, 2354–2365 (2020).
[14] Yao, T. et al. Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings. Nature Climate Change 2, 663–667 (2012).
[15] Immerzeel, W. W. et al. Importance and vulnerability of the world’s water towers. Nature 577, 364–369 (2019).
[16] Yao, T. et al. Recent third pole’s rapid warming accompanies cryospheric melt and water cycle intensification and interactions between monsoon and environment: Multidisciplinary approach with observations, modeling, and analysis. Bulletin of the American Meteorological Society 100, 423–444 (2019).
[17] Xu, X., Lu, C., Shi, X. & Gao, S. World water tower: An atmospheric perspective. Geophysical Research Letters 35 (2008).
[18] Yang, K. et al. Response of hydrological cycle to recent climate changes in the Tibetan Plateau. Climatic Change 109, 517–534 (2011).
[19] Shen, M. et al. Evaporative cooling over the Tibetan Plateau induced by vegetation growth. Proceedings of the National Academy of Sciences 112, 9299–9304 (2015).
[20] Yao, T. et al. Third pole environment (tpe). Environmental Development 3, 52–64 (2012).
[21] Cai, Y. et al. The Holocene Indian monsoon variability over the southern Tibetan Plateau and its teleconnections. Earth and Planetary Science Letters 335-336, 135–144 (2012).
[22] Sato, T. & Kimura, F. How does the Tibetan Plateau affect the transition of indian monsoon rainfall? Monthly Weather Review 135, 2006–2015 (2007).
[23] Sills, J., Immerzeel, W. W. & Bierkens, M. F. P. Asian water towers: More on monsoons–response. Science 330, 585–585 (2010).
[24] Yang, K. et al. Recent climate changes over the Tibetan Plateau and their impacts on energy and water cycle: A review. Global and Planetary Change 112, 79–91 (2014).
[25] Duan, A., Wu, G., Liu, Y., Ma, Y. & Zhao, P. Weather and climate effects of the Tibetan Plateau. Advances in Atmospheric Sciences 29, 978–992 (2012).
[26] Liu, T. et al. Teleconnections among tipping elements in the Earth system. Nature Climate Change (2023).
[27] Bi, M. et al. Effects of arctic sea ice in autumn on extreme cold events over the Tibetan Plateau in the following winter: possible mechanisms. Climate Dynamics 58, 2281–2292 (2022).
[28] Shaman, J. & Tziperman, E. The effect of ENSO on Tibetan Plateau snow depth: A stationary wave teleconnection mechanism and implications for the south Asian monsoons. Journal of Climate 18, 2067–2079 (2005).
[29] Wu, Z., Li, J., Jiang, Z. & Ma, T. Modulation of the Tibetan Plateau snow cover on the ENSO teleconnections: From the east Asian summer monsoon perspective. Journal of Climate 25, 2481–2489 (2012).
[30] Huang, J. et al. Summer dust aerosols detected from CALIPSO over the Tibetan Plateau. Geophysical Research Letters 34 (2007).
[31] Chen, D. et al. Assessment of past, present and future environmental changes on the Tibetan Plateau. Chinese Science Bulletin 60, 3025–3035 (2015).
[32] Zhu, Z. et al. Greening of the Earth and its drivers. Nature Climate Change 6, 791–795 (2016).
[33] Zhang, G., Zhang, Y., Dong, J. & Xiao, X. Green-up dates in the Tibetan Plateau have continuously advanced from 1982 to 2011. Proceedings of the National Academy of Sciences 110, 4309–4314 (2013).
[34] Shen, M. et al. Increasing altitudinal gradient of spring vegetation phenology during the last decade on the Qinghai-Tibetan Plateau. Agricultural and Forest Meteorology 189-190, 71–80 (2014).
[35] Liang, E. et al. Species interactions slow warming-induced upward shifts of treelines on the Tibetan Plateau. Proceedings of the National Academy of Sciences 113, 4380–4385 (2016).
[36] Gou, X. et al. Patterns and dynamics of tree-line response to climate change in the eastern Qilian Mountains, northwestern China. Dendrochronologia 30, 121–126 (2012).
[37] Chen, Y., Feng, J., Yuan, X. & Zhu, B. Effects of warming on carbon and nitrogen cycling in alpine grassland ecosystems on the Tibetan Plateau: A meta-analysis. Geoderma 370, 114363 (2020).
[38] Lee, E., He, Y., Zhou, M. & Liang, J. Potential feedback of recent vegetation changes on summer rainfall in the Sahel. Physical Geography 36, 449–470 (2015).
[39] Bégué, A., Vintrou, E., Ruelland, D., Claden, M. & Dessay, N. Can a 25-year trend in Soudano-Sahelian vegetation dynamics be interpreted in terms of land use change? A remote sensing approach. Global Environmental Change 21, 413–420 (2011).
[40] Le Houerou, H. N. The Rangelands of the Sahel. Journal of Range Management 33, 41–46 (1980). eprint 3898226.
[41] Diffenbaugh, N. S. & Giorgi, F. Climate change hotspots in the cmip5 global climate model ensemble. Climatic Change 114, 813–822 (2012).
[42] Bathiany, S., Claussen, M. & Brovkin, V. CO2-Induced Sahel Greening in Three CMIP5 Earth System Models. Journal of Climate 27, 7163–7184 (2014).
[43] Ouedraogo, I., Runge, J., Eisenberg, J., Barron, J. & Sawadogo-Kaboré, S. The Re-Greening of the Sahel: Natural Cyclicity or Human-Induced Change? Land 3, 1075–1090 (2014).
[44] Armstrong, E., Tallavaara, M., Hopcroft, P. O. & Valdes, P. J. North African humid periods over the past 800,000 years. Nature Communications 14, 5549 (2023). Number: 1 Publisher: Nature Publishing Group.
[45] Pausata, F. S. R. et al. The Greening of the Sahara: Past Changes and Future Implications. One Earth 2, 235–250 (2020).
[46] Nicholson, S. On the question of the “recovery” of the rains in the west african sahel. Journal of Arid Environments 63, 615–641 (2005).
[47] Brovkin, V., Claussen, M., Petoukhov, V. & Ganopolski, A. On the stability of the atmosphere-vegetation system in the Sahara/Sahel region. Journal of Geophysical Research: Atmospheres 103, 31613–31624 (1998).
[48] Vamborg, F. S. E., Brovkin, V. & Claussen, M. Background albedo dynamics improve simulated precipitation variability in the Sahel region. Earth System Dynamics 5, 89–101 (2014).
[49] Xue, Y. et al. The Sahelian Climate. In Kabat, P. et al. (eds.) Vegetation, Water, Humans and the Climate: A New Perspective on an Interactive System, Global Change — The IGBP Series, 59–77 (Springer, Berlin, Heidelberg, 2004).
[50] Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53–59 (2009).
[51] Scheffer, M. et al. Anticipating critical transitions. Science 338, 344–348 (2012).
[52] Livina, V. N. & Lenton, T. M. A modified method for detecting incipient bifurcations in a dynamical system. Geophysical Research Letters 34 (2007).
[53] Bak, P., Tang, C. & Wiesenfeld, K. Self-organized criticality: An explanation of the 1/f noise. Physical Review Letters 59, 381–384 (1987).
[54] Taubert, F. et al. Global patterns of tropical forest fragmentation. Nature 554, 519–522 (2018).
[55] Bunde, A. & Havlin, S. Fractals and disordered systems (Springer Science & Business Media, 2012).
[56] Huete, A. et al. Overview of the radiometric and biophysical performance of the modis vegetation indices. Remote Sensing of Environment 83, 195–213 (2002).
[57] Newman, M. E. J. & Ziff, R. M. Efficient Monte Carlo algorithm and high-precision results for percolation. Physical Review Letters 85, 4104–4107 (2000).
[58] Fan, J. et al. Universal gap scaling in percolation. Nature Physics 16, 455–461 (2020).
[59] Margolina, A., Herrmann, H. J. & Stauffer, D. Size of largest and second largest cluster in random percolation. Physics Letters A 93, 73–75 (1982).
[60] Liu, H. Q. & Huete, A. A feedback based modification of the NDVI to minimize canopy background and atmospheric noise. IEEE Transactions on Geoscience and Remote Sensing 33, 457–465 (1995).
[61] Stauffer, D. & Aharony, A. Introduction to percolation theory (Taylor & Francis, 2018).
[62] ben Avraham, D. & Havlin, S. Diffusion and Reactions in Fractals and Disordered Systems (Cambridge University Press, New York, 2005).
[63] Nijs, M. P. M. d. A relation between the temperature exponents of the eight-vertex and q-state Potts model. Journal of Physics A: Mathematical and General 12, 1857 (1979).
[64] Grossman, T. & Aharony, A. Accessible external perimeters of percolation clusters. Journal of Physics A: Mathematical and General 20, L1193 (1987).
[65] Cohen, R., Erez, K., ben Avraham, D. & Havlin, S. Resilience of the Internet to random breakdowns. Physical Review Letters 85, 4626–4628 (2000).
[66] Chen, H. et al. The impacts of climate change and human activities on biogeochemical cycles on the Qinghai-Tibetan Plateau. Global Change Biology 19, 2940–2955 (2013).
[67] Tan, K. et al. Application of the ORCHIDEE global vegetation model to evaluate biomass and soil carbon stocks of Qinghai-Tibetan grasslands. Global Biogeochem. Cycles 24 (2010).
[68] Wang, S. et al. Effects of warming and grazing on soil n availability, species composition, and ANPP in an alpine meadow. Ecology 93, 2365–2376 (2012).
[69] Shen, M., Piao, S., Cong, N., Zhang, G. & Jassens, I. A. Precipitation impacts on vegetation spring phenology on the Tibetan Plateau. Global Change Biology 21, 3647–3656 (2015).
[70] Delang, C. O. & Yuan, Z. China’s grain for green program (Springer, 2016).
[71] Deng, L. et al. Past and future carbon sequestration benefits of china’s grain for green program. Global Environmental Change 47, 13–20 (2017).
[72] Deng, L., Liu, G.-b. & Shangguan, Z.-p. Land-use conversion and changing soil carbon stocks in china’s ‘grain-for-green’ program: a synthesis. Glob Change Biol 20, 3544–3556 (2014).
[73] Li, C. et al. Drivers and impacts of changes in China’s drylands. Nature Reviews Earth & Environment 2, 858–873 (2021).
[74] Jensen, H. J. What is critical about criticality: in praise of the correlation function. Journal of Physics: Complexity 2, 032002 (2021).
[75] Chapin, F. S., Matson, P. A., Mooney, H. A. & Vitousek, P. M. Principles of terrestrial ecosystem ecology (Springer, 2002).
[76] Hilty, J. A., Lidicker Jr, W. Z. & Merenlender, A. M. Corridor ecology: the science and practice of linking landscapes for biodiversity conservation (Island press, 2012).
[77] Danchin, E., Giraldeau, L.-A., Valone, T. J. & Wagner, R. H. Public Information: From Nosy Neighbors to Cultural Evolution. Science 305, 487–491 (2004).
[78] Folke, C. et al. Resilience thinking: integrating resilience, adaptability and transformability. Ecology and society 15 (2010).
[79] Chapin III, F. S. et al. Consequences of changing biodiversity. Nature 405, 234–242 (2000).
[80] Yun, H. et al. Warming and increased respiration have transformed an alpine steppe ecosystem on the tibetan plateau from a carbon dioxide sink into a source. J Geophys Res Biogeosci 127, e2021JG006406 (2022).
[81] Bonan, G. B. Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of Forests. Science 320, 1444–1449 (2008).
[82] Ehlers, T. A. et al. Past, present, and future geo-biosphere interactions on the tibetan plateau and implications for permafrost. Earth-Science Reviews 234, 104197 (2022).
[83] Bowman, D. M. J. S. et al. Vegetation fires in the Anthropocene. Nature Reviews Earth & Environment 1, 500–515 (2020).
[84] Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. Terraclimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015. Scientific Data 5, 170191 (2018).
[85] Gorelick, N. et al. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment (2017).
[86] Cohen, R. & Havlin, S. Complex networks: structure, robustness and function (Cambridge University Press, 2010).
[87] Newman, M. Networks (Oxford university press, 2018).
[88] Fan, J., Meng, J., Ashkenazy, Y., Havlin, S. & Schellnhuber, H. J. Climate network percolation reveals the expansion and weakening of the tropical component under global warming. Proceedings of the National Academy of Sciences 115, E12128–E12134 (2018).
[89] Clauset, A., Shalizi, C. R. & Newman, M. E. J. Power-law distributions in empirical data. SIAM Review 51, 661–703 (2009).
[90] Alstott, J., Bullmore, E. & Plenz, D. powerlaw: A Python package for analysis of heavy-tailed distributions. PLOS ONE 9, 1–11 (2014).
[91] Mandelbrot, B. B. The fractal geometry of nature (Times Books, San Francisco, 1982).
[92] Stanley, H. E. Cluster shapes at the percolation threshold: and effective cluster dimensionality and its connection with critical-point exponents. Journal of Physics A: Mathematical and General 10, L211 (1977).