Fragmented tipping in a spatially heterogeneous world

Many climate subsystems are thought to be susceptible to tipping—and some might be close to a tipping point. The general belief and intuition, based on simple conceptual models of tipping elements, is that tipping leads to reorganization of the full (sub)system. Here, we explore tipping in conceptual, but spatially extended and spatially heterogenous models. These are extensions of conceptual models taken from all sorts of climate system components on multiple spatial scales. By analysis of the bifurcation structure of such systems, special stable equilibrium states are revealed: coexistence states with part of the spatial domain in one state, and part in another, with a spatial interface between these regions. These coexistence states critically depend on the size and the spatial heterogeneity of the (sub)system. In particular, in these systems the crossing of a tipping point not necessarily leads to a full reorganization of the system. Instead, it might lead to a reorganization of only part of the spatial domain, limiting the impact of these events on the system’s functioning.

Typically, tipping is illustrated and explained using simple, conceptual low-dimensional models, that have two alternative states and that can tip between them as climatic conditions change [1][2][3]20].In more complex, more detailed highdimensional models and in real-life data tipping is, however, often not as clear and pronounced [4,5,21].Tipping from one state to a completely differently structured state is hardly ever observed.Instead, partial restructurings occur more often.For instance, (large) parts of an ice sheet melt, instead of the whole sheet melting in one single tipping event [22].This suggests that low-and high-dimensional models behave differently.It could be that highdimensional models are tuned for stability too much, suppressing tipping behaviour [23].It could also be that the low-dimensional models are too restrictive in the number of physical processes, thereby exaggerating tipping behaviour [24,25].At least, the most simple models really only allow for two alternative states and nothing more.Adding complexity to these leads to more response options for the system, which might lead to less severe tipping events.For instance, adding more boxes to a box model [26,27], or incorporating spatial effects [21,25].
In this paper, we investigate the behaviour of conceptual models when spatial effects are incorporated: spatial transport and spatial heterogeneity.This setting has received only little attention in the literature [21,28] and a thorough theoretical understanding of such systems is still lacking, despite their omnipresence [29].In such models additional stable states called co-existence states can emerge, in which part of the domain resides in one state and the rest in another state, with a spatial interface separating these regions [21,24]-see figure 1 for real-life examples.Consequently, in these systems transitions can occur in which only in part of the spatial domain the system changes state, providing a more subtle, fragmented tipping pathway.
The rest of this paper is structured as follows.In section 2, we first review the classic theory of coexistence states in spatially homogeneous (gradient) systems.Subsequently, we detail how this theory changes in a spatially heterogeneous setting.We focus on the possibility of new equilibrium coexistence states and the different bifurcation diagrams these systems can have depending on the spatial heterogeneity.Then, in section 3, we illustrate the potential widespread relevance of coexistence states using several examples of climate subsystems on different spatial scales.For this, we use a variety of conceptual models, that have been proposed before in the literature, or spatially extended versions thereof.Finally, we end with a discussion in section 4.

Theory
We consider the evolution of a (single) state variable y.Depending on the particular system of interest, this could be for example temperature or vegetation coverage.Ignoring spatial effects for a moment, the local dynamics of y are described by an ordinary differential equation, i.e.
where f describes the evolution over time, depending on the value of the state variable y and the value of a (bifurcation) parameter µ, and V is the associated potential function (also referred to as the stability landscape) [1][2][3].The model (1) possesses tipping behaviour when there is a parameter region in which there is bistability of two different states-say, a state A and state Bwith tipping points (saddle-node bifurcations) to the other state at both ends of this region.This leads to the prototypical 'S'-curve bifurcation diagram as shown in figure 2(a).An explicit example of such system is the model (1) with In the rest of this section, we detail how the addition of spatial effects to the model (1) changes the bifurcation diagram.

Spatial transport
We now consider the evolution of the state variable y on a spatial domain with coordinate x.In addition to the local dynamics, spatial transport between locations now also plays a role.The most simple way to implement spatial transport is to add (linear) diffusion to the model (1).Thus, a new spatial model describing the evolution of y over time and space is given by the partial differential equation where D is the diffusion coefficient and ∆ the spatial Laplacian in x.
Because this type of model can be used to model phase segregation dynamics [31,32], its dynamics are well understood and the literature is quite extensive, including rigorous mathematical proofs [33,34].First, spatially uniform states exist in which the whole domain is in the same state (either state A or B).These are the natural spatially extended versions of the states in the local model (1).However, next to these, there are also coexistence states, in which parts of space reside in state A and the rest in state B, with a spatial interface (or 'front' in mathematical jargon) connecting these regions.These interfaces can be spatially localized, or extend over a large space depending on the strength of the spatial transport.
The dynamics of these coexistence states is intricate [33,35].In short: the interfaces typically migrate slowly over space depending on the distance to other interfaces and the domain boundary, but there are sometimes fast coalescence events when interfaces meet and separate regions merge.On the long term, a maximum of one interface can persist in the domain (the rest is slow transient behaviour).Such single interface moves with a speed that is proportional to the difference in potential between states A and B, i.e.V (A; µ) − V (B; µ).So, only in the degenerate case in which the potentials are equal, this leads to an additional equilibrium state of (3).The specific parameter value for which the potentials are equal is called the Maxwell point of the system.This means that the bifurcation diagram does not change much compared to the non-spatial system; see figure 2(b).However, the transient dynamics occurring after a tipping point has been crossed can be different in such systems: this can now occur via an invasion front which replaces one state with the other, which might lead to a slow tipping of the full system [36,37].

Spatial heterogeneity
The previous model (3) assumes that the local dynamics are the same throughout the whole domain, i.e. the model is spatially homogeneous.Reality is, however, spatially heterogeneous: local dynamics differ from location to location, for example due to topographical features or human activity.Hence, it is more realistic to consider a spatially extended model in which the local dynamics explicitly depends on the location x: We will assume that the spatial variation is not too much locally, i.e. ∂f ∂x ( y, x; µ) does not get very large compared to diffusion strength (as this would, potentially, give rise to different dynamics than described in this paper; see e.g.[38,39] for an exploration of those dynamics in case of spatially periodic forcings).
The presence of spatial heterogeneity in (4) does change the dynamics of the interfaces.Their movement now depends on the local difference between the potentials of state A and B. Because of this, the Maxwell point is different for different locations in space, making it no longer a degenerate case.Furthermore, coexistence states with multiple interfaces can now be equilibrium states of the model system.
A typical bifurcation diagram is given in figure 2(c), which is different compared to the classic one (figure 2(a)) because of the additional branch of stable coexistence states.In this case, the crossing of a tipping point does not necessarily lead to a tipping of the complete system, but can lead to a less critical shift of the system in which only part of the spatial domain undergoes a transition to the alternative state.
The precise structure of the bifurcation diagram depends on the specific heterogenity of the system.A few possible variations are shown in figure 3, illustrating the potential complexity of the bifurcation diagram for spatially heterogeneous systems.The results show that the parameter range for which coexistence states exist can vary much between systems with different heterogeneities.In particular, this range can be fully contained in the hysteresis loop of the noncoexistence states (e.g.(a) and (d)), which indicates that crossings of the tipping points will lead to a full tipping of the system despite the presence of coexistence states in part of the parameter regime.On the other hand, the spatial heterogeneity can also prevent overlap of the fully tipped branches (e.g.(b)), meaning that the full system cannot tip in one tipping event.It can also happen that the tipping process is fragmented in one direction, but not in other directions (e).Moreover, the spatial heterogeneity can lead to multiple tipping points (e.g.(c), (d) and (f)), in each of which another part of the spatial domain undergoes a transition to the alternative state.
A concise description of the mathematical theory and analysis of coexistence states can be found in appendix A.
The above described spatial interfaces between regions in different states can only emerge and persist if there is enough space available for them in the system.Hence, if the spatial domain is too confined, the above described coexistence states do not fit and the branch of coexistence states is absent in the bifurcation diagram (see figure 4).The precise minimum required length varies per model and depends on the typical length scale associated with the spatial transport.For instance, if the diffusion coefficient D in ( 4) is increased, a larger domain is needed to facilitate the coexistence states.See also appendix B for a mathematical treatment of the small domain limit.
The form of equation ( 4) is chosen for simplicity of presentation.Coexistence states can also occur in other models that for example have a diffusion coefficient that depends on space or have a different spatial transport mechanism altogether.Similar behaviour as described above can at least be expected from other spatially heterogeneous Lagrangian or gradient systems.

Earth system and ecosystem example conceptual models
The above described theory is quite generic, and there are only few assumptions.Hence, if a simple conceptual model only considers local dynamics and captures tipping behaviour, a spatially extended and heterogeneous model will behave as described above.In particular, such systems might have additional coexistence states, and tipping events in which only part of the domain undergoes transition.Many systems, and many previously proposed conceptual models, fit into this description.To illustrate this, next we give some examples of systems for which this theory might be relevant by looking at some illustrative conceptual spatially extended models for these systems.These are either models directly taken from the literature or simple spatial extensions of non-spatial models.Hereby we cover systems of many different spatial scales, ranging from global to more regional systems.

Global energy balance model
On a planetary scale, an example of coexistence states can be found in global energy balance models.These describe the evolution of Earth's temperature by considering the energy flux at the top of the atmosphere.These type of models were first introduced by Budyko [40] and Sellers [41], and have since served often as conceptual models to illustrate planetary climate shifts [42][43][44][45][46][47][48][49][50][51][52].
Here, we consider the following simple spatially one-dimensional variant: This models the longitudinally averaged temperature T (as function of x = sin(θ) where θ is the latitude) by looking at the absorbed incoming solar radiation (Q 0 (x) [1 − α(T)]), the outgoing Planck radiation (−εσ 0 T 4 ), the meridional heat transport , and the effect of atmospheric CO 2 on the energy budget (+µ).The model uses the following functional form for the temperaturedependent albedo (ice, with a high albedo, can only exist if temperatures are low enough): where α 1 is the albedo of ice, and α 2 the albedo of water.The sigmoid function ensures the change in albedo is smooth.The spatial heterogeneity in this models stems from the incoming solar radiation,which is latitude-dependent.Following [48], we have taken where Q 0 is the global average of the incoming solar radiation.As parameters we have taken The value for the diffusivity D is not known-and potentially not uniform over earth [44,45]-but for the illustrative purposes in this paper we have taken D = 0.30.Finally, since we are focussing on stationary states, the value for C T does not matter, but a typical value would be It has been shown that the model can reside in an ice-free state or a fully ice-covered 'Snowball Earth' state.When spatial effects are not modelled, these are the only states of the model.However, once spatial effects are incorporated, planet-scale coexistence states are also found with ice in only part of the earth (i.e. the poles) and no ice elsewhere [42][43][44][45][46][47][48][49][50][51].An example bifurcation diagram is given in figure 5(a) which shows the global mean temperature as function of the parameter µ.This diagram shows that tipping between a snowball earth state (blue branch) and a partially ice covered state (cyan branch) leads to large changes in global temperature, whereas tipping between coexistence states (cyan branch) and ice-free states (red branch) leads to much smaller changes in global temperature.It is also interesting to note that there is no bistability of snowball earth and ice-free states, implying that transitions between them always need to go via coexistence states.

Ocean circulation
In oceans, convection occurs at places where the ocean mixed-layer exchanges water with the deep ocean.This ocean convection plays an important role in the global ocean circulation.Using simple conceptual models, the possibility of bistability between a convective and a non-convective state has been found for a range of parameters [53][54][55][56].These results have been obtained in box models, where temperature T and salinity S of the mixed-layer are dynamic variables, which change through exchange with (static) atmosphere and deep ocean boxes.Here, we extend such model via the addition of diffusion as form of spatial transport.
To showcase coexistence states with convection only at certain regions, we apply a series of (standard) simplifications on the equations (details in appendix C).This leads to the following equation that describes the evolution of the local density difference ∆ρ(x, t) between the mixed-layer and the deep ocean: Here, ∆ρ A (x) represents the (static) density difference between the atmosphere box and the deep ocean box, and ρ 0 (x) is the local density in the deep ocean box.Further, k T gives the rate of exchange with the atmosphere box and κ(∆ρ) the rate of exchange with the deep ocean box.The latter is a function of ∆ρ, which is a non-negative increasing function [53].
That is, convection happens when the exchange rate κ(∆ρ) is large, which typically occurs if the mixedlayer is heavy enough compared to the deep ocean.Depending on the functional form of the exchange rate κ(∆ρ), certain parameter ranges can allow for bistability and coexistence states between convection and no convection.Not striving for the most realistic description, we take a continuous approximation of a step function similar to [53].Specifically, we take κ(∆ρ) = To illustrate this model, we take a spatially heterogeneous atmospheric reference density ∆ρ , where f acts as a bifurcation parameter.In ocean systems, such kind of spatial heterogeneity could for instance be caused by local differences in freshwater fluxes.In figure 5(b), an example bifurcation diagram for this model is given, showing the possibility of coexistence states in large enough ocean basins.In this system, the system does not tip fully between convective and non-convective states at a tipping point, but only part of the domain changes and the system transitions to a coexistence statefor chosen heterogeneity, in these coexistence states convection happens only at the edges of the spatial domain.

Atmospheric circulation
In the atmospheric circulation, there also seems to be the possibility of coexistence states.As a simple example, we look at the dynamics of the atmospheric boundary layer, in which heat exchange with the surface takes place.This boundary layer can be in a fully turbulent state or a quiescent, quasi-laminar state depending on the temperature difference between the boundary layer and the soil.In cold regions, or during cold nights, the layer is stable stratified, but at higher temperatures surface warming leads to a convective boundary layer.
Bistability and tipping between these two states has been explained using simple energy balance models [57,58].These describe the near-surface inversion strength ∆T, defined roughly as the boundary layer's temperature minus the soil temperature, which evolves according to the net imbalance of the energy fluxes in the layer [58].To the best of our knowledge, there is no standard way to include spatial heat redistribution.Hence, as a first effort to create a spatially extended model, we have modelled spatial heat re-distribution as a diffusive process.The model is given by Here, Q(x) − λ∆T is the linearized net long-wave radiation, and −C(x)∆Te −2α(∆T) is the turbulent sensible heat flux.The location x can be interpreted as a local coordinate of a larger region on Earth.In the numerical computations, we have rescaled the spatial domain to [−1, 1] and have taken a rescaled diffusion coefficient D = 0.01.In reality, this would correspond to a spatial domain that is much larger than the diffusion length scale.
To demonstrate the effect of spatial heterogeneity in this model, we take λ = 1, α = 3, Q(x) = 1 − 0.1x and C(x) = 0.1x + µ, modelling for example spatial differences in soil temperature or wind speeds (the latter influences the turbulent sensible heat flux); µ acts here as a bifurcation parameter, representing changes in wind speed.For these choices, an example bifurcation diagram is given in figure 5(c), which indicates the possibility of coexistence states with regional turbulence.However, for these heterogeneous conditions, the branch of coexistence states is almost fully contained in the hysteresis loop of the fully turbulent and the fully quasi-laminar states, suggesting that crossings of the tipping points might lead to a full tipping of the system despite the presence of a branch of coexistence states in the bifurcation diagram.

Ecosystems
In many ecological systems, bistability and transitions between alternative states have been found.For instance, in drylands, bistability between vegetation and bare soil occurs [59].On larger scales, the same environmental conditions can even support multiple ecosystem types, such as grasslands, savannas and tropical forests in the Tropics [7,60].In all of these systems, spatial effects are present: species disperse over space (e.g. via seed dispersal) and environmental conditions are spatially heterogeneous (e.g. a rainfall gradient).Hence, they can exhibit coexistence states in which different ecosystem types are separated by a local interface.
Especially in models of drylands, many patterned states have been analyzed, including coexistence states [37,61], although spatial heterogeneities are not often taken into account, with few exceptions [62,63].However, the models considered are typically more complicated than (4) and more advanced mathematical techniques are used to analyze them that are beyond the scope of this article.Hence, we illustrate the coexistence patterns in ecosystems using a different example.
Recently, the bistability between tropical forests and grasslands in a heterogeneous environment has been studied [64].In a spatially extended (rescaled) version of the model in [65], the local fraction of fireprone forest trees is given by F and its evolution is modelled as with the coordinate x representing e.g. a crosssection through a tropical forest region.Here, r(P) = 0.2 ( 1 − e 1.54−0.003P ) is the (logistic) growth rate depending on local precipitation P(x).Further, m(P) = 0.041 + e −2.15−0.008P is the precipitationdependent natural mortality rate, and f (F; P) is the mortality due to forest fires, which depends both on the local precipitation P(x) and the local tree coverage F (dry trees burn easier and fire spreads easier when there are less trees and more herbaceous vegetation).Specifically, f takes on the following functional form where Y C (P) models the assumed decrease in firespreading percolation threshold (see [65] for details), expressed in equations as Finally, sapling dispersal is modelled as diffusion.
The spatially heterogeneous term in this example is the local precipitation P(x).Following [64], we have modelled this as a linear precipitation gradient P(x) = P mean + 150x, and model climate change via the mean precipitation parameter P mean .An example bifurcation diagram is given in figure 5(d).Here, the red branch has solutions with non-zero forest tree fraction F everywhere, the blue branch corresponds to a no-forest state and the cyan branch to coexistence states with forest trees in part of the domain.In this diagram, the blue branch is attracting for all parameter values, and does not connect to the other branches via a saddle-node bifurcation (the lower unstable branch only connects to the blue branch in the limit P mean → ∞).From the bifurcation diagram it can be deduced that the spatial region with lowest rainfall (left in the insets) will die down when the fully tree-covered state tips into a coexistence state.However, as long as the system is in a coexistence state and there still are trees in part of the domain, restoration is possible if precipitation increases again.

Shallow lakes
On even smaller spatial scales, coexistence states could also play a role.As an example, we look at the turbidity of shallow lakes [15], which has been used as one of the prime examples of critical shifts in ecological systems.In these lakes, turbidity is closely related to the algae concentration A, which can change depending on the amount of nutrients in the lake.If there is a lot of nutrients, algae concentration is high and the lake is turbid.For lower nutrient concentrations, there are few algae and the lake is clear.Transitions between these states can take place if the nutrient concentration changes.
In [15], a simple non-spatial model is given for the evolution of algae A, depending on the nutrient concentration.Here, we have added spatial transport as a diffusive process to obtain the following spatially explicit model that models the local algae fraction over e.g. a cross-section of a shallow lake: Here, algae growth is modelled as logistic growth with competition coefficient c.The growth rate depends on the effect of nutrients μ(N, x) and the amount of vegetation V(A), which is larger if there are few algae.The value of p(x) represents the (local) shallowness of the lake.More information on the rescaling process can be found in appendix D.
To illustrate the effect of spatial heterogeneity in such lake systems, we take μ(x) = µ (1 − 0.1 cos(πx)) and p(x) = 3 + 0.1 cos(πx), representing e.g. higher nutrient concentrations and shallower waters near the edge of a lake.The other parameter values are r = 10, h v = 0.1, c = 1 and D = 0.01 (with scaled spatial domain [−1, 1], essentially modelling a lake large enough to allow for coexistence states).As bifurcation parameter we take the spatial average nutrient effect µ.An example bifurcation diagram is given in figure 5(e).This shows the possibility of coexistence states in such systems.In particular, for this specific heterogeneity, the results indicate that a clear lake can tip to a turbid state directly, whereas the transition in the other direction goes via coexistence states in which part of the lake is turbid.

Discussion
Classically, tipping in complex systems has been illustrated using simple conceptual models with two alternative states [1][2][3]20].When tipping occurs in these systems, that leads to a reorganization of the full system, with potentially large impacts on its functioning.In this study, we have investigated a prototypical tipping system and extended it by making it spatially explicit, via the incorporation of spatial transport and spatial heterogeneity.The resulting spatially extended model system revealed a more intricate bifurcation structure with additional branches of stable states compared to the classic 'S'-curve bifurcation diagram.This is due to the presence of stable coexistence states with part of the spatial domain in one state and part in another.Further, it revealed the possibility of transitions in which the system only undergoes reorganization in part of the spatial domain, limiting the impact of these events on the full system's functioning.Therefore, this study further specifies how tipping of spatially heterogeneous systems happens; they do not necessarily tip fully in one event, but can have a more fragmented tipping in smaller steps, each corresponding to a restructuring in only part of the spatial domain.These events are therefore less severe than full tipping events, and lead to smaller hysteresis loops.On top of that, while part of the domain is still in its original state, restoration might also be easier and can happen more gradual: as climatic conditions would improve again, the spatial interface between states can move, slowly recovering the system.This detail of the tipping process might be relevant for many natural systems.Here, we have illustrated this using a few conceptual example systems on different spatial scales, but there are many more.In principle, if a model or real system would allow for bistability locally, spatial effects could lead to the formation of coexistence states-and thus tipping in these systems might be fragmented.One of the prominent examples is the occurrence of multiple convection patterns in ocean general circulation models [66,67].Such behaviour was often attributed to numerical artifacts (due to the coarse horizontal resolution in such models).However, the different convection pattern could in the framework here be interpreted as different stable coexistence states.Indeed, the vertical coupling is very strong and follows basically the box model example in section 3.2, whereas the horizontal diffusive coupling is rather weak.
Coexistence states in this study have been characterized by a spatial interface that spatially separates alternative states.Such spatial structures have been observed in reality for systems in which these spatial interfaces are distinct and clearly visible by eye.Examples include the ice grounding line [22] and the tropical forest-savanna boundary [68][69][70].In particular, empirical studies of the latter have made connections to the here described mathematical theory [68,70].Together with the results on conceptual models presented in this paper, these suggest the relevance of coexistence states for many real climate subsystems.It is difficult to study the associated multistability and resulting fragmented tipping empirically, because of the required temporal and spatial resolution.Nevertheless, at least in more realistic models, the behaviour of e.g. an ice grounding line seems consistent with a fragmented tipping scenario: the interface moves slowly most of the time, with some sudden events in which interface changes and a larger part of the domain undergoes transition [22].In addition, there are some early studies suggesting some historic climatic changes, such as atmospheric warming, could be step-like [71], which could be consistent with a fragmented tipping scenario.
As illustrated, the precise importance and relevance of coexistence states depends a lot on the size of the system and the specific heterogeneity.For instance, if a system is too confined, coexistence states cannot form at all, and spatial heterogeneities change the bifurcation structure and the severity of critical transitions.Consequently, for some systems these insights might be more important than for others.Future research should distinguish between climate subsystems that can show fragmented tipping and those that do not.This difference is important to make, as a fragmented tipping scenario indicates an imminent tipping event does not lead to a full system collapse but still leads to similar problems locally in the regions that change state.Of particular interest are the implications on a planetary scale, as it has been hypothesised planetary tipping might happen in the future [17][18][19].If tipping on a planetary scale is fragmented its effects might be limited globally-though they can still be harsh locally.
One of the most important tasks remaining is therefore to distinguish between the type and severity of a tipping event before it happens.Current generic early warning signs seem focused mostly on predicting when an imminent system change is going to happen [28,72,73], but not on what is happening then.If it is only a reorganization in a minor part of the spatial domain, it is not necessarily worrisome.Further, in spatially-extended systems early warning signs, like critical slowing down, might only be visible in the part of the domain in which change is going to happen, while it might be absent in the global response.Both of these issues might be tackled by inspection of the spatial structure of the destabilizing perturbation [25].
In this paper, we have deliberately used conceptual models that are very basic and relatively simple to showcase the general phenomenon and argue its generality.In more complex and realistic models, including non-gradient systems, coexistence states and fragmented tipping should also be present in general.However, dynamics of such models can be much more difficult.For instance, models with 2D or 3D spatial domains, the spatial interfaces can have complicated structures, and in systems with multiple components might also undergo bifurcations [37,61,74,75].
Finally, we have focused on stationary states, but also the transient behaviour of coexistence states is relevant.When not in equilibrium, a spatial interface slowly moves towards its equilibrium position, hereby converting part of the domain from one state to another-but slowly [36,37].This is in contrast with the tipping events, in which change is fast.Together, this leads to a potentially complicated transient behaviour, which can consist of long periods of slow behaviour interspersed with few fast tipping events [25,33].If we would have a better understanding of this transient behaviour of the spatial interface, and how it is influenced by climate change and (human-made) heterogeneities, it might be possible to create circumstances that allow these systems to slowly restore themselves naturally more easily.
Summarising, the take-home messages of this work are as follows.First, the bifurcation structure of any spatially heterogeneous model can be severely more complicated than that of a non-spatial model because of additional branches of coexistence states.Second, the possibility of these coexistence states indicates that tipping in these systems can be fragmented-and thus more gradual-with multiple smaller steps instead of one large tipping event.Third, restoration of a system that has tipped to a coexistence state is easier than from a fully tipped state, because restoration can happen gradually by movement of the spatial interface and because hysteresis loops are smaller.Fourth, the specifications of both the spatial size and the precise spatial heterogeneity determine the precise response of a spatially extended system.This includes to what extent coexistence states and the described fragmented tipping scenario are relevant for the system under consideration.Fifth, examples of many climate subsystems, ranging from global to local scales, suggest coexistence states can emerge on many spatial scales.All this together implies that tipping in many climate (sub)systems might be fragmented and that these systems might be more resilient than indicated by non-spatial models.Level curves for a typical Hamiltonian system for various values of parameter µ.The circles correspond to fixed points of (A.1), and hence to spatially uniform solutions of (3).Blue circles denote the lower stable state, red circles the higher stable state and black, non-filled circles the unstable uniform state.The thick magenta lines in (c) indicate the spatial interfaces (also called fronts, or heteroclinic connections) that are the building blocks of the coexistence states.leading order construction per region.We note that a rigorous extension of the ideas formulated here would require higher-order terms and matching of the regions.See e.g.[77, chapter 2] for more general explanation on the here employed matched asymptotics method.
In the outer region, (A.4) becomes at leading order the equation with the added constraint that y needs to be smooth.That is, in the outer regions, the diffusion term does not play a role at leading order.Hence, this leads to y(x) = y A (x) or y = y B (x).
In the inner regions, we use the coordinate change ξ := x−xF ε , where we zoom in on the position x = x F where the interface is located.We then obtain it is now clear that an interface can only be located at x = x F when ∆H(x F , µ) = 0.

A.3. Spatially heterogeneous equation-stability of equilibrium solutions
For any given stationary state y = y * , the (linear) stability can be analyzed by setting y(x, t) = y * (x) + e λt ȳ(x) in (4), which yields the linear stability problem Again, this equation can be analyzed per region.

Appendix B. The small domain limit
If the spatial domain is too small (compared to the diffusion strength), spatial interfaces do not fit in the domain and coexistence states cannot form.If the domain is really small, the spatial heterogeneity becomes irrelevant and dynamics are essentially similar to a model that has no spatial effects.To see this, consider the following spatially heterogeneous model with one spatial dimension (see also e.g.[77, chapter 5] for a more general treatment of socalled homogenization methods):

Appendix C. Details of the ocean convection model
To model the possibility of bistability between convective and non-convective states, we use a standard box model that models the temperature T and salinity S of the mixed-layer as dynamic variables, which change through the exchange with static atmosphere and ocean boxes [53][54][55][56].We have extended the model by the incorporation of a spatial dimensions and the inclusion of spatial diffusion: where ∆ρ := ρ − ρ 0 is the density difference between the mixed-layer and the deep ocean.Further, T A (x) and S A (x), respectively T 0 (x) and S 0 (x), are the spacedependent (reference) temperature and salinity of the atmosphere and deep ocean box, respectively.k T and k S give the rate of exchange with the atmosphere box and κ(∆ρ) the rate of exchange with the deep ocean box.The latter is a function of ∆ρ, which is a nonnegative increasing function [53].That is, convection happens when the exchange rate κ(∆ρ) is large, which typically occurs if the mixed-layer is heavy enough compared to the deep ocean.Since the equation for ∆ρ does not depend on ∆µ, it can be solved without taking ∆µ into account.Hence, for the purposes here, only the equation for ∆ρ is needed.

Appendix D. Details of the shallow lake model
The shallow lake model that has been created in [15] is non-spatial and dimensional.It is given by the following ordinary differential equation where With spatial effects added (and dropping the hats), this becomes the model in (14).

Figure 1 .
Figure 1.Examples of visibly observable coexistence states and the spatial interfaces between the different states in real systems.(a) Spatial interface between tropical forest and savanna ecosystems (Reproduced from Google Earth.Image stated to be in the public domain.© 2021 Maxar Technologies; Gabon, 1 • 19 ′ 47.15 ′′ S, 13 • 52 ′ 48.66 ′′ E).(b) Spatial Interface between two types of stratocumulus clouds (RAMMB/CIRA SLIDER [30]; 9.45 • S, 73.62 • W, 8 September 2017 16:30:36 UTC).(c) Spatial interface between sea-ice and water in the Eltanin Bay in the Bellinghausen Sea (73 • 43 ′ S, 83 • 49 ′ W, 2 March 2015.Reproduced from NASA's Earth Observatory, NASA.Image stated to be in the public domain).(d) Algae bloom in part of Lake St. Clair on 28 July 2015 (Reproduced from NASA's Earth Observatory, NASA.Image stated to be in the public domain).

Figure 2 . 1 100 ∂ 2 y
Figure 2. Bifurcation diagram for (a) a non-spatial model, (b) a spatially homogeneous model and (c) a spatially heterogeneous model.Solid lines denote stable equilibria and dotted lines unstable ones.The vertical axes in (b) and (c) indicate the spatial average of the state variable.Insets in (b) and (c) show the spatial structure of equilibria for the different branches.The green line in (b) indicates the Maxwell point of the bifurcation parameter.The cyan branch in (c) consists of stable coexistence states.For (c), the model used is ∂y ∂t = 1 100

Figure 3 . 1 100 ∂ 2 y
Figure 3. Showcase of some variations of the bifurcation diagram for spatially heterogeneous models with different spatial heterogeneities.The insets show the spatial structure of equilibria for the different stable branches as a function of the spatial coordinate x.The model used is ∂y ∂t = 1 100

Figure 4 .
Figure 4. Bifurcation diagrams for different diffusion strengths, but the same spatial heterogeneity and domain size, demonstrating the change of structure of the diagram as the domain becomes too small to fit in the spatial interfaces (as more space is required when D is larger).The model used is ∂y ∂t = D ∂ 2 y ∂x 2 + y(1 − y 2 ) + µ + 1 2 cos(πx) on [−1, 1] with no-flux boundaries.

Figure 5 .
Figure 5. Bifurcation diagrams for the example systems in section 3.In all figures, the red and blue branches denote equilibria in which the whole system is in the same state.The cyan branch denotes the stable coexistence states, in which part of the spatial domain is one state and the rest in the other.The insets show the spatial structure of solutions along the stable branches; axes of different insets of a figure are the same and the dashed lines in the insets in (a) correspond to T = 273.15K (i.e. the melting point of ice).(a) Global energy balance model (GEBM).(b) Ocean convection model.(c) Atmospheric boundary layer model.(d) Tropical forest model.(e) Model for turbidity in shallow lakes.

Figure A. 1 .
Figure A.1.Level curves for a typical Hamiltonian system for various values of parameter µ.The circles correspond to fixed points of (A.1), and hence to spatially uniform solutions of (3).Blue circles denote the lower stable state, red circles the higher stable state and black, non-filled circles the unstable uniform state.The thick magenta lines in (c) indicate the spatial interfaces (also called fronts, or heteroclinic connections) that are the building blocks of the coexistence states.

Figure A. 2 .
Figure A.2. Sketch of a coexistence state in spatially heterogeneous models.Part of the domain resides in state A (blue), and part in state B (red), with a spatial interface (cyan) connecting these regions.

) A − cA 2 , 2 )
We define μ := N N+hN and rescale Â = A/h A , ĉ = h A c to obtain the scaled model ≫ 1.Hence, solutions can be approximated using perturbation techniques.Specifically, we can sety(x, t) = y 0 (x, t) + Owhere the integral over the spatial derivatives vanishes because of the no-flux boundary conditions.Since y 0 is uniform in space, the remaining integral is only the spatial average over the heterogeneity.That is, d dt ⟨y⟩(t) = f (⟨y⟩, t) := 1 2 ˆ1 −1 f (⟨y⟩, x; µ) dx (B.8)which has dynamics similar to a non-spatial model.