Network dynamics of drought-induced tipping cascades in the Amazon rainforest

Tipping elements are nonlinear subsystems of the Earth system that can poten- tially abruptly and irreversibly shift if environmental change occurs. Among these tipping elements is the Amazon rainforest, which is threatened by an- thropogenic activities and increasingly frequent droughts. Here, we assess how extreme deviations from climatological rainfall regimes may cause local forest- savanna transitions that cascade through the coupled forest-climate system. We develop a dynamical network model to uncover the role of atmospheric moisture recycling in such tipping cascades. We account for the heterogeneity in critical thresholds of the forest caused by adaptation to local climatic con- ditions. Our results reveal that, despite this adaptation, increased dry-season intensity may trigger tipping events particularly in the southeastern Amazon. Moisture recycling is responsible for one-fourth of the tipping events. If the rate of climate change exceeds the adaptive capacity of some parts of the for- est, secondary effects through moisture recycling may exceed this capacity in other regions, increasing the overall risk of tipping across the Amazon rain- forest.


1
Tipping elements are nonlinear subsystems of the Earth system that can poten-4 tially abruptly and irreversibly shift if environmental change occurs. Among 5 these tipping elements is the Amazon rainforest, which is threatened by an-6 thropogenic activities and increasingly frequent droughts. Here, we assess how 7 extreme deviations from climatological rainfall regimes may cause local forest-8 savanna transitions that cascade through the coupled forest-climate system. 9 We develop a dynamical network model to uncover the role of atmospheric 10 moisture recycling in such tipping cascades. We account for the heterogeneity 11 in critical thresholds of the forest caused by adaptation to local climatic con-12 ditions. Our results reveal that, despite this adaptation, increased dry-season 13 intensity may trigger tipping events particularly in the southeastern Amazon. The Amazon rainforest is the most biodiverse terrestrial ecosystem and plays a fundamental role 20 in regulating the global climate 1,2,3 . However, human-induced impacts and climatic extremes 21 are increasingly threatening the forest's integrity and the services it provides 4,5,6 . Furthermore, 22 changes might not be gradual, but could be rather abrupt due to nonlinear interactions, as sug-23 gested by simulation studies 7,8 , data-based approaches 9,10 , conceptual models 11,12 and long-24 term experiments 13 . Parts of the Amazon rainforest may be bistable, meaning that they could 25 tip to an alternative state of low tree cover 9,10 . Indeed, the Amazon has been suggested to be 26 a tipping element in the Earth system 14 and might be at risk of approaching or exceeding its The Amazon is not a uniform forest as trees can adapt to local climatic conditions, for in-38 stance through variable rooting strategies 24,25 . This can lead to different absolute forest mortal-39 ity thresholds with respect to precipitation and drought conditions on local to regional scales. 40 Forest adaptation can therefore ensure that plants will operate close to their physiological max- 41 imum, but this creates vulnerabilities when the climate changes faster than the ecosystem can ration) and by directly re-evaporating precipitation from their leaves (interception evaporation). 48 The total amount of moisture recycling accounts for up to half of the precipitation over the 49 Amazon basin and moisture is recycled up to six times 28,29,30,31 . Thus, the rainforest depends on 50 itself, and precipitation and evapotranspiration cycles promote cascading forest growth 29 . The 51 positive interplay between the forest and regional precipitation implies that local perturbations 52 can propagate through the system via reduced moisture recycling. In other words, the Amazon 53 rainforest can be considered a network of local tipping elements that are connected via moisture 54 recycling.  Here, we integrate for the first time in a dynamical network model the tipping behaviour of 69 the Amazon forest, atmospheric moisture flows from evapotranspiration to precipitation and the 70 adaptation of the forest to annual precipitation and droughts (Fig. 1). Specifically, we combine 71 a dynamical system model to represent empirically obtained forest tipping points with regard to

88
Tipping due to drier conditions. To investigate a range of drought intensities and precipitation 89 anomalies, we study the extent of the tipped area with respect to Z-scores, which represent how 90 many standard deviations the conditions are away from the mean across 1974-2003. We find 91 a close correlation between Z MCWD and the tipped area, where a higher index reflects a larger 92 tipped area (see Fig. 2a Figure 1 | Nonlinear effects and moisture recycling network in the Amazon rainforest. a, Dynamical property of each 1˝x1˝cell of the rainforest depicted as state of rainforest cell (forest cover) versus MAP value. The state of the rainforest is limited by full forest cover (1.0) and no forest cover (-1.0). Between these two stable states, there is a tipping process as soon as the MAP value has fallen below its adaptation specific MAP crit value. Since we are focussing on drought triggered tipping events from forest to non-forest states in this study, each cell can only exist on the occupied states (brown), but not on the unoccupied states (grey). The blue arrow depicts a potential reduction in precipitation that is sufficient to trigger a tipping event in this specific cell. b, Same as in a for MCWD. c, Exemplary moisture recycling network: the rainforest cells are interconnected via a moisture recycling network due to precipitation and evapotranspiration. Through this mechanism effects of reduced tree cover can be promoted and tipping cascades are possible. d, Moisture recycling network for the hydrological year 2014 thresholded for links above 10 mm/yr to remain visibility. In the simulation results, links above 1 mm/yr are used. The dominant flow direction comes from the Atlantic ocean through easterly winds, reaches the Andes, and is then bend southward along the Andes.   We separate tipping events into primarily induced tipping events from MAP or MCWD and 105 secondary events from network effects (tipping cascades). Our model shows that between 10% 106 and 60% of the tipping is due to the cascading effects from the moisture recycling network  We also compared these results with the results of an only MAP-based normalised drought in-115 dex Z MAP analogous to Eq. 2, but find no correlation between the tipped area and the MAP 116 based index (see Supplementary Fig. 3). ell has an area of approximately`111 km 2˘. b, The additional tipped area due to network effects for each year is shown in percentage of the tipped area in panel a. This shows the effects of cascading transitions which are on the order of 10% to 60% depending on the evaluated hydrological year. The same analysis has been performed for a MAP based index, but no correlation was found (see Supplementary Fig. 3).
shows increased patterns of vulnerability (see likelihood of vulnerability in Fig. 3a). This region 120 is located in the southeastern Amazon, and caused by the combination of MCWD anomalies  Amazon region (see Fig. 4a). year to year over the study period (see Supplementary Fig. 5).

155
The network effects are especially strong close to the region of direct MCWD-induced tipping 156 and downwind from that (Fig. 3c, d). MCWD is the most important reason for tipping events 157 in the southeast, whereas MAP is not responsible for many tipping events (see Fig. 3a, b).

158
Overall, the region in the southeast is vulnerable with respect to MCWD since this region has 159 a relatively low interannual variability (standard deviation) of MCWD, while the intra-annual 160 variability (mean) MCWD value is high (see Supplementary Fig. 7c, d). nual rainfall. With a potential increase of future extreme drought events 36,37,44 , the average 173 regional climate will be drier and some parts of the rainforest might thus be set under imminent 174 risk of instability and could transgress into a less or non forest-covered state. We uncover that 175 tipping events occur most frequently in southeastern Amazon (Fig. 3). This is also the region 176 that is affected by three other factors. First, extended tipping cascades can be expected due to 177 local interaction structures and reduced downwind moisture transport (Figs. 3 and 4). Second, it 178 is also one of the regions located along the "arc of deforestation" and therefore already suffers In our study, we also find that potential future extreme drought conditions with a higher MCWD 184 anomalies show a considerably larger tipped area. Cascading tipping events are more pro-185 nounced under these circumstances (Fig. 2). These are the drought conditions that can be where k is the number of the month in the hydrological year. We make use of the actual mea- Computation of the Z-score The Z-score is used to find the ranges of future conditions that we 235 are simulating in this work. We simulate ranges from current conditions up to extreme droughts 236 that are 3.5 standard deviations away from the mean (see Fig. 2). The MCWD based Z-score is 237 computed by Here, µ MCWD and σ MCWD are the average and standard deviation of the calibration period from β px, a, bq " pσ upper´σlower q¨x a´1 p1´xq b´1 Here, we use σ upper " 3.0 and σ lower " 1.0 for the upper and lower bounds. We choose 300 a " b " 2.5 which ensures that, on average, 75% of all values lie between 1.5 to 2.5 stan-301 dard deviations and 12.5% lie between 1.0 to 1.5 or between 2.5 to 3.0 standard deviations,

302
respectively. This means that 75% lie in the central interval and 25% outside (75-25-rule). We Network of coupled nonlinear differential equations. We use a combination of nonlinear 309 differential equations together with a complex network to describe the state of the rainforest 310 cells and their interactions. We use this approach instead of a threshold approach since we 311 want to be able to account for partial changes in the state and their effects on the network.

312
For instance such changes can be critical for the tipping of cells that are not coupled directly,

316
In the differential equation approach in this work, we model the main hydrological parameters 317 and the stability of the rainforest, but no further parameters such as biotic variables. The main 318 hydrological properties are the precipitation (MAP), the MCWD and the moisture recycling.

319
Following the reasoning above, we describe the mathematical details in the remainder of this 320 section.

322
Each 1˝x1˝cell is represented by a differential equation of the form where x i stands for the state of the rainforest cell and can be interpreted as the fraction of the tree 324 cover. The shape of this function can be see in Supplementary Fig. 11. Furthermore, Eq. 6 has 325 the normal form of a saddle-node bifurcation and is a simple form of a differential equation with or completely treeless. It is not possible for a cell to have lower tree cover values than 0% or 331 values higher than full forest cover such that the state x i is limited to the interval r´1.0, 1.0s.

332
The advantage of choosing state limits of´1.0 and`1.0 is that the critical value then remains 333 analytically representable and has the specific value C crit " a 4{27 (see Supplementary Fig. 11).

334
This value is derived from the discriminant of the polynomial of Eq. 6 and more details can be

342
In our case, the rainforest cells are not independent, but interact via moisture recycling such that 343 Eq. 6 becomes Here, the entries of the critical matrix M ji p∆MAP ji , ∆MCWD ji q represent the strength of the 345 moisture recycling link between two grid cells from j to i. The state x j must be divided by 346 2 since the distance from minimum to maximum state is 2. Similar forms of the network and 20 Although both equations (Eqs. 8 and 9) are in principle not limited, we restrict them to the 359 interval r0, C crit s since the critical value for tipping of Eq. 6 is reached at C crit such that higher 360 values are not necessary to tip a certain cell.

361
Then, the critical function F crit pMAP i , MCWD i q is computed as Again, the values of Eq. 10 are restricted to the interval r0, C crit s since a state change occurs as 363 soon as the upper value of the interval, i.e. C crit , is reached. The first term of Eq. 10 is sufficient 364 to determine the critical function F crit pMAP i , MCWD i q if F crit pMAP i q or F crit pMCWD i q are 365 smaller than zero or larger than C crit . In case F crit pMAP i q and F crit pMCWD i q are larger than 366 zero, but smaller than C crit , both terms of Eq. 10 are required. The second term takes the addi- before tipping point in the sketch in Fig. 1a, b). This is an advantage of a fully dynamic model 371 such as this, while threshold-only models would not be capable of doing this.

372
An example could be that F crit pMAP 1 q " F crit pMAP 2 q " 1 2¨C crit due to respective MAP values 373 for two cells at the same time. Then it makes sense that the state of these two cells that have 374 exactly this critical value with respect to MAP is not the same in case they have a different value 375 with respect to their MCWD values. Let us assume that cell 1 has F crit pMCWD 1 q " 1 4¨C crit and 376 cell 2 has F crit pMCWD 2 q " 1 16¨C crit . Then, the second term of Eq. 10 takes these differences 377 between the cells 1 and 2 into account shifting cell 1 a bit closer to its tipping point than cell 2 378 such that the reduction effect on the respective outgoing moisture recycling links is stronger for 379 cell 1 than for cell 2.

381
Computation of the critical matrix. In analogy to Eqs. 8 and 9, we define the critical matrix 382 for MAP as where ∆MAP ji represents the difference of the mean annual precipitation arising from the 384 moisture recycling link δ ji from cell j to cell i. Thus: ∆MAP ji " ∆MAP pδ ji q " δ ji, MAP .

386
For MCWD we have where ∆MCWD ji " ∆MCWD pδ ji, MAP q is the potential increase of MCWD in response to if F crit pMAP i q ă F crit pMCWD i q.  Fig. 3a with Supplementary Fig. 12a, b, c). Thus, the 404 qualitative pattern is the same. The absolute values also show a close quantitative match within 405 their standard deviations for resolutions of 1˝x1˝or coarser (see Supplementary Fig. 13a, b).

406
The finer the resolution is, the higher the tipped area tendentially is. This is due to the fact that Nonlinear effects and moisture recycling network in the Amazon rainforest. a, Dynamical property of each 1°x1° cell of the rainforest depicted as state of rainforest cell (forest cover) versus MAP value. The state of the rainforest is limited by full forest cover (1.0) and no forest cover (-1.0). Between these two stable states, there is a tipping process as soon as the MAP value has fallen below its adaptation speci c MAPcrit value. Since we are focussing on drought triggered tipping events from forest to non-forest states in this study, each cell can only exist on the occupied states (brown), but not on the unoccupied states (grey). The blue arrow depicts a potential reduction in precipitation that is su cient to trigger a tipping event in this speci c cell. b, Same as in a for MCWD. c, Exemplary moisture recycling network: the rainforest cells are interconnected via a moisture recycling network due to precipitation and evapotranspiration. Through this mechanism effects of reduced tree cover can be promoted and tipping cascades are possible. d, Moisture recycling network for the hydrological year 2014 thresholded for links above 10 mm/yr to remain visibility. In the simulation results, links above 1 mm/yr are used. The dominant ow direction comes from the Atlantic ocean through easterly winds, reaches the Andes, and is then bend southward along the Andes. Moisture recycling links based on separate months can be found in Supplementary Figs. 1 and 2 comparing the year 2014 with the extreme drought year 2010.

Figure 2
Vulnerability of the rainforest against MCWD-based drought intensity. a, The total tipped area is shown over the course of the normalised drought index based on the MCWD Z-score. The tipped area represents the number of tipped cells in the model where each 1°x1° cell has an area of approximately (111 km2).
b, The additional tipped area due to network effects for each year is shown in percentage of the tipped area in panel a. This shows the effects of cascading transitions which are on the order of 10% to 60% depending on the evaluated hydrological year. The same analysis has been performed for a MAP based index, but no correlation was found (see Supplementary Fig. 3).

Figure 3
Vulnerable regions and tipping reason. a, The likelihood of tipping as an average over all ensemble members and all evaluated years from the hydrological years 2004 to 2014. The southeastern region is more vulnerable than other regions. In Supplementary Fig. 4, the yearly resolution results? can be found.
b, Overall tipping reason averaged over the entire Amazon basin with error bars as the standard deviation over all years and all 100 ensemble members. A version separated into the future drought conditions from 2004 to 2014 can be found in Supplementary Fig. 5 for all these potential future drought scenarios. MAP does not contribute to tipping events (probability is less than 0.1%) and is thus omitted from this gure. c, Tipping reason map: MCWD, d, Tipping reason map: Network (Cascading effects of the moisture recycling network). Note that panel a is the sum of the panels c and d.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. Supplement.pdf