Multiple drivers and extreme events collapse social-ecological systems sooner


 The world’s ecosystems are undergoing unprecedented changes due to the impact of climate change and local human activities. A major concern is the possibility of tipping points where ecosystems and landscapes change abruptly to undesirable states. We consider what happens to the timing of tipping points when current stresses strengthen whilst systems experience additional stresses and/or extreme events. We run experiments on four mathematical models that simulate tipping points in lake water quality, the Easter Island community, the Chilika lagoon fishery, and forest dieback. We show that the strongest impacts occur under increasing levels of primary stress, but additional and more extreme stresses in all four models bring the tipping points significantly closer to today. Translating the results to the real world underlines the need for humanity to reduce damaging disturbances and global warming, and to be vigilant for signs that natural systems are degrading more rapidly than previously thought.


31
For many observers, UK Chief Scientist's John Beddington's argument that the world faced a 'Perfect 32 Storm' of global events by 2030 1 has now become a prescient warning. Recent mention of 'ghastly 33 futures' 2 , 'widespread ecosystem collapse' 3 , and 'domino effects on sustainability goals' 4 tap into a 34 growing scientific consensus that the Earth is rapidly destabilising through connected tipping points driver 20 . Taking a network approach, Krönke et al. found that systems with large average clustering coefficients and spatially structured networks (e.g. the Amazon forest) are more vulnerable to tipping 72 cascades than more disordered network topologies 21 . Such arguments were used to explain Cooper on ecosystems 28,29 , we are unaware of published modelling studies that simulate the interaction of 104 fast drivers, multiple drivers and system noise on the timing of tipping points in real world systems.

106
The evidence for increasing rates of change in key drivers at the global scale is strong. Despite 107 decreases in global birth rates and increases in renewable energy generation, the general trends of 108 population, greenhouse gas concentrations and economic drivers, such as gross domestic product, is 109 upwards 8,30 , often with acceleration through the 20th century. Similar non-stationary trends for 110 ecosystem degradation 31 imply that unstable sub-systems are common. Furthermore, there is strong 111 evidence globally for the increased frequency and magnitude of heatwaves and heavy precipitation 112 events 10 . Examples include the sequence of European summer droughts since 2015 32 , and fire-113 promoting phases of the tropical Pacific and Indian ocean variability 33 . It also seems certain that the 114 number of extreme events that cause fatalities or normalised financial costs above a range of 115 threshold values has increased since 1980 8 . Human influence has also likely increased globally the 116 frequency of concurrent heatwaves and droughts, already implicated in reduced crop yields 34 and 117 regional fires and flooding 10 . Overall, global warming even at 1.5°C will increase the frequency of 118 unprecedented extreme events, raise the probability of compound events (e.g. 35 ), make multiple 119 system failures more likely 36 , and make rarer events more common 10 .

121
Thought experiments informed by these observations suggest that stronger interactions between 122 systems may be expected to increase the numbers of drivers of any one system, change driver 123 behaviour and generate more system noise. For example, for any particular element (e.g. the Amazon 124 forest) it is possible to envisage a time sequence that starts with one main driver (e.g. deforestation), 125 then multiple drivers (e.g. deforestation plus global warming), more noise through extreme events 126 (e.g. more droughts and wildfires), with additional feedback mechanisms that enhance the drivers multiple drivers and connectivity 25 . As an alternative we turn attention to systems models that can be 139 externally manipulated to simulate internal emergent tipping dynamics at local and regional scales - phosphorus: lake phosphorus concentration; TRIFFID: tree coverage). As expected, when the strength 167 of the primary slow driver in each model is increased, the modelled systems collapse sooner (as 168 defined by a statistical breakpoint in their temporal trends; see Methods Section 3; Figure 1). We show 169 that in three of the models (Easter Island, TRIFFID and Lake Phosphorus), changes in primary driver 170 strengths have comparatively larger reductions in the absolute time before the tipping point is 171 reached at low driver levels (i.e., when the primary driver is relatively weak), when compared to the 172 same change at higher values of the primary driver, although the percentage value by which the 173 tipping point is brought forward is relatively consistent (Figure 2). In other words, the absolute 174 advance in the tipping point may be reduced as driver strength increases, but the percentage shift in  For each of the scenarios (i.e. driver combinations), the boxplots depict the 2.5%, 25%, 50%, 75% and 97.5% breakpoint date percentiles, in each of the following ranges along the x-axis: 0.25-0.35, 0.45-0.55, 0.65-0.75. Subplots: (A) Chilika model, primary slow driver: fisher population growth, secondary driver: climate change, tertiary driver: fish price; (B) Easter Island model, primary slow driver = tree clearance, secondary driver: agricultural carrying capacity, tertiary driver: tree mortality; (C) TRIFFID model, primary slow driver: temperature change, secondary driver: disturbance rate; (D) Lake Phosphorus model, primary slow driver: phosphorus external input, secondary driver: phosphorus recycling rate, tertiary driver: phosphorus sedimentation rate. The dashed vertical lines denote the minimum primary driver strength associated with a breakpoint.
Increasing the strength of multiple (i.e., secondary and tertiary) drivers further reduces the breakpoint 183 date (Figure 2: coloured points). As with a single driver, this effect is more noticeable on unit time at 184 low driver levels, but the percentage effect is relatively consistent. For example, at a normalised driver    Figure 2D). With all 208 additional drivers, 12.3% of tipping points observed in the Lake Phosphorus model occurred at primary 209 driver strengths below the minimum threshold required to result in a tipping point when the primary 210 driver is acting in isolation (Lake Chilika: 1.2%; Easter Island: 14.8%; TRIFFID: 7.7%; Table S1).

212
Next, for each of the four models, the trajectories of the primary slow drivers were randomly 213 perturbed by the addition of noise (Methods Section 2.3). Noise was generated within the system 214 dynamics software used to run the models (STELLA 51 ) by randomly sampling per timestep from a  However, the addition of high noise (normalised σ values > 0.666) highlights that increasing the Island model ( Figure 3B), the addition of high noise brings the median tipping point forward from 228 timestep 1179 to timestep 848 (69.6% reduction); the same levels of noise at a 0.7 normalised baseline 229 driver strength advance the date of collapse from 663 to 290 (56.4% reduction). In turn, the equivalent 230 advances in the breakpoints of the other three models under high noise are 11 (34.4%) and 6 (24.0%) 231 years for Lake Chilika, 47 (19.7%) and 18 (16.5%) years for TRIFFID, and 104 (12.3%) and 61 (14.5%) 232 timesteps for Lake Phosphorus, for 0.3 and 0.7 normalised baseline driver strengths respectively.

233
Combined, addition of low, mid and high noise levels resulted in 2.5%, 6.6%, 2.0% and 2.6% of 234 modelled tipping points occurring at primary driver strengths below the minimum threshold required 235 to result in a tipping point when acting in isolation for Lake Chilika, Easter Island, TRIFFID and Lake

236
Phosphorus respectively (Table S2). These results are consistent regardless of whether the noise is coupled to the magnitude of the primary slow driver or not ( Figure S2; Table S3).

252
The effects outlined above are additive, and so combining multiple drivers with noise further reduces

263
(54.2% reduction) for Lake Phosphorus. Across all combinations of noise and multiple drivers, 1.7%, 264 7.5%, 6.6% and 8.9% of modelled tipping points occurred at primary driver strengths below the 265 minimum threshold required to result in a tipping point when acting in isolation for Lake Chilika, Easter 266 Island, TRIFFID and Lake Phosphorus respectively (Table S4). These results are consistent when noise 267 is coupled to the strength of the primary slow driver ( Figure S3; Table S5).

287
Previous findings have strongly supported the idea that Earth subsystems may interact to the extent 288 that an abrupt shift in one raises the probability that a shift may occur in another 20,21,52 . In this paper 289 we have explored through four social-ecological models how these interactions may bring forward the timing of tipping points through the effects of strengthened drivers, multiple drivers and higher 291 internal variability or noise. The potential effects are significant with combinations of strengthened 292 main driver, an additional driver and noise giving at least 38-81% reductions in the future date of a 293 predicted tipping point compared to estimates for a non-interacting system with a constant single 294 driver and no noise. Importantly, the effect per unit time on bringing forward a tipping point is 295 greatest at low driver trajectories, which further compounds the suggestion that abrupt Earth system 296 changes may occur sooner than we think. Our findings also show that 1.2-14.8% of tipping points can 297 be triggered by additional drivers and/or noise below the threshold of driver strengths required to 298 collapse the system if only a single driver were in effect.

300
Overall, we find that as the strength of a main driver increases, the systems collapse sooner. Adding 301 multiple drivers further brings forward collapses, as does adding noise, and the two effects can be 302 additive. However, the relative importance of these changes varies across systems. For the Chilika 303 fishery, the most influential driver is captured as the primary driver and so additional drivers have 304 limited influence. However, the addition of noise in the primary driver brings the breakpoint date 305 much further towards the present. For Easter Island, Lake phosphorus, and TRIFFID, the opposite is 306 true -the addition of high levels of noise in the primary driver advances the date of system collapse 307 far less than additional drivers. Thus, while the earliest collapses in all the systems are found when 308 both additional drivers and noise are applied, the precise importance of individual drivers and noise 309 may be system-dependent.

311
These results have research implications for developing and applying tipping point models. Whilst our 312 findings derive from models based on real-world systems, the superior complexity of reality may limit 313 the transferability of our results. For instance, it is plausible that more complex systems will have 314 stronger regulating mechanisms that stabilise the system through sets of balancing feedback loops 53 , 315 constraining the more extreme of our findings. Model development is therefore required to better 316 capture the diversity of system elements, interactions and feedbacks for more complex systems.
317 Furthermore, each model in this study was originally created to study the impact of a primary driver, 318 presumably perceived as the most impactful. Our results show that this may not be the case (e.g.

319
Easter Island) and models should include a range of plausible drivers if they are to avoid  or common driver to many systems, and which may be expected to strengthen before it subsides.

353
Clearly, climate economics need to incorporate these additive and cumulative effects that are 354 occurring at local and regional scales into larger scale models where they are currently lacking 62,63 . The worrying. Similarly, the implication for regions experiencing more extreme events is that a tipping 358 points may occur even before the main driver has ramped up.

360
The commonality of findings across four well-studied social-ecological systems has potentially 361 profound implications for our perception of future risks associated with the climate and ecological

387
The Lake phosphorus model is a simplified version of the original 'lake response to P input and 388 recycling' model 65 , as modified by Wang et al. 48 . Lake water phosphorus concentration is driven by 389 external phosphorus input. In turn, lake water phosphorus is recycled back into the system as a Phosphorus is also removed from lake waters via sedimentation, where the volume removed in sediment is proportional to the phosphorus concentration of the lake. Therefore, on any given

416
is driven by a temperature dependent growth term and an externally set disturbance rate: Where v is the vegetation coverage, T f is the temperature forcing parameter (Methods Section 2.3), g 421 is the vegetation growth rate, g 0 is the maximum growth rate (2/year), y is the disturbance rate surface bare soil and forest (5°C). Therefore, the growth term is assumed to be parabolic with the local 425 temperature (Equation 2b), meaning that once the local temperature increases beyond the optimal 426 temperature, negative tree growth ensues (i.e. additional tree mortality) 37

467
Baseline outputs were generated with the Primary driver active AND the Secondary and Tertiary driver 468 remaining at its default value AND the Noise level remaining at zero (Table S7). In turn, the Monte

469
Carlo sensitivity analysis function in STELLA randomly samples a future change trajectory for the 470 primary slow driver per simulation (as plotted on the horizontal axes of Figures 2-4). The primary 471 trajectory is sampled between the lower and upper limits of uniform distribution bounds, meaning 472 that there is a uniform likelihood of selecting any given trajectory between the bounds (Table S7) Table S7.

476
The built-in STELLA 'TIME' function generates future scenario trajectories that change linearly over     (Table S7); in turn, the 512 TIME function linearly increases the value of the driver from its default value to its sampled trajectory 513 value by the final timestep of the model horizon.

515
In order to produce one secondary trajectory per simulation in the Chilika model, the sampling of 516 rainfall and temperature trajectories are connected along a linear gradient, ranging from no change 517 to a combination of +30% rainfall change and +4.5°C temperature change by 2081-2100 relative to 518 1986-2005 (as per RCP8.5 projections for the east coast of India 66 ). In STELLA, this is operationalised 519 by the model variable 'climate change trend', with Monte Carlo sensitivity analysis in STELLA randomly 520 sampling a value between 0 and 1 per simulation. As an illustration, if a value of 0.6 was to be sampled, In order to investigate Experiment #3 and Experiment #4, the value of each primary slow driver is 526 perturbed per timestep by randomly generated noise. We simulate a standard Wiener process to 527 generate noise, equal to √ × (0,1), where 'dt' equals change in time and 'N(0,1)' is a normal 528 distribution with a mean of 0 and standard deviation of one. In turn, the product of the Wiener process 529 is multiplied by the scaling factor 'sigma' (σ), providing an overall level of noise to be added to the 530 value of the primary driver on any given timestep. As per the future trajectories, the magnitude of 'σ' 531 is randomly sampled once per simulation from uniform distributions, with lower and upper limits 532 shown in Table S7.          (Table S1). The gradient trajectories representing 30% and 70% of the maximum possible trajectories have been identified in orange and red, respectively. The same approach is also used for all secondary and tertiary drivers, evolving along randomly selected linear gradients between 'no change from default' and 'maximum change' by the model horizon. Experiment #1 [1] : change in the trajectory of the primary driver only; Experiment #2 [2] : changes in the trajectories of the primary driver and additional driver trajectories; Experiment 3 [3] : changes in the trajectories of the primary driver and the addition of noise to the primary driver; Experiment #4 [4] :

SI2 -Supplementary
change in the trajectories of the primary driver and additional driver trajectories, plus the addition of noise to the primary driver. Experiments #3 and #4 here include both coupled and uncoupled noise.