Oscillations in a simple climate–vegetation model

. We formulate and analyze a simple dynamical systems model for climate–vegetation interaction. The planet we consider consists of a large ocean and a land surface on which vegetation can grow. The temperature affects vegetation growth on land and the amount of sea ice on the ocean. Conversely, vegetation and sea ice change the albedo of the planet, which in turn changes its energy balance and hence the temperature evolution. Our highly idealized, conceptual model is governed by two nonlinear, coupled ordinary differential equations, one for global temperature, the other for vegetation cover. The model exhibits either bistability between a vegetated and a desert state or oscillatory behavior. The oscillations arise through a Hopf bifurcation off the vegetated state, when the death rate of vegetation is low enough. These oscillations are anharmonic and exhibit a sawtooth shape that is characteristic of relaxation oscillations, as well as suggestive of the sharp deglaciations of the Quaternary. Our model’s behavior can be compared, on the one hand, with the bistability of even simpler, Daisyworld-style climate–vegetation models.


Background
Climate has an important effect on vegetation. Plant growth is affected by temperature, carbon dioxide (CO 2 ) levels and nutrient availability. The space available for growth matters as well: ice-covered parts of land are not suitable for it. The interaction also works the other way around, though: vegetation plays an important role in climate regulation. Many different effects are at play, with the albedo being one of the most important ones: vegetation is darker than bare ground or ice and therefore absorbs more solar radiation and warms the planet.
This vegetation-albedo feedback appears to be important in semi-arid regions (Otterman, 1974), where it interacts with the hydrological cycle. J. G. Charney and colleagues (Charney, 1975;Charney et al., 1975) were the first to include in a model this biogeophysical feedback, as he called it; many others have followed since (Claussen et al., 1999;Zeng et al., 1999;Zeng and Neelin, 2000). The vegetation-albedo feedback also matters in certain high-latitude regions (Otterman et al., 1984), where boreal forests mask snow in winter, causing an effective warming of the surface Bonan, 2008).
The uptake of CO 2 by plants is another major effect: it attenuates the greenhouse effect and thus cools the surface. This is an example of a biogeochemical feedback. In addition to the obvious effects of albedo change and CO 2 uptake, there are many other mechanisms through which vegetation influences climate. Examples are the effect of phytoplankton on cloud formation -called the CLAW hypothesis, after the initials of the original authors (Charlson et al., 1987;Ayers and Cainey, 2007) -evapotranspiration or more exotic feedbacks, such as the so-called lightning-biota feedback (Shepon and Gildor, 2008). Meir et al. (2006) review several important mechanisms by which terrestrial ecosystems affect climate, while Claussen (2009) provides an introduction to vegetation-climate interactions on paleoclimate timescales.
Although vegetation plays an essential role in the climate system, it has only been rather recently included as an active player in climate models. The hierarchy of climate models (Schneider and Dickinson, 1974;Ghil and Robertson, 2000;Ghil, 2001) ranges from simple, conceptual ordinary differential equation (ODE) models -through intermediate models of varying complexity -all the way up to full-scale general circulation models or global climate models (GCMs). Across this whole range, vegetation can be included to better explain various climatic phenomena and trends. In some cases, the predictions of models that couple atmosphere, ocean and vegetation dynamics -often referred to as Earth system models (ESMs) -differ radically from models that exclude the vegetation (Meir et al., 2006). It follows that it is of the essence to include vegetation in our models to obtain a better understanding of climate evolution and variability.
The simplest climate models are conceptual models, which are usually governed by a small number of ODEs. These models do not claim to be realistic in the sense of making precise quantitative predictions, but allow one to study basic underlying mechanisms. They are also useful for exploring qualitative changes in the climate system's behavior, commonly known in dynamical systems theory as bifurcations (Ghil and Childress, 1987;Dijkstra and Ghil, 2005), and related more recently to the broader and fuzzier concept of "tipping points" . Bifurcations or tipping points correspond to situations in which a small parameter change can have a large effect on the behavior of the whole system.
Studying conceptual models can also provide guidance in interpreting results from larger, more detailed models (Brovkin et al., 1998Ghil and Robertson, 2000). Early applications of dynamical systems theory and of numerical bifurcation studies to the climate system included H. Stommel's bistability study of the thermohaline circulation (Stommel, 1961) and E. N. Lorenz's study of chaos in a simple convection model (Lorenz, 1963). Another such application was to explain Quaternary glaciation cycles as climatic-oscillator solutions (Källèn et al., 1979;Le Treut and Ghil, 1983;Saltzman, 1983;Ghil and Childress, 1987;Ghil, 1994;Crucifix, 2012).
In the latter idealized climate-oscillator models, vegetation is not usually included. There are, however, some simple models that explore the interaction between climate and vegetation, and we will briefly review some noteworthy ones in the next section.

Simple ODE models for climate-vegetation feedbacks
A pioneering model dealing with vegetation and climate was Daisyworld (Watson and Lovelock, 1983). This model illustrates how vegetation may act to regulate planetary temperature through the albedo feedback, for a wide range of parameter values. The varying parameter is in this case the incoming solar radiation. The model has been thoroughly studied and extended; see the review of Wood et al. (2008). von Bloh (1996, 1997) introduced another set of simple, spatially zero-dimensional (0-D) models for vegetation-climate interactions. These highly simplified models include an ODE for temperature evolution, absent from Daisyworld, but look at only one type of vegetation, whereas Daisyworld has two. In their two-ODE model, Svirezhev and von Bloh (1996) find multiple steady states. Such bistability seems to occur across a hierarchy of climatevegetation models, from the simplest (Dekker et al., 2007;Janssen et al., 2008;Aleina et al., 2013) to more complex (Brovkin et al., 1998;Claussen, 1998;Irizarry-Ortiz et al., 2003;Renssen et al., 2003) ones, although it is, of course, harder to ascertain in the latter. In models across the hierarchy, vegetation and temperature are often coupled with precipitation, which provides an additional feedback mechanism.
Beyond the possibility of multiple steady states, it is of interest to examine in the simplest conceptual models the possibility of internal, self-sustained oscillations. Oscillatory behavior has been observed in Daisyworld-like models, for example when an explicit temperature equation is added (Nevison et al., 1999;Fernando Salazar and Poveda, 2009), or when delays are introduced (De Gregorio et al., 1992).
So far, though, simple, 0-D climate-vegetation models seem not to have included an ocean component. Earth's oceans constitute about 70 % of the area of the planet, and are a very important factor in determining its climate (Stommel, 1961;Saltzman, 1983;Ghil, 1994;Dijkstra and Ghil, 2005). The global ocean is often included in 0-D paleoclimate models for glacial cycles (Källèn et al., 1979;Le Treut and Ghil, 1983;Ghil and Childress, 1987;Gildor and Tziperman, 2001;Crucifix, 2012). Our model will include an ocean and the associated sea ice, and will combine certain features of Daisyworld-like models and climatic-oscillator models.
The purpose of this paper is to explore the solution space of such a model and to relate its behavior to other models and to observations. In the next section, we formulate the model, and the numerical results are presented in Sect. 3. The model's behavior and its relation to more detailed models and to observations are discussed in Sect Schematic representation of the model planet's surface, including a fraction p of land and 1 − p of ocean. The fraction of the ocean that is covered by sea ice depends on the global temperature T , while the land is covered by a fraction A of vegetation.

Model description
The climate system contains several subsystems, all working together to produce highly nonlinear behavior through its many feedback mechanisms. Some of the simplest and most important feedback effects act via the planetary albedo: darker areas -like the vegetated ones -absorb more solar energy and thus warm the planet, while lighter areassuch as those covered by snow and ice -tend to cool it. The ice-albedo feedback was included in energy balance models (EBMs) of climate long ago (Budyko, 1969;Sellers, 1969), and a vegetation-albedo feedback regulates Daisyworld's behavior. The latter will constitute the main mechanism in our present model. The model's two governing equations are given below: The variable T denotes global average temperature, while A denotes the fraction of land that is covered by vegetation. Temperature T changes in Eq. (1a) as a result of the balance between incoming and outgoing energy. The parameter Q 0 is the incoming solar energy, which is equal here to 342.5 Wm −2 . The function α = α(T , A) denotes the dependence of the albedo on T and A, and it is given by where p is the fraction of the planet that is land, taken to be 0.3 like on Earth. The values α o , α v and α g denote the albedo of the ocean, of the vegetation and of bare ground, respectively. The latter two are constant, and an essential ingredient of our model is that α v < α g , i.e., that forests and savannas are darker and absorb more energy than bare ground.
The albedo of the ocean will be taken as a function of global temperature T , to account for the possible presence of sea ice. We will use a ramp function, as in Sellers (1969) and Ghil (1976). This function is given by the following equation: Here α max = 0.85 for the ice-covered ocean, α min = 0.25 for the ice-free ocean, T α, = 263 K and T α,u = 300 K. The value of T α,u is rather high -almost 27 • C -which means that a tiny bit of sea ice will be present even for very high global temperature. Figure 1 shows a schematic representation of our planet's surface features and a plot of the ocean albedo versus temperature is shown in Fig. 2. The function R o (T ) denotes the energy outgoing from the planet. Several EBMs have used a modified version of the quartic Stefan-Boltzmann law (Sellers, 1969;Ghil, 1976), while other such models have used a linearization thereof (Budyko, 1969). We will also use here a linear dependence of outgoing radiation on temperature, but radiation increases in it more slowly with temperature than in the Stefan-Boltzmann case. This slower increase takes into account the fact that increasing temperature entrains increasing CO 2 levels and thus an increased greenhouse effect, which tends to decrease the outgoing radiation.
The form of R o (T ) we choose is where B 0 and B 1 are constants, while T opt is the optimal growth temperature for the vegetation. There is a huge uncertainty in the values of these parameters, especially in B 1 , the linearized effect of temperature on radiative forcing (IPCC, 2013), since it all depends on which effects are taken into account and which are not (Zaliapin and Ghil, 2010). We will not attempt to estimate the most realistic values of these or other model parameters -since that is not the role of such simple, conceptual models (Ghil, 2001) -but for definiteness we will use B 0 = 200 and B 1 = 2.5. Their exact value does not play an important role in affecting the model's behavior; still, we note that -over the range of temperatures we are interested in -they correspond rather well to a linearization of an outgoing energy term of the form gσ T 4 with g = 0.55, where g is a constant denoting the grayness of the atmosphere (see Fig. 3). Such a parametrization of outgoing energy -the Stefan-Boltzmann law with a multiplicative factor -is often used, but usually g is allowed to depend on another variable, such as temperature (Sellers, 1969) or CO 2 level (Svirezhev and von Bloh, 1997), rather than being merely a constant, as in Fig. 3 here. Vegetation cover A in Eq. (1b) grows logistically, with a temperature-dependent growth rate β(T ). The shape of the function β follows Daisyworld and related models, to wit, This means that the growth rate is zero, except for a certain interval, in which the dependence is parabolic with a maximum at T opt , which is 283 K in our model. The parameter γ is the mean death rate of the global vegetation A, and its value is, at first, set to 0.1 yr −1 . Table 1 lists the definition of the model parameters and their nominal values. In Sect. 3.2, γ will serve as a control parameter in our exploration of model behavior.

Multiple equilibria and their stability
We start, as usual, by looking at the fixed points of the system (Eq. 1), and will use, again as usual, interchangeably the terms steady states, stationary solutions and equilibria. One obvious fixed point is given by A = A 0 = 0, with temperature T 0 given by the solution of Q 0 (1 − α(T 0 , 0)) = R o (T 0 ). For our parameter values, a solution T 0 only exists when T < T α, and it equals T 0 243 K. This result implies that a planet without vegetation will necessarily drop into a very cold "snowball-Earth" state. Note that this result differs from that of EBMs with no vegetation, in which two stable fixed points co-exist for sufficiently high Q 0 (North et al., 1981;Ghil and Childress, 1987).
A quick check shows that, for p = 0, i.e., for an "aquaplanet", three fixed points are found for A = 0, as in classical EBMs: (i) a snowball Earth -called a "deep freeze" in the 1980s, before geological evidence for it was found; (ii) a high global temperature, comparable to Earth's present climate; and (iii) an intermediate one, between the former two, and considerably closer to the present climate than to the subfreezing one. Hence, it is the inclusion of the relatively high ground albedo for p > 0 that pushes our climate-vegetation model to a low temperature, if no vegetation is present. In the presence of vegetation, A = 0, a fixed point (T + , A + ) must satisfy according to the vegetation equation (Eq. 1b). From the temperature equation (Eq. 1a), we see furthermore that, at a fixed point, one also has The two curves in the phase plane (T , A) of the model that are given by Eqs. (6) and (7) are called the null isoclines -or nullclines, for short -of the system, and they are plotted in Fig. 4. They intersect at two points, aside from the previously discussed bare-land solution; hence, we have two additional equilibria for a vegetated planet.
To study the stability of these fixed points, we consider the Jacobian matrix J = J(T , A) for the system Eq. (1), where For the no-vegetation fixed point A = 0, the eigenvalues of J(T , 0) are a and β(T ) − γ . These eigenvalues are both negative, since T is so low at this fixed point that ∂α/∂T = 0 and β(T ) = 0. Hence, this non-vegetated state is always linearly stable. Let us now look at the other two fixed points, shown in Fig. 4, where A > 0. Note that since A = 1 − γ /β(T ), the temperature must be in the range where β(T ) > 0. With our choice of parameters, this range also lies in the interesting range between T α, and T α,u , within which we have sea ice feedback; it corresponds to the decreasing part of the temperature nullcline, given by Eq. (7), and shown as the dashed red line in Fig. 5.
The leftmost of these two fixed points, where A > 0, is a saddle. This can be seen from numerical computations, but we can also deduce it analytically. In the following, we will use β instead of ∂β/∂T , since the function β depends on temperature alone, and no confusion can arise. We have the following proposition.
The proof is given in Appendix A, and it does not depend on any particular parameter values, only on their relative signs; nor is the exact parametrization of the growth rate β used.
The stable and unstable manifolds of the saddle were also computed numerically, and they are plotted in Fig. 5. The stable manifold acts as a separatrix of the phase-plane flow, and it is almost vertical. This means that initial states with lower temperatures than that of the saddle will lead to the stable steady state (T 0 , A 0 = 0) that is characterized by freezing temperatures and no vegetation, while initial states with higher temperatures will lead to the steady state with at least some vegetation present.
For the stability of the latter, rightmost fixed point, Proposition 1 above says nothing, since β (T ) < 0 there. Numerical computation shows that the eigenvalues of the Jacobian J(T , A) at this point are complex conjugate and equal to λ ± = −0.006 ± 0.028 i. Hence, the fixed point is a focus and the real part λ r = tr(J)/2, where tr(J) is the trace of the Jacobian matrix, is negative; thus the focus is stable.
We can compute the trace of J analytically: The term a on the right-hand side is larger than 0, but the term x = γ − β(T ) compensates and, for the standard parameter values in Table 1, tr(J) is negative, implying the stability of the focus. One suspects, however, that, as some of these parameter values change, the focus may lose its stability. We investigate next this loss of stability and the more complex model behavior to which it leads.

Analytical results
We study now in greater detail the stability of the fixed point with non-zero vegetation (T + , A + ). The result is stated in the following proposition and the proof thereof is given in Appendix A, after the one for Proposition 1.
Proposition 2. A model fixed point (T + , A + ) that lies in the region with active sea ice feedback (T α, < T + < T α,u ) is a focus if and only if the following inequality holds: (1 − y) 2 < 4η(y)y, where y = (γ − β(T + ))/a, η(y) is (up to a sign change) the ratio of the derivatives of the two nullclines and a is given by Eq. (9). Given the inequality (Eq. 11), the focus is stable if y < −1.
The proof in the Appendix relies simply on rewriting the classic conditions on the Jacobian.
We know that y < 0, since x < 0 and a is positive, and also that η(y) is always negative. Since for all negative y we have a necessary condition to have a focus is that |η| > 1, i.e., that the ratio of the derivatives of the nullclines should be large enough in absolute value.
The proof of Proposition 1 in the Appendix implies, however, that this condition is already necessary for a fixed point to not be a saddle. Since the temperature nullcline is a piecewise linear curve, its derivative is a constant and a focus can only occur when the vegetation nullcline is steep enough. This condition, in turn, corresponds to the vegetation cover A + being small enough, as can be seen in the phase plane plotted in Fig. 5. It follows that a Hopf bifurcation and, therefore, oscillations, can only occur around a model state with  Figure 6. Bifurcation diagram, with the vegetation A plotted versus the death rate γ . The branches were computed using the XPPAUT software package (Ermentrout, 2002). low vegetation. Since we did not use the actual value of any of the parameters, this result is quite general.
As shown in the caption of Fig. 5, the temperature value for the vegetated fixed point is T + = 297.9. The corresponding value of y is −1.96. The ratio η(y) of the nullcline slopes equals −9.6; hence |η(y)| > 1, and thus the point could be a focus. Substituting in Eq. (11) yields 8.75 < 75.17, which shows that the inequality is satisfied and that this point is a focus. Since y = −1.96 < −1, the focus is stable.

Numerical results
As announced at the end of Sect. 2, we now vary the death rate γ of the vegetation and follow the change in model solutions. For decreasing γ , i.e., for longer-lived vegetation, the focus loses its stability and the associated Hopf bifurcation gives rise to a limit cycle, as shown in the bifurcation diagram of Fig. 6. This bifurcation occurs for γ = 0.02572, i.e., when the overturning time 1/γ of the vegetation is about 40 years, which corresponds definitely to trees rather than grasses.
We can also see in Fig. 6 that, when γ is higher than 0.41, non-zero vegetation can no longer be a steady-state solution. At γ 0.41, a saddle-node bifurcation gives rise to a stable (blue curve) and an unstable (red curve) branch of steady states, in addition to the stable branch A ≡ 0 studied in the previous subsection and not shown in the figure.
We mentioned in Sect. 3.1 above that the parameter change needs to make tr(J) positive in order to change the stability. This also works if, instead of decreasing γ , we decrease the thermal heat capacity C T . It is the product γ C T that determines the behavior of the system, not the two parameters separately.
This product is a ratio of timescales, with C T a typical timescale for temperature adjustments, while 1/γ is a typical timescale for vegetation. The oscillations occur, therefore, when the ratio of these two timescales satisfies a certain inequality. Such a condition was found to hold for other Hopf bifurcations, where the timescales need to match. For instance, in the temperature-ice sheet oscillator of Källèn et al. (1979), the characteristic time of ice sheet evolution C L should be large enough in comparison with the timescale C T associated with the thermal inertia of the oceans; see Ghil and Tavantzis (1983). Nevison et al. (1999) also obtained certain results that are somewhat similar to ours. Their paper explores a Daisyworld with an equation for temperature evolution. The system exhibits self-sustained oscillations for a certain range of C T and γ , just like in our model. The albedo feedback in this model is due to two types of daisies, black (absorbing) and white (reflecting) ones. The bimodality of the albedos in both models leads us to suspect that the role of sea ice in our model corresponds roughly to that of the white daisies in the model by Nevison and colleagues.
This analogy should not be pushed too far, however. White daisy growth is favored by warmer temperatures, whereas sea ice appears at low temperatures. A comparison between our Fig. 7 and Fig. 2 in Nevison et al. (1999) shows that the behavior of black daisies and temperature is analogous to our vegetation and temperature, while increases in sea ice and white daisies have a different phase. The former occurs with decreasing temperature, while the latter happens as temperature increases.
Another interesting difference between the two models lies in the value of γ for which oscillations are obtained. We have 1/γ 40 yr, while Nevison and colleagues obtain the oscillations when 1/γ is of the order of a few months, with a similar value of C T . The difference probably resides in the nature of the feedback: in our model, the ocean and sea ice play a role, with both operating on seasonal and interannual timescales. In the Nevison et al. (1999) model, it is two kinds of daisies on an all-land planet that interact, both with much shorter timescales.
Oscillatory behavior is also observed in the Daisyworld variant studied by De Gregorio et al. (1992). In their paper, an explicit delay is introduced into the Daisyworld system. When the delay is long enough, oscillations are observed. This again shows the importance of timescales in these oscillations. The mechanism by which the oscillations arise in the delayed system of De Gregorio and colleagues is also shown to be a Hopf bifurcation, like in delayed, spatially dependent EBMs (Bhattacharya et al., 1982;Roques et al., 2014) and in certain simple and intermediate models of the El Niño-Southern Oscillation (ENSO) (Jin et al., 1994;Tziperman et al., 1994;Ghil et al., 2008, and references therein).
This again shows the importance of timescales: only when the effect of temperature on the growth rate of the vegetation is substantially delayed does the system exhibit internal oscillations. The discussion on timescales in Daisyworld is the subject of a recent paper by Weaver and Dyke (2012), where the persistency of regulatory behavior under a separation of timescales is examined. Figure 7 (upper panel) shows plots of temperature (red curve) and vegetation cover (blue curve) versus time. The oscillations have an amplitude of a few degrees for the temperature, but are very large for the vegetation, going from almost no vegetation to 70 % of the land being covered by plants. In addition, the vegetation plot shows a sawtooth-like shape, with the vegetation cover almost vanishing for long times, after which it shoots back up.
Such sawtooth-like behavior characterizes the ice volume evolution during Quaternary glaciations; see, for instance, Ghil and Childress (1987, Ch. 11, and references therein). This behavior has been modeled using relaxation oscillations; cf. Le Treut and Ghil (1983), Ghil (1994) and Crucifix (2012), among others. In these studies, the sawtooth shape is due to a quick melting and slow build-up of the ice sheets, whereas in Gildor and Tziperman (2001), the sea ice plays a key role. In the intermediate-complexity model of Ganopolski and Calov (2011), like in the early conceptual model of Imbrie and Imbrie (1980), it is the strongly nonlinear response of the climate-cryosphere system to orbital forcing. In the present model, it is the surge in vegetation growth accompanying a temperature rise that provides the characteristically asymmetric nature of the oscillation.
Given the importance of sea ice feedback in our model, sea ice extent is plotted in Fig. 7b. In the temperature interval in which the sea ice feedback is active, the sea ice extent depends linearly on T , so the shape of the curve in this panel is similar to that of temperature in the upper panel, with the directions of waxing and waning reversed. There is a slow expansion of sea ice and rapid melting, which resembles the behavior of land ice volume during glacial-interglacial cycles. This highly nonlinear oscillation in our model indicates that it might be of interest to include vegetation in "ice age oscillators", possibly combining it with a simple carbon cycle. As far as we know, this has not been done for ODE models. Certain results of more detailed models, however, do show that vegetation plays an important role in glacial cycles (Meissner et al., 2003;Horton et al., 2010).

Summary
We described a simple dynamical systems model for climate-vegetation interaction. The model planet has a large ocean, which can be covered by sea ice, and a land area, which can be covered by vegetation. The system variables are temperature and vegetation cover, and the coupling between those is given by the growth rate of the latter and by the albedo, both of which are temperature-dependent.
The model is similar to Daisyworld and related models, the main difference being the inclusion of an ocean. Our model is also related to EBMs and, in this respect, the novelty lies in the inclusion of vegetation as a main variable.
The model exhibits two stable states, one with and one without vegetation. The vegetated state can lose its stability through a Hopf bifurcation and give rise to a limit cycle. This happens when the typical timescale for vegetation overturn becomes high enough. The influence of the model's sea ice in reinforcing the albedo feedback is essential to the presence of oscillatory behavior. The oscillations observed are anharmonic and have a sawtooth shape, reminiscent of deglaciations during the Quaternary. We have obtained analytical results on the role of the feedbacks in determining the system's behavior and compared the resulting oscillations with other models.
Although some parameter values in our highly idealized model are not entirely realistic, the results add to the evidence that vegetation, in combination with other feedback effects, can play an important role in affecting climate. The model studied herein is also interesting because it is one of the simplest ODE models for climate-vegetation interactions that exhibits oscillatory behavior. This intrinsic variability is a basic manifestation of how vegetation affects climate, and constitutes an example of how complex behavior arises in the Earth system even at the lowermost levels of the modeling hierarchy. Once demonstrated and understood at such a level, one can ask next whether similar behavior does persist in more detailed models, and whether it does reflect the actual behavior of the natural climate-vegetation system.

Discussion of the results
While the sea ice albedo feedback plays a crucial role in our model's oscillatory behavior, the amplitude of the oscillations in sea ice extent is quite small, from 0 to a mere 6 %. Gildor and Tziperman (2001) first proposed sea ice as a determining factor in the inception of glacial cycles, but our model's oscillations have a much shorter period than the dominant 100 000 year cycle of the latter. Note that our climate-vegetation model differs from that of Gildor and Tziperman (2001) by the absence of the much slower and more massive continental ice sheets, and it differs from the model of Svirezhev and von Bloh (1996) by including the ocean and the sea ice, as well as by the less crucial parametrization of the outgoing radiation. The latter authors actually prove the absence of limit cycles in their climatevegetation model by using the Bendixson-Dulac criterion (e.g., Andronov et al., 1966). This proof reinforces our argument that sea ice plays an import role in intrinsic climate variability on long timescales.
Because of Proposition 1 and the shape of the growth curve, it is essential that the fixed point for which the Hopf bifurcation occurs have a temperature higher than T opt , since β (T ) needs to be negative, and lower than T α,u , so that the sea ice feedback be on. The former condition implies that the optimal growth temperature for plants needs to be relatively low, while the latter condition requires that sea ice be present in polar regions for relatively high global temperatures.
Such distinctions between local and global variables are hard to make in highly idealized, ODE models, and it is clear that there is a need to pursue this type of analysis in more detailed, spatially dependent models, either with several climate zones or with continuous dependence on latitude. Such an extension would allow the model, for example, to have vegetated taiga and tundra interact with sea ice at higher latitudes and desert vegetation interact with atmospheric and oceanic temperatures at lower latitudes.
Spatially dependent versions of Daisyworld have already been studied and yielded interesting results. Adams et al. (2003) observed pattern formation in a spatially onedimensional Daisyworld. In their model, vegetation in different latitude belts interacts with global temperature. Patterns are also seen by Biton and Gildor (2012), whose onedimensional Daisyworld is forced by seasonal variations in solar radiation. Both of these studies suggest that pursuing our model analysis in higher dimension, while still including ocean and sea ice, would be of interest.
Earth system models of intermediate complexity provide some evidence that vegetation plays a role in Quaternary climate variations. In particular, the albedo changes that accompany shifts in vegetation act to amplify orbital forcing (Claussen, 2009). In addition, more detailed models, such as those by Zeng et al. (1999) and Zeng and Neelin (2000) conclude that vegetation in the Sahel region interacts strongly with global temperatures on shorter timescales.
The concept of a hierarchy of models -already fairly well accepted in climate modeling (Schneider and Dickinson, 1974;Ghil and Robertson, 2000;Ghil, 2001) -suggests precisely this kind of further exploration of more detailed models. The present paper clearly shows that intrinsic variability is possible in a climate-vegetation system, at least in its simplest possible form -either as bistability or as selfsustained periodic oscillations -and that the combined effect of vegetation-albedo and ice-albedo feedbacks can give rise to an oscillation-producing mechanism.
The conditions β (T ) < 0 and ∂α/∂T < 0 on these two feedbacks are necessary in our simple model for the vegetated fixed point to become unstable and lead to intrinsic variability. We note that limit cycles are also observed when changing for example the parametrization of outgoing radiation into a modified Stefan-Boltzmann law, or if we change the parameters T α,u and T opt . The important thing is that the conditions above on the derivatives be satisfied.
In our model, the feedback that produces oscillations is negative: increasing temperature lowers vegetation cover, which in turn lowers energy absorption and decreases temperature. Due to the large number of interacting processes in nature, the effective sign of the feedback is not known directly from observations, and it might be different in different areas (Claussen, 2009;Meir et al., 2006). One can think of different models, too, in which another component or physical process than sea ice and the associated feedback provide the requisite oscillatory mechanism.
One such process is the hydrological cycle. Aleina et al. (2013) and Brovkin et al. (1998), for instance, have included it already in conceptual climate-vegetation models but observed no oscillations. A study by Fernando Salazar and Poveda (2009) did observe oscillations in a model where Daisyworld is extended with clouds and a basic hydrological cycle. The authors concluded that interaction between clouds and vegetation can improve the circumstances for life on the planet. Clouds are known to play an important role in changing the radiation balance of the planet, but their exact influence is poorly known (IPCC, 2013).
In addition to the influence of vegetation on the hydrological cycle, there is a whole range of other feedback effects between vegetation and several components of the climate system. Of particular importance here is the biogeochemical feedback, as opposed to the biogeophysical effects, such as albedo change, that we have explored. Biogeochemical effects include the influence of vegetation on the carbon cycle, and therefore on the greenhouse effect and the radiation balance of the Earth. A more detailed model than the one studied here could include a simplified version of these effects. Lovelock and Kump (1994), for instance, have examined the regulatory effect of vegetation, both on land and in the ocean, under different conditions. They conclude that vegetation is involved in regulation for low temperatures, but that the regulatory system collapses at high temperatures.
Another interesting model extension that comes to mind is the inclusion of oceanic vegetation, i.e., phytoplankton. Plankton interacts with climate in different ways, through its possible effect on cloud formation (Ayers and Cainey, 2007), but also through its albedo.

Broader context
The results of this paper have to be viewed in the broader perspective of the hierarchy of climate models already mentioned repeatedly in its preceding sections (Schneider and Dickinson, 1974;Ghil, 1994Ghil, , 2001Ghil, , 2015Dijkstra and Ghil, 2005;Ghil and Robertson, 2000). Recall that this hierarchy ranges from simple, conceptual ODE models (e.g., Stommel, 1961;Lorenz, 1963;Källèn et al., 1979), like the one formulated and analyzed herein -through intermediate models of varying complexity (e.g., Claussen et al., 2002;Ganopolski and Calov, 2011) -all the way up to full-scale GCMs (e.g., IPCC, 2013, and references therein).
Within this hierarchy, the role of the simple models, sometimes referred to as "toy" models, is to provide insight and help understand the behavior of the more complex models, as well as of the climate system itself. The role of the intermediate models is to refine these insights and bridge the gap between the toy models and the GCMs (Ghil, 2001;Claussen et al., 2002): on the one hand, they are still simple enough to allow a fairly thorough analysis of their behavior; on the other, they may be detailed enough for a direct comparison with the GCMs and with increasingly more plentiful and accurate observational data sets.
Finally, GCMs allow an extensive comparison with the observations and can thus help invalidate (Popper, 1959) the theories suggested by toy models, whether "validated" (i.e., not refuted) by intermediate ones or contradicted by the latter. Still, when GCM results are at variance with those of simpler models, it is not always the former that are correct, and it is a careful analysis of the observations that ultimately decides which results are closer to the elusive truth; see the discussion in Dijkstra (2007) and Ghil (2015).
The results of our toy model suggest that vegetation might play a larger role in climatic variability -whether in apparent jumps between two or more types of near-stationary states or in oscillatory behavior -than heretofore suspected. In particular, it might contribute to more-or-less regular, but not necessarily simply periodic variability.
It is clear, for instance, that the clouds' contribution to planetary albedo is larger than that of vegetation (IPCC, 2013). It is also clear that our model's change in land albedo between its bare and vegetated state is unrealistically large; cf. Crucifix and Hewitt (2005). But, the role of clouds has been explored over the last few decades across the entire hierarchy of climate models. While much remains to be done to gain a complete understanding of cloud-radiation interactions, greater attention to the role of vegetation in the evolution of planetary albedo seems to be well worth the effort. In particular, since the timescale of changes in vegetation is considerable slower than that for changes in cloud cover, the former might play a greater role in low-frequency climate oscillations.
Recall also that, in the early days of energy balance models (EBMs; see Ghil, 2001Ghil, , 2015, and references therein), even larger and more unrealistic albedo differences between low-and high-temperature surfaces were used in simple, albeit infinite-dimensional models. Still, the EBMs' suggestion of multiple equilibria being possible in the climate system on long timescales has led to a rich literature on bifurcationsmore recently and excitingly called "tipping points"  -and their potential role in both past and future climate evolution.
While some Earth models of intermediate complexity do indeed show multiple equilibria, these appear to be mostly of local relevance, for instance in the Sahel. Brovkin et al. (2003), though, found no support for the co-existence of multiple equilibria at northern high latitudes. Moreover, when global oscillations do appear in such models, they tend to be attributed to cyclicity in the ocean circulation.
Our paper is only trying to make a case for the possibility of vegetation playing a more important role than contemplated heretofore and does not claim in the least to have definitively proven that this is so. A similar argument about local versus global effects has been made with respect to the oceans' thermohaline circulation. Recall that the Stommel (1961) paper -much quoted recently in the context of multiple equilibria and symmetry breaking in the meridional overturning of the Atlantic or even global ocean -was originally written to explain seasonal changes in the overturning of "large semi-enclosed seas (e.g. Mediterranean and Red Seas)"; see, for instance, Dijkstra and Ghil (2005).
There is no better way of concluding this broader assessment of our toy model's results than by citing Karl Popper: "Science may be described as the art of systematic oversimplification" (Popper, 1982). It might be well to remember this statement, given an increasing tendency in the climate sciences to rely more and more on GCMs, to the detriment of simpler models in the hierarchy.

Appendix A: Proof of Propositions 1 and 2
In this appendix, we prove the propositions from the main text. Recall the definition a = −C −1 T (Q 0 ∂α/∂T + ∂R o /∂T ) , and let b = −C −1 T Q 0 ∂α/∂A. From now on, we focus on the temperature interval where β(T ) is positive and, therefore, where the sea ice feedback is active. For these values of T , both a and b are positive constants, because of the linear dependence of the albedo on T and A and of the outgoing energy R o on T .
This proves Proposition 1, since, when β (T ) > 0, the above inequality is always satisfied. Therefore, det(J) < 0, which implies that the fixed point is a saddle and thus unstable.
We will now use the same notation to prove Proposition 2, which we repeat here for the reader's convenience.
Proposition 2. A model fixed point (T + , A + ) that lies in the region with active sea ice feedback (T α, < T + < T α,u ) is a focus if and only if the following inequality holds: (1 − y) 2 < 4η(y)y, where y = (γ − β(T + ))/a, η(y) is (up to a sign change) the ratio of the derivatives of the two nullclines and a is given by Eq. (A1). Given the inequality (Eq. A4), The focus is stable if y < −1.
Proof. In order to have a focus, the Jacobian matrix J has to have eigenvalues with a non-zero imaginary part. This translates into the condition tr(J) 2 − 4det(J) < 0.
We substitute the values into J above, while using our new constants a and b and the variable x, and get the inequalities where y = x/a and η(y) = β bγ /(β 2 a). To show that the function η actually is, up to a change of sign, the ratio of the slopes of the nullclines, is a matter of computation: a quick calculation shows that the slope of the T -isocline is equal to −a/b, cf. Eq. (7), and that of the A-isocline is γ β (T )/β 2 (T ), cf. Eq. (6). The statement about the stability of the focus follows from the fact that a sufficient condition for the focus to be stable is that tr(J) < 0. This translates into a + x < 0, or y < −1.