On the stability of a climate model for an Earth-like planet with land-ocean coverage

The discovery of potentially habitable Earth-like exoplanets orbiting around hosting stars increased the interest in investigating their climate systems. In this perspective, we study the stability properties of a simple energy balance climate-vegetation model. It describes the interaction between the planetary surface, that can be partially covered by vegetation, a large ocean and a land surface. The model allows to investigate changes in response to variations in the hosting star’s irradiance (e.g. due to mutual distance variations). The stability properties of the system change according to variations of the control parameters, as for example the effective growth rate of vegetation γ and the stellar distance a. In particular, the model presents both stable and unstable fixed points, which can exhibit several kinds of stability properties depending on the distance of the planet from the hosting star.


Introduction
A complete description of the Earth's climate system in terms of the full set of atmospheric variables (e.g., temperature, humidity, pressure) is quite complex, due to their variability on several timescales [1][2][3][4][5][6]. Energy Balance Models (EBMs) have been proposed to investigate the climate system in a global way by modeling the basic atmospheric mechanisms [7][8][9]. The mathematical formulation of EMBs, based on the radiative equilibrium between the incoming solar radiation and the outgoing infrared radiation, is particularly simple since governed by a set of Ordinary Differential equations (ODEs) depending on parameters related to the dynamical properties of the system.
The recent discoveries of Earth-like exoplanets [10][11][12][13] increased the interest of the scientific community towards these astronomical objects, also fostered by the high probability that several more of them will be identified in the near future. Climate models can provide insights on one of the more interesting issues concerning exoplanets: their habitability. The first climate model used to investigate exoplanet climate was those by [14] including a single column atmosphere. More recently, also general circulation models (GCMs) have been used to characterize the possible habitability [15][16][17][18][19][20][21][22][23]. Although containing strong simplifications, EBMs allow to explore a large set of variables inside the parameter space and, in particular, they include the effects of vegetation that distinguishes them from the more complex and complete GCMs. This is particularly important in studying the climate of terrestrial exoplanets since it is possible to take into account the possibility that life affects climate. Moreover, since the knowledge of physical properties of exoplanets is limited, EBMs are suitable to gain new insight into the global properties of their climate [24].
In the framework of the Budyko-Sellers models [7,8], a simple equation governs the time-evolution of the Earth's global mean surface temperature T where C T is the heat capacity of the planet, E in is the incoming solar radiation flux, which depends on the solar constant S and the albedo α(T), and E out is the infrared outgoing planetary radiation flux.
For the Earth-like case, two stable fixed points, corresponding to a glacial and a desert state, and an unstable one, similar to the present climate state, are found [25]. However, it does not take into account some components of the climate system which can affect its evolution, such as the atmosphere-biosphere interaction, a not-constant solar irradiance, and the greenhouse effect. The first attempt to include the climate-vegetation interaction was made in the framework of the well-known Daisyworld model [9]. This model is governed by two ODEs for the time-evolution of two types of vegetation (white and black daisies) coupled with equation (1) for the global surface temperature. It allows to investigate how vegetation regulates the temperature via the albedo feedback when the incoming solar radiation changes (see [26] for a complete review).
A simple climate-vegetation model, based on the ice-albedo feedback characterizing a planet with a large ocean, has been recently investigated in [27]. This model is governed by two nonlinear and coupled ODEs for both temperature and vegetation coverage. The stability properties of this model have been investigated in [27] and three fixed points, two stable and one unstable saddle point, have been found. Results show that different stability states, bistability between a desert and a vegetated state as well as an anharmonic oscillatory behavior manifesting thought a Hopf bifurcation, are present. However, the stability properties can change according to variations of the death rate of vegetation, with the occurrence of Hopf or saddle-node bifurcations. To capture effects of spatial hetereogenity 1-D models have been developed by including latitude dependence and longitudinal diffusion has been introduced [28]. In particular, a constant diffusion term, with a control parameter, was used to investigate the stability of multiple steady states [29]. Recently, in [30] a modified Daisyworld model has been proposed which includes the spatial dependence on the latitude, a non-constant heat diffusion term, that takes into account the differences between the equator and poles, and, for the first time, a parametrization of the greenhouse effect, through a grayness function [8]. They showed that: (i) the system selfregulates even for solar incoming radiation values far from the present Earth's conditions; (ii) the non-constant diffusivity breaks the symmetry of the system with respect to the equator; and (iii) the grayness function contributes to self-regulate the planet's climate, modifying its stability properties, producing a formation of stripped patterns for low luminosity values. Since a direct knowledge of the latitude dependence is not possible for exoplanets the use of 0-D models is preferred. In particular, to make the model in [27] more suitable to describe the global properties of an exoplanet climate we propose a modified version that includes a varying solar incoming radiation flux.
In this paper we present the model and its analytical, in the simple case when the vegetation drops off, and numerical, for the vegetated case, solutions for the stability conditions. We show that the level of incoming radiation, together with the vegetation growth/death rate and the fraction of land covering the planet, are fundamental to set fixed points and the stability properties of the system described by this model.

The model
The Daisyworld-like model considered in this work describes an Earth-like planet covered by land and ocean, with the possible presence of vegetation. It is based on two equations for the global average temperature T and the fraction of land A covered by vegetation (see also [9,[27][28][29][30]) where α(T, A) is the planetary albedo, S is the mean incoming radiation flux, R(T) is the flux outgoing from the planet, β(T) and γ are the vegetation growth and death rates, respectively (see [9,27,30] for more details). To investigate the stability of the model as regards an Earth-like planet climate, the incoming radiation flux is assumed to be S Q a 2 = - , where Q  is the mean solar radiation flux at the Earth's orbit modulated by the parameter a d d =  , where d is the stellar distance and d  is the Sun-Earth distance. The albedo of the planet depends on the fraction of covered land where p is the fraction of land on the planet, while α o , α v and α g represent the albedos for ocean, vegetation and bare ground, respectively, with α v <α g . The albedo of the ocean is not fixed, rather it is allowed to depend on the temperature, to account for the presence of a given fraction of ice. The dependence is given in terms of a ramp function that introduces a nonlinearity in equation (2) leading to multiple equilibria Here α max refers to an ocean completely covered by ice, while α min refers to an ice-free ocean. Generally, the outgoing energy from the planet is described by a black-body radiation process R(T)=σ T 4 , σ being the Stefan-Boltzmann constant. However, a linear functional shape is widely used in literature as alternative to the black-body radiation, according to the Budyko-Sellers framework (see, e.g., [7,8,27,29]). Since we are interested in investigating the linear stability of the model, we use the linear approximation where the phenomenological constants B 0 and B 1 , which are tuned on the optimal temperature T opt , implicitly represent the effects of the present gaseous components of the planetary atmosphere (i.e., their values are set to take into account an Earth-like greenhouse effect). The growth-rate of vegetation is described by a logistic equation depending on the temperature, according to the Daisyworld model, that is, the growth is zero except in an interval of temperatures where it is described by a parabolic function of temperature, with a maximum at the optimal temperature T opt . The death-rate γ is assumed to be constant for simplicity in performing linear stability analysis (as also in [9,27,28,30]), although several models consider a temperature-dependent death-rate (see, e.g., [31]). The parameters usually used in the model to describe the Earth-like climate are reported in table 1. It is worth to remark that in the framework of Daisyworld-like models [9,27,28,30,32], the vegetation (or if preferable the biota) only acts as a feedback in regulating the thermal equilibrium of a planet, without considering any effect due to the carbon cycle or in changing the greenhouse gases dynamics. Moreover, Daisyworld-like models also consider land to be ice-free since no hydrological cycle is considered, although if the land albedo is temperature dependent due to the growth/death of vegetation. This will add complexity to the model, since an additional term, describing land-ocean coupling, should be considered in equation (2). On the other hand, this might quantitatively alter the results by changing the fixed points but any effects should be observed on their stability properties. As also previously reported [9], the mechanisms through which life (e.g., vegetation) affects temperature evolution are not themselves important but it is only required that the biota influence the temperature [32]. In particular, the vegetation-albedo feedback is one of the main mechanisms affecting Daisyworld-like models [9,27,28,30].

Multiple steady-states
The fixed points T  and A  of the system described by equations (2)- (3), satisfy the following equations Table 1. Parameters used in the model (see also table 1 in [27]).
Since β(T) is a piecewise quadratic function, we expect at most three fixed points representing three different climate states of the planet. Looking at the equation for A, we immediately realize that two main classes of fixed points can be supported by the system, namely climate states in the absence or presence of vegetation. In the following we distinguish between these different situations.

Steady-states in the absence of vegetation
In the absence of vegetation, A=0, a fixed point T T  = can be obtained from the classical equation ) , which reduces to the linear equation when we use the linear approximation (7). The equilibrium point then depends on the definition of T o  a ( ) and on the value of the parameters a and p, by keeping fixed the other parameters Q 0 , B 0 , B 1 and T opt . Note that this is a very peculiar condition for EBMs. As an example, let us consider the case of a planet without ocean, say p=1. In this case the equilibrium solution does not depend on the ocean albedo α o , being The simple EBM requires that this equilibrium point is possible (i.e., T 0  > ) for a 'rocky' planet orbiting at a distance a from its hosting star given by The values of B 0 and B 1 are phenomenological constants derived in [7] to fit the departure from a black-body function for the Earth's infrared emission in space.
For an Earth-like planet where p 0 1  < , the equilibrium solution depends on α o and can be easily calculated from equation (10) as In figure 1 we report the contours of the equilibrium temperature in the plane (p, a). An Earth without vegetation lies at the border of the minimum and intermediate albedos, while an Earth with ocean covered by ice can exist even without vegetation.

Steady-states in the presence of vegetation
By solving equations (9a)-(9b) for the vegetated system, the steady-states T A ,   ( ) can be obtained from the pair of equations First of all, it can be noted that when , the only steady-state solution is the absence of vegetation, corresponding to  (5)).
and derived from equation (20). This corresponds to the solution inside the growth range of the vegetation (i.e., T T T , low high  Î ( ]). Conversely, the steady-state solution corresponding to the absence of vegetation outside this range is obtained from equation (21), with a temperature T 242 ~K [27]. On the contrary, a solution with A 0, 1 and T  Î ( T low ,T up ]. In the following we investigate this situation. The solutions of equations (20)-(21) identify two curves, called null isoclines, in the phase plane (T, A) [27], and their intersections define the fixed points of the climate system. The first one, obtained from equation (20), has a parabolic shape, depending on the functional form of T  b ( ), while the second one, derived from equation (21), is a linear function of T  . This suggests that the number of fixed points for T  Î ( T low ,T up ] moves from zero up to two. The null isoclines are reported in figure 2 in the phase space (T, A), for different values of the parameters γ and a and two values of p. The steady states correspond to the intersections of the two curves defined in equations (20)- (21). It can be seen that, in general, low values of the vegetation death-rate γ correspond to two climate states, and one single climate state for distances a<1 with a low equilibrium temperature. As γ increases, some steady states tend to disappear (i.e., for a<1), while two steady states survive when a is close to the Earth's orbit (a∼1−1.2). Higher values of γ correspond to the absence of vegetation and steady states tend to disappear. The same behavior can be observed for different values of p, although with different ranges of equilibrium temperatures and albedos.

Stability analysis of climate states
According to the usual dynamical system theory, the nature and stability of a fixed point can be investigated by studying the eigenvalues λ k (in our case k=2) of the Jacobian matrix J(T, A), its determinant Det(J) and the discriminant Δ(J)=Tr(J) 2 −4Det(J). Table 2 shows the different nature of fixed points, according to the sign of Det(J), Δ(J), and Tr(J).

Stability of the climatic states in the absence of vegetation
In the absence of vegetation A 0  = , the stability of the climate system of the planet can be analytically investigated. In fact, in this case the eigenvalues of the Jacobian matrix are For both extreme cases, namely ice-free and ice covered ocean, the ocean albedo does not depend on temperature, and T 0 . In this case the two eigenvalues λ 1,2 of the matrix are λ 1 =−B 1 /C T and λ 2 =−γ, which correspond to a stable node, whose stability is independent on the various parameters of the model.
In the case when the ocean is partially covered by ice, the albedo depends on temperature and the situation changes. In fact, in this case the node is stable only when which corresponds to a planet not too close to its hosting star with an equilibrium temperature clearly different from T opt , namely which corresponds to In the case of Earth-like parameters, the above conditions result in a<1.25 and T 263, 268 298, 300 ]}, which means that the climate of an Earth-like planet with soil and ocean cannot be stable in the absence of vegetation. When the above conditions are not satisfied, the fixed point can become an unstable node, when both signs of equation (26) are reversed, or a saddle when only one is reversed.

Stability of the vegetated steady-state solution
The stability of the vegetated steady-state solution cannot be analytically investigated and a numerical solution of equation (24) is required. Figure 3 shows the stability diagrams for two different values of p (i.e., p=0.3 and p=0.8 as in figure 2) and three different values of a (a 0.9, 1, 1.1 Î { } ). The upper panels of figure 3 (i.e., (p, a)=(0.3, 0.9) and (p, a)=(0.8, 1)) correspond to the case in which only one fixed point characterizes the vegetated state, while the lower panels of figure 3 (i.e., (p, a)=(0.3, 1) and (p, a)=(0.8, 1.1)) refer to the situation in which two stable fixed points are found for the vegetated state (see also figure 2). As pointed out before, several steady states can be identified, according to the different nature of the eigenvalues.  Node ) asymptotically to the trajectory associated with the negative eigenvalue, then turn and flow out to infinity asymptotically to the trajectory associated with the positive eigenvalue.
Conversely, when a 0.98, 1.05 Î [ ], two fixed points can be found from which a bifurcation diagram is obtained (see figure 3 for p=0.3, a=1). As it is possible to see, a saddle-node bifurcation occurs for γ∼0.41 from which stable and unstable solutions branch off. Particularly, the fixed point characterized by a lower temperature (see figure 2) can be classified as an unstable saddle. On the other hand, the fixed point associated with higher temperature is a stable node, with the Jacobian matrix T A J , node node   ( ) which has two distinct real negative eigenvalues ( R , 1 2 l l Î -). This suggests that all phase-space trajectories flow into the stable fixed point ), initially asymptotically to the trajectory associated to the smaller eigenvalue, then turn and flow asymptotically to the trajectory associated to the larger eigenvalue.
When p increases (i.e., right panel in figure 3 with p=0.8), some differences/similarities are found in the nature of the steady-state points obtained with different values of a. In this case we have that: (i) no fixed points are found when a<0.96 and a>1.25, (ii) only one fixed point is found when a ä [0.96,1.05), and (iii) two fixed points are found when a ä [1.05, 1.25]. When only one fixed point is found, it is always a stable fixed point, classified as a focus for γ<0.34 and as a node when γ>0.34. This is different with respect to the previous case (i.e., for p=0.3), when an unstable saddle point was found, due to the fact that the slope of the nullcline obtained from equation (21) changes in sign with p. This changes the positive eigenvalue to a negative one, leading a stable fixed point whose nature depends on γ. In this case, the only found fixed point has a temperature in the range T 290, 300  Î ( ) K, similar to the larger temperature fixed point obtained for p=0.3 and a ä [0.98,1.05]. This focus-node bifurcation occurs when a hyperbolic equilibrium is found. It is characterized by eigenvalues having real parts of the same sign, moving from real to complex values (i.e., the discriminant changes its sign from positive to negative). This suggests that both p and a parameters change the nature of the fixed points, with the concurrent effect of the death-rate γ.
Conversely, similar results to the Earth-like land distribution (p=0.3) are found when two fixed points are obtained. Indeed, also in this case, they can be classified as an unstable saddle, for the fixed point with smaller temperature T A , saddle saddle   ( ), and a stable node when temperature is larger T A , node node   ( ). Figure 4 shows the phase space nullclines and invariant manifolds, for the case a=1 and γ=0.2, in which the manifolds of the unstable saddle point are reported T A , 282.9, 0.78   = ( ) ( ). Indeed, the vertical black line is the stable manifold, while the unstable one is represented by the black curve characterized by a parabolic-like shape. This behavior is related to the ice-albedo feedback which changes the stability properties of the unstable fixed point which can move along the unstable manifold. system to orbital forcing [42] or to internal 'seesaw' mechanisms [43]. [27] also showed that, since the shapes of the temperature and vegetations curves are similar but time-shifted, this could resemble the behavior of land ice volume during glacial-interglacial cycles. For these reasons, the vegetation dynamics needs to be included in iceoscillatory models to describe the observed nonlinear oscillations.

Conclusions
In this paper we investigated the stability properties of a simple energy balance climate-vegetation model, consisting of a planetary surface partially covered by vegetation, a large ocean and a land surface, which can change in response to variations in stellar irradiance as, for example, related to its variable distance with respect to the hosting star. We showed that the nature of the fixed points and their stability change according to variations of the control parameters, say the effective growth rate of vegetation γ and the stellar irradiance a. We note that the climate of an Earth-like planet (e.g., with similar land-ocean areal proportions and similar stellar distance as the Earth) with soil and ocean cannot be stable in the absence of vegetation, showing unstable fixed points referring to unstable nodes or saddles. Conversely, in the presence of vegetation the Earth-like case (p=0.3, being p the fraction of land covering the planetary surface) shows: (i) unstable saddle points when a ä [0.87, 1.05], and (ii) stable nodes for a ä [0.98, 1.05) (being a the ration between the stellar distance and the Sun-Earth distance). When p increases (as for example when p=0.8 in figure 3 left panels) a stable fixed point, classified as a focus for γ<0.34 and as a node when γ>0.32, is found (being γ the death-rate of vegetation). Moreover, when from equations (20)-(21) two fixed points are obtained, they can be classified as an unstable saddle and a stable node.
Finally, we also show that, since the solar irradiance is not constant but depends on the control parameter a, we can have a Hopf bifurcation for different combinations of the γ and a parameters (i.e., p=0.3, a ä [0.98,1.02], and γ ä(0.01,0.03)). The observed oscillatory behavior needs to be investigated with more accuracy because it can describe quasi-periodic behaviors which were observed in climatic time series on Earth and could occur also on other Earth-like planets.
Including the effect of variable stellar irradiance on the planet's climate represents a significant addition to the model proposed in [27] and, therefore, it's important to understand the stability properties of the model with respect to its control parameters. Notwithstanding its simplicity, the EBM model proposed here represents, indeed, a useful tool to study the climate of Earth-like exoplanets. Of course, the model has some substantial limitations with respect to more detailed and complete models which include, for example, a description of the planet's atmosphere, such as GCM models [15][16][17][18][19][20][21][22][23]. On the other hand, taking into account that the physical information available on exoplanetary systems is typically limited, the present model represent a useful complementary tool for the study of exoplanets, as it allows to investigate in a simple way a large set of possible cases obtained for different irradiance values, ocean-land fractions, and vegetation coverages. This is also confirmed by the recent application in [24] of the model to the study of the planets of the TRAPPIST-1 system.