The topology of non-linear global carbon dynamics: from tipping points to planetary boundaries

We present a minimal model of land use and carbon cycle dynamics and use it to explore the relationship between non-linear dynamics and planetary boundaries. Only the most basic interactions between land cover and terrestrial, atmospheric, and marine carbon stocks are considered in the model. Our goal is not to predict global carbon dynamics as it occurs in the actual Earth System. Rather, we construct a conceptually reasonable heuristic model of a feedback system between different carbon stocks that captures the qualitative features of the actual Earth System and use it to explore the topology of the boundaries of what can be called a ‘safe operating space’ for humans. The model analysis illustrates the existence of dynamic, non-linear tipping points in carbon cycle dynamics and the potential complexity of planetary boundaries. Finally, we use the model to illustrate some challenges associated with navigating planetary boundaries.


Introduction
Increasing human pressure on ecosystems coupled with global change has raised concerns about human societies approaching boundaries or 'tipping points' that, once crossed, may induce fundamental shifts in Earth System dynamics [1,2]. The Earth System has undergone dramatic shifts in the past [3] and some of these appear to have been critical transitions associated with the crossing of a planetary boundary [4]. Motivated by the possibility that human activities may induce a critical transition in Earth System dynamics, Rockström et al [5] have identified nine planetary Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. boundaries and provided an estimate of the position of the present Earth System relative to these boundaries. The authors suggest that three boundaries (rate of biodiversity loss, biogeochemical flows of nitrogen, and atmospheric CO 2 concentration) may have already been transgressed and that we are nearing several others.
The notion of planetary boundaries is a powerful starting point in defining a safe operating space (SOS) for humanity. In the seminal treatment of Rockström et al [5], the planetary boundaries are static estimates. In reality, these boundaries interact dynamically: approaching one boundary may cause others to shift. How these boundaries interact and respond to changes in the Earth System state is a difficult and important question. The work presented here is an initial effort to generate insights into the fundamental features of interacting planetary boundaries.
Models are an essential tool in understanding climate dynamics. Presently, much effort is directed towards the development of highly detailed, complex computational climate models in order to predict the impact of specific human activities on the future climate e.g. [6]. There have been earlier attempts to estimate boundaries for global systems using such models [7] and the computational and data demands are substantial. The challenges associated with estimating even a single boundary suggests that complex models may not be a practical starting point for exploring how boundaries interact. On the other hand, one-dimensional models such as Budyko-Sellers type energybalance models [8][9][10] that have been used to clearly elucide thresholds in the climate system and the mechanisms behind them are too simple to explore the interaction of planetary boundaries. For our purposes, we need a model that produces sufficiently rich dynamics to gain insights into interacting boundaries and that is simple enough to allow rigorous analysis.
In their discussion of climate models in The Fourth Assessment Report of the IPCC, Randall et al [11] advocate the use of a spectrum of models of differing levels of complexity, each being suited for answering specific questions. We adopt this strategy to develop a relatively simple heuristic model specifically designed to explore the relationship between the atmospheric and land system change planetary boundaries. As such we consider only the most basic interactions between terrestrial, marine and atmospheric carbon stocks. The goal of the model is not to predict but, rather, to construct a conceptually and qualitatively correct model of a feedback system like the Earth System. We use the model to characterize the topology of basins of attraction for stable climate configurations and assess their sensitivity to what Rockström et al refer to as controlling variables.
These basins of attraction for stable climate configurations define SOSs where humanity will not experience unacceptable climate change. It is important to note that we are not attempting to review the entire ensemble of boundaries identified by Rockström et al. Specifically, while Rockström et al identify nine Earth System processes: (1) climate change, (2) rate of biodiversity loss, (3) phosphorous and nitrogen cycles (which are interlinked), (4) stratospheric ozone depletion, (5) ocean acidification, (6) global freshwater use, (7) change in land use, (8) atmospheric aerosol loading, and (9) chemical pollution along with their associated thresholds (boundaries), we address only those in italics. Although the planetary boundaries are described in terms of individual thresholds and separate processes, they are, in fact, tightly coupled. For example, if we think of just climate change and rate of change in land use, Rockström et al provide estimates for just two points along the boundary of the safe operating space in these two dimensions. Our goal is to develop a model to connect these points and explore the topology of the full boundary. As such, we do not use the model to conduct a direct analysis of the interactions between these boundaries such as how much transgressing one boundary affects the location of another. Rather, we use the model to conduct Earth System analysis sensu Schellnhuber [12,13] to understand some general features of the boundaries of the SOS related to land-use change (percentage of global net ecosystem production appropriated for human use), climate change (atmospheric carbon dioxide concentration), and ocean acidification (the ability of the ocean to take up carbon dioxide).
Our main contribution is to provide a general, topological picture of the dynamics of three planetary boundaries using a model that, in its entirety, is documented and transparent.
The key advantage of our heuristic model over its more complex cousins is that it allows us to use bifurcation analysis to rigorously and concisely analyze and visualize the dynamic interaction of planetary boundaries. The model can be used to help organize questions, explore interactions, and identify qualitative features of planetary boundaries that can eventually be addressed using large data sets or more realistic models as well as to inform more complex models [14,15,11].

The model
Our Earth System model, similar in spirit to that of Lenton [14], captures interactions between net ecosystem production (NEP) and carbon stocks. Our modeling exercise departs somewhat from more empirically based efforts (see e.g. [14]) in that we emphasize only the essential shape (i.e. topology) of empirical relationships. This allows us to preserve the qualitative dynamics while focusing on implications of basic relationships that are robust to a wide range of parameter choices. All essential details of the model are summarized in figure 1. The model is based on the most basic theoretical considerations, with parameters chosen to broadly match observed empirical patterns (see the appendix for further theoretical and empirical details).
The model is a simple mass balance that tracks carbon stocks (colored boxes in figure 1) and fluxes (solid arrows) using the system of ordinary differential equations (1)-(3) in figure 1 where c t (t), c m (t), c a (t), and C(t) represent terrestrial, marine, atmospheric and total available carbon stocks at time t, respectively. By available carbon, we refer to that which can move via biogeochemical processes between c t (t), c m (t), and c a (t). C(t) is increased by the release of geological carbon, possible only through human activities represented by the industrialization process, I(t) (equation (3), figure 1), which the Earth System then distributes via biogeochemical feedbacks that control fluxes (dashed arrows in figure 1). For this definition of C(t), conservation of mass implies c t (t) + c m (t)+c a (t) = C(t). We can thus represent the dynamics with three state variables, e.g. c t (t), c m (t), and C(t), and the fourth, e.g. c a , will be determined by the mass conservation relation.
The change in c t (t) (equation (1), figure 1) depends on NEP and human activities, H(t), that affect natural processes related to terrestrial carbon release (e.g. fires, tilling soil, land-use change, etc). NEP depends on the per-unit plant biomass photosynthesis and respiration rates, p(T, c a ) and r(T), as well as c t (related to the plant biomass available to photosynthesize and respire). The functions p(T, c a ) and r(T) depend on global mean temperature, T, as shown by the Figure 1. Complete Earth System model (stock data from [16]). Colored boxes indicate stocks, solid arrows indicate flows, dashed lines indicate feedback effects on carbon flows (i.e. the RHS of the differential equations at the top right). Feedback relationships are shown graphically (three graphs below the stock and flow diagram) and mathematically on the right. The variables c t (t), c m (t), c a (t), and C(t) represent terrestrial, marine, atmospheric, and total available carbon stocks at time t, respectively. Parameters marked with † are varied in the analysis. See appendix for further discussion of the theoretical and empirical basis for the functional forms. solid and dashed curves in graph A, figure 1, respectively. The dependence of T on c a (graph C, figure 1), completes the feedback loop (through graphs C and A) between atmospheric carbon, r(T), and p(T, c a ). The direct dependence of p(T, c a ) on c a is due to fertilization effects (feedback through graph B). When c a = 0, p(T, c a ) must be zero, independent of T, as there is no carbon in the atmosphere for plants to extract. The photosynthetic rate increases with c a but saturates as other limiting factors become important. This leads to the functional form in graph B. The photosynthetic rate is then computed by multiplying the temperature effect (graph A) by the fertilization effect (graph B) as indicated by the '⊗' symbol in figure 1. Note, all these rates have been rescaled as proportions of some maximum and thus range between 0 and 1.
Terrestrial carbon is partitioned between aboveground organisms and soil stocks. Carbon in the soil is primarily derived from detritus from aboveground organisms and fine root inputs belowground and leaves the soil mainly through the decomposition of organic matter [17]. On the time scales of interest for our model, aboveground and belowground carbon pools equilibrate relatively quickly, and we approximate this partitioning process as a simple proportion. That is, carbon is partitioned between aboveground and belowground pools according to the proportion α ag ∈ [0, 1] so that aboveground biomass is given by α ag c t . This simplifying assumption reduces the number of state variables and makes the analysis more clear. We discuss implications of relaxing this assumption in the conclusions. Given this assumption, carbon accumulation is, by definition, NEP = r g [α ag c t ] [p(T, c a ) − r(T)] where r g is an ecosystem-dependent conversion factor between metabolic rates and standing carbon, i.e. the growth rate. We know that photosynthesis and respiration depend on other factors in addition to air temperature and atmospheric carbon. For example, plants compete for resources such as water and soil nutrients. To capture this fact, we assume that the growth rate decreases with increasing competition, i.e. is density dependent. The simplest, biologically meaningful and widely accepted way to add density dependence is to assume that the growth rate decreases linearly with biomass density, i.e. r g (c t ) = r tc [1 − α ag c t /k] where k is the terrestrial aboveground biomass carrying capacity and r tc is the maximum growth rate. Combining the effects of temperature, atmospheric carbon, and competition for other nutrients yields  figure 1. The ocean 'solubility pump' removes atmospheric carbon as air mixes with and dissolves into the upper ocean. This dissolved inorganic carbon is then fixed by phytoplankton. We capture the movement of carbon across the atmosphere-ocean interface using a simple model based on Le Chatelier's principle by which diffusion processes will tend to push the carbon dioxide concentrations towards the equilibrium specified by Henry's Law. A simple expression that captures this dynamic is D(c a , c m ) = r d [P(c a ) − K H,pc C d (c m )], where P(c a ) is the partial pressure of CO 2 in the atmosphere as a function of c a , C d (c m ) is the concentration of dissolved CO 2 in the upper ocean as a function of c m , K H,pc is the appropriate Henry's constant, and r d is the concentration-difference dependent diffusion rate. A stoichiometric argument and unit analysis allows us to write P(c a ) = β 1 c a and K H,pc C d (c m ) = β 2 c m where β 1 and β 2 are conversion factors. This leads to the expression D(c a , c m ) = r d [β 1 c a − β 2 c m ]. This can be further simplified to D(c a , c m ) = a m [c a − βc m ] where a m and β are lumped parameters (see the appendix for details). This is the simplest possible representation that correctly captures the essential elements of the process by which carbon moves between the atmosphere and upper ocean. Note that we consider carbon dynamics only in the upper ocean because carbon dynamics between the upper and lower ocean occur on very slow time scales relative to those of interest in our model (again, see the appendix for further discussion). Finally, it is important to note that the rate constant, a m , could decrease with increasing c m due to ocean acidification: as c m increases, the ocean is becoming more acidic which, in turn, may reduce its capacity to take up carbon. This has important implications for planetary boundaries (see section 4).

Analysis I: pre-industrial Earth System dynamics
Our analysis consists of using numerical bifurcation techniques to rigorously compute the boundaries of the basins of attraction for different long-run attractors of the model. Our analysis builds on similar efforts that have focused on how feedbacks between key processes affect the overall stability of different long-run attractors of the Earth System. This includes the energy-balance models mentioned in the introduction as well as more recent work that explicitly considers the feedbacks between biological productivity and weathering that generate two stable Earth System states; one with and one without complex life [18,19]. Although the equilibria and processes that generate them are different than those in the models of Lenton et al [18] and von Bloh [19], our model will also converge to one of two attractors of interest: c t = k (desirable) and c t = 0 (an undesirable, global desert state) that depend on choices of parameter values and initial conditions. We use the model to explore how two features of global change affected by human activities interact to affect the boundary between these attractors: (1) the release of trapped geological carbon through the industrialization process, I(t), and (2) effects on c a and c m through the release of terrestrial carbon stocks. The simplicity of the model allows us to present the results comprehensively and compactly in figure 2. This is the main strength of our approach. Finally, we use the model in a scenario mode to illustrate the difficulty of adaptively responding when planetary boundaries are approached.
To begin, we treat C as a parameter (i.e. I(t) = 0) and explore how varying C and human terrestrial biomass offtake rate, α, affects the topology of the Earth System model. This reduces the system to a two-dimensional model in (c t , c m ) phase space, allowing for a clearer representation of model topology. In section 4 we relax this assumption and allow C to vary over time. We rescale C so that its default value of '1' corresponds to roughly 4500 Gt of carbon with 750, 2850, and 900 in the atmosphere, terrestrial systems, and surface ocean, respectively [6,16]. These choices of units have no effect on the dynamics and are meant only to orient readers. Likewise, the temperature scale is arbitrary-it is just an intermediate quantity that links atmospheric carbon stocks to photosynthesis and respiration rates. The temperature scale shown in graph A, figure 1 was selected based on empirical patterns reported in the literature. Given the uncertainty around the relationship between c a and T [20], we assume it is linear with slope a T , intercept b T (graph C, figure 1). Finally, the parameters r tc , a m , and α simply set the time scale; it is only their relative values, especially with regard to α (which we vary in our analysis) that are important determinants of the dynamics. See figure 1 for a summary of parameter values used in the analysis. Figure 2 summarizes the dynamics of the model. Panel (A) is a phase plane diagram in c m -c t state space. The green and red lines are the c m and c t -isoclines, respectively. The slanted c t -isocline is the locus of points where the photosynthesis and respiration curves intersect: to the right of this line, the terrestrial system accumulates carbon and to the left, it loses carbon. The vertical c t -isocline is located at the terrestrial carbon carrying capacity. The c m -isocline is the locus of points for which D(c a , c m ) = 0, or equivalently, c m = c a . Two trajectories are shown for different initial conditions. Panels (D) and (E) show time trajectories for phase plane trajectories 1 and 2, respectively. There is one, globally stable equilibrium where the isoclines intersect (black circle in panel (A)) where c t = 0.6 (2700 Gt), c a = 0.2 (900 Gt), and c m = 0.2 (900 Gt) so the system can always return to this equilibrium (roughly present conditions) from any perturbation. However, large perturbations (trajectory 2), lead to long periods (panel (E)) with almost zero terrestrial carbon (sustained hot desert) before recovering in contrast to the quick recovery in panel (D). The region in which the system recovers quickly (the SOS) is shaded light gray. Panels (B) and (C) in figure 2 illustrate how increasing human appropriation of terrestrial carbon (α) causes the SOS to shrink (the size of the gray region and the horizontal distance from the equilibrium to the lower branch of the c t -isocline both decrease). For α = 0.3, although a trajectory outside the SOS (type 2) would take a long excursion through a hot desert phase with c t near zero, the system will eventually recover (as in panel (E)) because the equilibrium point (c t , c m , c a ) = (0, 0.5, 0.5) is an unstable saddle. When α increases to 0.5 (panel (C)), however, this equilibrium becomes a stable node (a saddle-node bifurcation occurs between α = 0.3 and 0.5). Now there are two stable equilibria, one with c t > 0 and one with c t = 0 along with their respective basins of attraction (labeled 'PE' and 'ZE' in panel (C), respectively). In contrast to panels (A) and (B) where all trajectories eventually converge to the c t > 0 equilibrium and the SOS boundary distinguishes between those that converge quickly or after a long delay, in panel (C), the SOS boundary distinguishes between trajectories that converge quickly to the c t > 0 equilibrium or those that never do (type 2 trajectories are absorbed into the c t = 0 hot desert state). Combining these two scenarios, we can define the SOS as the collection of initial conditions in c m -c t state space for which the system does not take a very long excursion (possibly infinitely long in the case of panel (C)) through a hot desert phase through which humans would likely not survive. Although it is unlikely that humans could drive c t to zero, other biological factors may generate thresholds beyond which terrestrial systems cannot recover, locking the system in the hot desert state. (type 2 paths as in figure 2(C)). Thus, there is a second layer to the SOS boundary that operates on a slower time scale. The SOSs in panels (A)-(C) are slices of a higher-dimensional object. On shorter time scales, the SOS is defined by the interior of these slices as described above. On longer time scales, the SOS is defined by the interior of the hypervolume in the space defined by c t , c m , α and C (moving along the increasing α direction in this hypervolume causes slices in the c t -c m plane to shrink). Although we exploit the separation of shorter (state space) and longer (parameter space) time scales so we can visualize the SOS more easily, this analysis highlights the fact that the SOS boundary can be an extremely complicated object in reality.
If we believe this model and our parameterization, our present state probably lies somewhere in the white circle-quite far from the total carbon-terrestrial carbon appropriation rate boundary. However, the location of the boundary is sensitive to parameter choices. Can we detect when we are approaching a boundary? Figures 2(B) and (C), in which type 1 and 2 trajectories start very close to one another and diverge at a later time, suggest this can be devilishly difficult. Given the time scales involved, the system could be on a type 2 trajectory for a very long time before it would become clear, based only on measurements of the system state, that the system was in the PE region. Figure 2 illustrates how boundaries relevant on longer time scales (parameter space boundaries) interact and influence boundaries relevant on shorter time scales (state space boundaries).

Analysis II: Earth System dynamics with industrialization
The challenge of remaining in the SOS through an industrialization processes is knowing where the boundaries are. Because we do not know the actual topology of the basins in practice, we must somehow detect boundaries based on measurements of the dynamics of the state of the system. Our analysis suggests that the model resembles a ball rolling on a 'climate sheet' in the c t -c m plane with two basins, one good (PE), one bad (ZE). Releasing geological carbon and clearing land changes the shape of the basins and causes the ball to roll. This analogy is useful for two reasons: (1) we can imagine riding the ball and detecting boundaries by how its velocity (=carbon flux in our model) is changing and (2) the ball has momentum, so we have to stop changing the basin topology before we reach the boundary to avoid transgressing it. That is, the climate system is sensitive both to the amount and rate of geological carbon release. We illustrate these points using three model scenarios. Figures 3(A) and (B) show a scenario in which C(t) mimics Earth's industrialization since 1900. Assuming industrialization is the only source of additional carbon, we have C(t) = C o + C r (t) where C o is the total carbon stock when industrialization begins and C r (t) is the carbon released by industrialization. The dynamics of C r (t) are modeled as a logistic-i.e. I(t) = dC r (t)/dt = r i C r (1 − C r /c max ) where r i sets the pace for industrialization and c max sets the total carbon released when a global agreement curbing emissions is reached (see caption for parameter values). The scenario begins at a pre-industrial equilibrium with c m = c a = 621 Gt, and c t = 3257 Gt (0.138 and 0.724 in non-dimensional variables). We assume society curbs carbon emissions by 2100 when c a (t) stabilizes at roughly twice the pre-industrial level ( figure 3(B), c max = 0.4). Figure 3(A) shows how the marine and terrestrial systems process this signal in terms of net carbon fluxes. These scenarios are a good qualitative match with those generated by a range of other more detailed models [21,22]. Early in the scenario, increasing c a increases photosynthesis through the temperature and fertilization effects and, in turn, increases terrestrial carbon uptake (terrestrial system is a net sink). As temperature increases over time, however, plants become less efficient (around 2025) and flux eventually returns to zero at about 2075. The ocean continues to act as a sink for more than a century beyond that point, enabling the Earth System to return to the PE equilibrium.
The essential dynamic in the model is the relative rates at which the ocean can absorb atmospheric carbon and the terrestrial system loses its capacity to absorb carbon. For the parameter choices for the scenario shown in figure 3(A), we can increase c max up to 0.71, tripling c a and c m before the system moves into the PE region. However, although recent trends in carbon fluxes suggest the Earth is far from a planetary boundary in atmospheric carbon, several processes may affect parameter values. Examples include ocean acidification (decreases a m ), land-use change that releases carbon more rapidly (increases α), and carbon could be released at a faster rate than in the scenario shown in figure 3(A) as large economies develop rapidly before total emissions are curbed (increases r i ). Figures 3(C)-(F) show the potential impacts of these factors.
For example, with a 50% higher industrialization rate (r i = 0.06), all else being equal, society must curb emissions at a 3.5% lower level to avoid moving into the PE region than in the slower industrialization scenario. Although this 'momentum' effect is small, it could become important if the system is near a boundary; the rate at which the boundary is approached causes it to move. If more rapid industrialization is combined with ocean acidification (a m reduced from 0.04 to 0.01) and increased land clearing (α increased from 0.1 to 0.15), the SOS is much smaller. Figures 3(C) and (D) shows a scenario in which society curbs emissions at 30% above present stocks (c max = 0.3). As in the first scenario, as c a increases, both the terrestrial and marine systems become sinks. However, because the marine system takes up carbon more slowly, c a increases beyond the intersection of the photosynthesis and respiration curves in figure 1 (around 2050) when the terrestrial system becomes a net source. The marine system remains a sink, eventually pulling down c a , allowing the terrestrial system to accumulate carbon again (about 2080), which restabilizes the system with c a about double its pre-industrial level. Remaining in the SOS through the industrialization process hinges on the marine system pulling c a down in time for the terrestrial system to recover its capacity to accumulate carbon. (D)-after losing carbon while the marine system absorbs it, the terrestrial system recovers quickly when temperature conditions permit and reduces c a sharply around 2100. For c max = 0.409 26 (dashed), on the other hand, the terrestrial system never recovers. Rather, after losing carbon at a stable rate between roughly 2050 and 2100, it begins to rapidly lose carbon just after 2100 causing both c a and c m to increase rapidly and c t to decrease to zero. Two features might be used to detect the boundary in real time using terrestrial flux data: (1) a long period of stable flux to the atmosphere (curvature near zero, and 2) a change in curvature from convex down to convex up. Far away from the boundary, on the other hand, the curvature is strongly negative during periods of positive flux. Although clear in a deterministic model, noise and complexity in the real system could make the boundary very difficult to detect.

Conclusions
Our goal has been to expand the planetary boundaries concept from a static notion of a collection of fixed boundaries to the context of coupled, non-linear dynamics where these boundaries interact fluidly. Using a stylized Earth System model, we illustrate planetary boundaries using the complementary concept of the safe operating space (SOS). The SOS is equivalent to a region in state space of the dynamic model on short time scales (figures 2(A)-(C)) and to regions in parameter space on longer time scales (figure 2(F)). We have invested considerable effort to capture the main processes and feedbacks in a climate system while keeping the model as simple as possible for two reasons. First, we provided a fully documented model (figure 1) in which all assumptions are transparent. Second, the simple model allowed us to rigorously compute boundaries, explore how they change as parameters are varied, and present the results compactly.
A general feature of non-linear systems is the potential for sensitivity to initial conditions. Even with a complete understanding of the model, without perfect knowledge of the initial conditions (capacity for perfect measurement precision in reality), long-run behavior is impossible to predict. Thus, although we can compute boundaries rigorously in theory, their utility is diminished by measurement limitations in reality. By combining bifurcation analysis and scenarios, we illustrated how interactions among boundaries can generate multidimensional tipping points (figures 3(E) and (F)) that are difficult to detect. This suggests that future work on a more detailed analysis of whether, in more complex models, interacting processes may create critical fragility regions containing infinitely many, impossible-to-detect tipping points would be fruitful (e.g. fractal boundaries). A critical question is whether a single, detectable, global tipping point might emerge from the aggregate of many local and regional tipping points. The global tipping point suggested for loss of biodiversity might indeed be an example of such a global property emerging from the aggregate of many lower level non-linear changes [23].
By showing how multiple tipping points move with changes in the state of the system in highly non-linear ways, the analysis adds clarity to our intuition that the notion of clearly defined tipping points based on a single variable of interest may be misleading. We also illustrate that there are two types of boundaries that are controlled by fast and slow variables, respectively. The first boundary type defines a SOS in state space that is more relevant to short-term fluctuations in carbon stocks related to, for example, rapid land-use change or large-scale fires. The second boundary type defines a SOS in parameter space that is relevant to slow variables such as atmospheric carbon concentration and the long-run average of human use of terrestrial carbon. Finally, we show how movement within the SOS in parameter space changes the shape of the SOS in state space and that crossing the planetary boundary in parameter space causes the SOS in state space to vanish.
Perhaps the most important qualitative feature of planetary boundaries revealed by the model is that there are layers of boundaries that interact in very complex ways. This generates a situation in which both the state of the Earth System and the shape of the SOS boundary are changing concurrently, making the detection of the crossing of a boundary very difficult in practice. Further, given that it is well-known that boundaries of basins of attraction can be extremely complex and fractal even in simple, low-dimensional models, we can guess that planetary boundaries would be very difficult to detect in high-dimensional (e.g. the nine dimensions identified by Rockström et al [5]) reality. Moreover, the analysis suggests that planetary boundaries are not likely to be precisely definable. Rather, they are bands, analogous to mine fields, that separate different regions of state space and that are very dangerous places in which to be. The original planetary boundary analysis of Rockström et al [5] accommodated this difficulty by defining a 'zone of uncertainty' within which a threshold was likely lie, even though its precise position could not be determined. In effect, the zone of uncertainty is equivalent to the planetary boundary as a band, and is analogous to the mine field.
Our analysis suggests that the Earth System is far from a tipping point in the photosynthesis/respiration dynamics of land systems, key processes that influence atmospheric carbon concentration and with it, global average temperature. However, a key point from our analysis is that non-linear feedbacks cause thresholds to move. For example, although we illustrate the critical importance of the capacity of the ocean to absorb carbon, we consider only the relatively small carbon pool in the upper layers of the ocean. Over longer time frames (centuries or longer) the transfer of carbon with deeper ocean layers becomes important. This underscores the need to build on recent work (see e.g. [24]) and develop a better understanding of the long-term capacity of the ocean carbon sink to keep up with the rising rate of human emissions. Likewise, a more detailed analysis of the dynamic partitioning between aboveground and belowground terrestrial carbon stocks may shift boundaries, and thus warrants further work. Further, we did not address the potential release of large amounts of methane from melting permafrost [25] or changes in albedo associated with changes in the extent of the cryosphere, two positive feedback processes that also have a major influence on the global average temperature and could move the photosynthesis/respiration dynamics boundary and shrink the SOS. As this research proceeds, other features of the Earth System we have not addressed such as potential changes in thermohaline circulation, biogeophysical impacts of land-use change, and the dependence of H(t) and I(t) on technological and social change must be addressed.
We conclude with the suggestion that future research employ an iterative process using a combination of simple models like that presented here to identify key feedbacks that impact boundaries followed by the use of more complex models to explore these and related feedbacks in more detail to improve the robustness of the results [26]. This iterative process should be part of an ongoing interaction between the policy and research communities to periodically reassess the individual boundaries as their interactions become better understood and as human activities change the position of key control variables on the individual boundaries.

A.1. Temperature relationships
The key relationships considered are between air temperature, photosynthesis and respiration (at the ecosystem level) and between air temperature (without considering complicating geophysical factors) and carbon concentration in the atmosphere. Photosynthesis and respiration are affected by many factors but at the most basic level, studies suggest that both respiration and photosynthesis rates will increase with temperature when temperatures are low and then begin to decrease beyond some physiologically set threshold. These relationships are formalized as follows. Photosynthesis. The basic hump-shaped relationship between photosynthesis (net carbon assimilation rate in µmol m −2 s −1 ) and temperature has been demonstrated for a number of plant species [27][28][29]. Here we use a simple function, to capture this relationship. The parameters a, b, and c, are chosen to fit the general features of temperaturephotosynthesis relationships from the literature. Thus, we let p(T) = f (T; a p , b p , c p ) where p(T) is photosynthesis as a function of temperature, T. Note that this relationship can be specified in terms of leaf temperature (e.g. [28]) or air temperature in a controlled volume during photosynthesis (e.g. [29]). For our purposes, this must be adjusted to an annual mean temperature for consistency with temperature-respiration relationships as discussed below. An example of p(T) is shown by the solid curve in figure A.1(A).
Respiration. Mahecha et al [30] suggest that at the ecosystem level, temperature sensitivity of respiration rate is constant, i.e. Q 10 ≈ 1.4, globally. Thus respiration increases exponentially with increasing temperature. However, over shorter time scales (minutes to hours), Q 10 has been shown to decrease linearly with temperature [31][32][33]. In this case, the temperature-respiration relationship would have the shape of the dashed curve in figure A.1(A) for T < 18. Cf [34]. Although the basic chemical reaction that drives respiration (e.g. [14]) may allow for respiration rates to increase indefinitely, at some point high temperatures may affect the physical system that supports the reaction thereby reducing the efficiency with which the reaction can be carried out. This suggests that the temperature-respiration would bend over at high temperatures as in the dashed curve in figure A.1(A) for T > 18. Recent work by Yuan et al [35] provides evidence that this is indeed the case. Again, the function f (T; a, b, c) defined by (A.1) is sufficiently general to capture this relationship, and we let r(T) = f (T; a r , b r , c r ). Where r(T) is the respiration rate per unit of plant biomass. The dashed curve in figure A.1(A) shows an example in which a r , b r , c r , and the temperature scaling have been chosen to match the basic relations reported in [35]. The temperature axes should be interpreted as mean annual nocturnal air temperature. It is important to note that the details of these functional forms is not critical to the qualitative dynamics for the system. The important feature is that for some temperature range p(T) > r(T). Otherwise, terrestrial ecosystems could never accumulate carbon (biomass), i.e. NEP would always be negative (see colored regions in graph A in figure 1 in the main document for a graphic depiction).

A.2. Atmospheric carbon and NEP
CO 2 affects NEP in two ways: indirectly through temperature (both photosynthesis and respiration) and directly through a fertilization affect. Given the uncertainty associated with changes in radiative forcing associated with CO 2 concentrations and changes in global and annual mean temperature [20], we make only the most basic assumption: annual mean temperature (in the absence of feedbacks between different carbon stocks) increases linearly with CO 2 .
If we define c a (t) as the atmospheric carbon concentration at time t, then we have where, again, a T and b T are arbitrary parameters. Substitution of (A.2) into the expressions for p(T) and r(T) leads to a family of mappings between c a , photosynthesis, and respiration, i.e. Regarding the fertilization effect of atmospheric carbon, the simplest assumption we can make as that increasing carbon increases photosynthetic activity, recognizing the fact that other limiting factors will eventually begin to reduce the effects of additional CO 2 in the atmosphere. This implies that ∂p/∂c a > 0 and ∂ 2 p/∂c 2 a < 0. A very simple function that exhibits such behavior is where 0 < b f < 1 and c f > 0 are parameters. This functional form has the basic shape associated with Michaelis-Menten kinetics assumed in Lenton [14]. Assuming a multiplicative relation between temperature and carbon fertilization effects (again, as in [14]) and writing

A.3. Terrestrial carbon dynamics
Terrestrial carbon dynamics involve multiple stocks in plants and animals aboveground and in soil organisms below ground. Further, carbon dynamics are intimately related to other basic nutrients such as water, phosphorous, nitrogen, etc. Here, we provide additional detail to complement the discussion in the main text. Specifically, we consider the interaction between terrestrial and atmospheric carbon stocks alone to illustrate the importance of other nutrients and density dependence for a reasonable model. Equations (A.7) and (A.8) define carbon flow rates per unit of plant biomass. Because we are working at the global scale and on time scales on the order of a century, we neglect the finer temporal and spatial scale complexities of the terrestrial carbon cycle associated with the movement of carbon between the aboveground and belowground stocks. Assuming that on the relatively large temporal and spatial scales of interest here aboveground plant biomass is roughly proportional to the total terrestrial carbon stock, by definition we have (A.9) as described in the main text. Obviously our model is a gross simplification-both photosynthesis and respiration would depend on ecosystem structure-i.e. the way c t is partitioned across species. On the other hand, our model is operating at the global level, so we assume that r g is a global average conversion factor across all ecosystems. Assuming no human impacts for the moment, the simplest model for the change in terrestrial carbon stores is then This simple model, however, is too limited to produce realistic results. To see why, consider for the moment a system with only terrestrial and atmospheric stocks. Conservation of mass implies that c a + c t = C where c t and C are the terrestrial and total carbons stocks, respectively. Substituting C − c t for c a in (A.10) yields the following first-order, autonomous, non-linear ordinary differential equation Figure A.2 shows the right-hand sides of (A.11) with total carbon normalized to 1. The different colors correspond to different parameter values as discussed below. The main point is that for the values of c t where this curve is below zero, dc t d t < 0, and the terrestrial carbon stock decreases, and vice versa (shown by arrows for red curve, figure A.2(B)). Thus the system tends to move away from the roots shown in squares-i.e. those equilibria are unstable. For initial conditions to the right of the unstable equilibria, the system will converge to the equilibria shown with circles-a state with high terrestrial carbon concentrations, in a cold climate, with a low carbon atmosphere, and very low rates of photosynthesis and respiration (e.g. perpetual tundra with lots of peat bogs and hydrocarbon reserves). For those to  This dichotomy of either frozen tundra or burning desert is obviously unrealistic. NEP is influenced by many other factors than atmospheric carbon alone. Most important among these is competition for other scarce resources such as water, space, and soil nutrients. Describing such processes in an ecologically satisfying way increases model complexity extremely rapidly-and is well beyond the scope and objective of this letter. The most basic elements of the population growth process in ecology are (1) growth rates are limited by physiology when resources are abundant and (2) growth rates eventually become limited by competition as the population grows and resources become scarce. The most basic, widely accepted model in ecology that captures these two processes assumes growth rate decreases linearly as population density increases (i.e. the logistic model with density-dependent growth) as described in the main text. Here, we illustrate the importance of density dependence on the model by analyzing the two-dimensional case given by Note, as mentioned in the main text, we can easily eliminate the parameter α ag by choosing units. We can use units of aboveground biomass or total terrestrial carbon stock to measure the intrinsic growth rate and carrying capacity (because one is proportional to the other just as are pounds and kilograms). Thus, we can replace r tc α ag and k/α ag with r tc and k, respectively while making a note that we are using, for example, kilograms. Examples of dc t d t generated by (A.12) are shown in figure A.3(A). A comparison of figures A.3 and A.2 makes the effect and importance of density dependence clear: it introduces a new stable equilibrium at c a = k. Without biological processes of density-dependent competition for resources other than carbon (or something like this process), the model would not exhibit an intermediate, present-day-Earth-like state between the extremes described above. The resulting Earth System dynamics are illustrated with arrows on figure A.3(A). Now there are four equilibria, two stable and two unstable. With this formulation, we can describe the idea of a planetary boundary in terms of land use. This simple model does not capture the complexity of land use, of course, but may (or slightly more complex models like it) correctly describe the issue of moving boundaries, and an SOS defined by feedbacks. The boundary and corresponding SOS shown correspond to the blue curve (strongest fertilization effect). The system can tolerate reductions in terrestrial carbon (and corresponding increases in atmospheric carbon) to roughly 0.41 (corresponding to 0.59 in atmosphere). Given that the equilibrium is 0.7, we can compute the size of the basin of attraction, or specified resilience as 0.29. Now are now in a position ask how human activities affect the boundary and the resilience of the system.
Finally, we introduce human activity through a terrestrial carbon offtake term, H(t), i.e. clearing, burning, or farming techniques that reduce the carbon sequestration capacity of terrestrial systems. We simply assume that humans release (or reduce uptake by) a proportion of the terrestrial carbon stock in each time period and let H(t) = αc t (t) where αc t is the proportional rate. To illustrate the impact of H(t) we can modify the model in (A.12) to read Figure A.3(B) shows the effect of human carbon offtake. As α increases, the SOS shrinks and eventually vanishes. We can compute when this occurs (peak of curve falls below horizontal-e.g. red curve, figure A.3) using bifurcation techniques. If α exceeds 26% of terrestrial carbon stocks, the system enters the 'hot desert' (HD) basin.

A.4. Marine-atmospheric carbon dynamics
As mentioned in the main text, our model is based on the fact that carbon will diffuse across the atmosphere-ocean interface whenever concentrations deviate from equilibrium. Equilibrium concentrations are determined by Henry's law where P(c a ) is the partial pressure of CO 2 in the atmosphere as a function of c a , C d (c m ) is the concentration of dissolved CO 2 in the ocean as a function of c m , and K H,pc is the appropriate constant of proportionality. If P(c a ) > K H,pc C d (c m ), carbon will tend to diffuse from the atmosphere to the ocean and vice versa. As discussed in the main text, we use a linear model to relate the rate of diffusion to the size of the deviation from equilibrium-the simplest model possible that captures the essence of this dynamic. We compute P(c a ) by assuming that the atmosphere behaves as an ideal gas. Thus, P(c a ) = xP T where x is the mole fraction of carbon dioxide in the atmosphere and P T is the total pressure. Assuming P T = 1 atmosphere, we have P(c a ) = x. By definition x = n CO 2 n CO 2 + n oth (A.15) where n CO 2 is the number of moles of CO 2 in the atmosphere and n oth is the total number of moles of other gases in the atmosphere. Because n CO 2 n oth , we can approximate (A.15) as x ≈ n CO 2 n oth .
(A. 16) Stoichiometry dictates that n CO 2 = n C which allows us to write n CO 2 = c a × 10 15 MW c (A. 17) where c a is the mass of carbon in the atmosphere measured in gigatons (the 10 15 term converts gigatons to grams) and MW c is the molecular weight of carbon. Given that the mass of the atmosphere is approximately 5.15 × 10 21 g we have n oth ≈ 5.15 × 10 21 MW oth (A. 18) where MW oth is the weighted average molecular weight of the atmospheric gas mixture (≈29 g mol −1 ). Thus, x ≈ c a × 10 15 × 29 5.15 × 10 21 × 12 = 5 × 10 −7 c a . (A. 19) These calculations show that P(c a ) can be approximated reasonably well, as claimed in the main text, by P(c a ) = β 1 c a (A.20) where β 1 ≈ 5 × 10 −7 . This approximation is consistent with the data-i.e. there are 760 Gt of carbon in the atmosphere, so the mole fraction (and, by assumption, the volume fraction) of CO 2 is 760 × 5 × 10 −7 = 380 × 10 −6 = 380 ppm. The relationship between dissolved CO 2 and total carbon stock in the upper ocean is significantly more complex than that between total atmospheric carbon and partial pressure of CO 2 . To capture the relationship faithfully would require a model with several carbon pools by which carbon is tracked as it moves from a dissolved inorganic state, taken up by phytoplankton and moved into the food chain, respired back to the water column, or moved to the deep ocean as falling dead material. We can make two assumptions to simplify this story. (1) The population dynamics of phytoplankton (and thus the movement of carbon from dissolved CO 2 to organisms) occurs on weekly time scales. On time scales of decades and longer of interest here, this process can be viewed as instantaneous. (2) We are interested in the capacity of the ocean to 'pump' carbon out of the atmosphere, and more specifically the rate at which it can do so. Thus, the total amount of carbon in the marine stock (e.g. both upper and deep ocean) is less important than the concentration in the upper ocean (25-200 m) where the key mechanism of mixing and uptake by phytoplankton occurs (which controls the pumping rate). Finally the movement of carbon into the deep ocean is slow compared to our time scale of interest.
Based on assumption 2 we set the net flux from the upper ocean to the deep ocean to zero. Based on assumption 1, we represent the partitioning of carbon between the inorganic and organic states as a simple proportion. This allows us to write: where α di ≤ 1 is the fraction of the total carbon in the upper ocean that is dissolved CO 2 , and V is the volume of the upper ocean. Given that the area of the ocean is roughly 361 × 10 12 square meters, we have Given that K H,pc is order 10, and alpha is order 10 −1 (or smaller), then β is order 1. In order to keep the model analysis a clear as possible, we simply assume that β = 1.