Abrupt Climate Change in an Oscillating World

The notion that small changes can have large consequences in the climate or ecosystems has become popular as the concept of tipping points. Typically, tipping points are thought to arise from a loss of stability of an equilibrium when external conditions are slowly varied. However, this appealingly simple view puts us on the wrong foot for understanding a range of abrupt transitions in the climate or ecosystems because complex environmental systems are never in equilibrium. In particular, they are forced by diurnal variations, the seasons, Milankovitch cycles and internal climate oscillations. Here we show how abrupt and sometimes even irreversible change may be evoked by even small shifts in the amplitude or time scale of such environmental oscillations. By using model simulations and reconciling evidence from previous studies we illustrate how these phenomena can be relevant for ecosystems and elements of the climate system including terrestrial ecosystems, Arctic sea ice and monsoons. Although the systems we address are very different and span a broad range of time scales, the phenomena can be understood in a common framework that can help clarify and unify the interpretation of abrupt shifts in the Earth system.


Supplementary Information
Overdamped Duffing oscillator As a conceptual model for illustratory purposes we use the overdamped Duffing oscillator given by Eq. 1 and Eq. 2 in the main text and repeated here: To perform the two cases of a saved and collapsing system ( Fig. 1) we first run the system into equilibrium for constant D = -0.1. Thereafter, we let D make one single cycle with the following parameters:   In all above cases, we solve the system numerically with a Runge-Kutta scheme of fourth order, using 50.000 time steps per cycle.

Monsoon model
The monsoon model by Levermann et al. 8 describes the energy and water balance of the atmospheric column over land. The monsoon flow from the ocean advects moisture over land which then condensates in the rising air. The condensational heat is radiated to space.
− ∆ + = 0 (Eq. S1) Eq. S1 describes the balance between condensational heat, temperature advection and radiation in the atmospheric column. Eq. S2 describes the steady state water balance (moisture convergence = condensation). Eq. S3 describes the momentum equation (temperature gradient drives circulation). Eq. S4 sets a minimum moisture threshold below which no precipitation is possible (monsoon failure). The four equations can be combined to yield one algebraic equation that is cubic in P and therefore has three solutions. Two of these solutions however, are discarded as unphysical due to a wrong direction in the circulation (W<0) or negative precipitation (P<0) 8,9 . The remaining solution describes the "on" state of the monsoon. In a parameter regime where only imaginary or unphysical solutions exist, the monsoon is thought to be in its "off" state, and precipitation is then set to 0.
The model is only valid for summer conditions and cannot describe any precipitation not associated with the summer monsoon. We therefore set P to zero in other seasons, knowing that no monsoon rainfall occurs during that time. Moreover, the model ignores the sensible heat flux at the surface, SH, because it is small in summer. Here we do not make this assumption but lump together SH and R to become R* = R+SH. The equations then do not change, except that R is replaced by R* with a different meaning. R* and the specific moisture over the ocean, qo, are the two external parameters that we assume to oscillate with a period of one year.
In order to obtain an asymmetric annual cycle with an abrupt onset like in observations, we apply two strategies: 1. We follow the approach by Schewe et al. 9 by introducing a small bistable regime at the critical threshold (the shaded region with a width of 0.5 g/kg at the boundary).
2. We introduce a phase lag of one month between the two cyclic drivers, making oceanic moisture qo lag the seasonal cycle in R*. This approach may be justified by the different memories of ocean and atmosphere: while the atmosphere can adjust to the solar insolation within days, the oceanic temperature which controls specific humidity lags the insolation cycle by many weeks.

Arctic sea ice model
The Arctic sea ice model by Eisenman and Wettlaufer 5 in the version of Eisenman 68 calculates an energy balance for a well-mixed box of ocean water, covered with ice of a single thickness. The incoming shortwave and long-wave radiation at the surface are prescribed as harmonic oscillations. The only dynamic state-variable in the model is enthalpy E. When ice is present, the enthalpy is negative and proportional to the ice thickness hi. In the absence of ice, E is proportional to the ocean temperature To: Li is the latent heat of melting/freezing, co the water's heat capacity, Ho the mixed layer depth, and Tm the mixed layer temperature.
The evolution equation of the model is a simple energy budget: where B*T is the linearized blackbody emission and A is the sum of all temperature-independent fluxes: = ( + 2 tanh ( ℎ )) ( − cos 2 ) − [ + cos 2 ( − )] (Eq. S7) The bifurcation parameter Lm represents the annual mean outgoing long-wave radiation budget at the surface (not including temperature feedbacks) and becomes smaller when the climate warms. For values of Lm between approx. 65 W/m 2 (which roughly represents the present-day climate) and 50 W/m 2 , two stable solutions can be found: a seasonally ice-covered ocean and an ice-free ocean (Fig. 5a).
The temperature T in Eq. S6 represents the difference between the surface temperature and the freezing temperature of ice. In the presence of sea ice, the temperature distribution within the ice layer is assumed to be linear. During melting, T is at the melting point: In our experiments without annual cycle, we set time t in Eq. S7 to a constant day of the year, but still integrate Eq. S6 over time. Under permanent winter conditions, sea ice then grows to an infinite thickness. In such cases we limit E to an arbitrary lower bound of -7 GJ/m 2 (corresponding to a thickness of approx. 23 m).

Vegetation model
The model by Zeng et al. 87 describes the interaction of vegetation V and precipitation P in North Africa. V can be interpreted as the fraction of leaves relative to their maximum coverage. The model is a nondimensionalised version of earlier, very similar concept models 6 . The only dynamic equation describes the evolution of V, based on a relaxation approach: As rainfall P adjusts much faster to changes in V than vice versa, P is determined by an algebraic equation.
The sinusoidal signal captures the low-frequency fluctuations that arise from sea surface temperatures in an idealised way.    S3. Dependence of the example system on initial conditions for different parameters Dm (horizontal axis), amplitude Da (vertical axis) and period T (subfigures). The color shows the time mean of the absolute difference between the system's solutions. This difference is 0 (blue) if there is only one solution (monostable regime) and > 0 if the system is bistable; the boundaries between these regimes hence consist in bifurcation points. The model then calculates a vegetation coverage whose time mean is shown on the vertical axis. As in Fig. S1-S3, two stable states exist for small amplitude forcing, and the transition to the remaining state is more gradual for fast forcing (low T) than slow forcing (large T). ©American Meteorological Society. Used with permission.