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

T Alberti, F Lepreti, A. Vecchio, V. Carbone

► To cite this version:

T Alberti, F Lepreti, A. Vecchio, V. Carbone. On the stability of a climate model for an Earth-like planet with land-ocean coverage. Journal of Physics Communications, 2018, 2 (6), pp.065018. 10.1088/2399-6528/aacd8d. hal-02304635

HAL Id: hal-02304635

https://hal.science/hal-02304635

Submitted on 12 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

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

To cite this article: T Alberti et al 2018 J. Phys. Commun. 2 065018

View the article online for updates and enhancements.

## You may also like

- Organosilicon-Based Electrolytes with Superior Thermal and Electrochemical Stability to Enable High Energy Lithium Ion Batteries

Monica Lee Usrey, Adrian Pena Hueso, Michael Pollina et al.

- Caveat of High Temperature Accelerated Stability Test for Anion Exchange Membranes

Eun Joo Park, Sandip Maurya, Michael Hibbs et al.

- Summary of research on stability of independent microgrid based on high permeability distributed generation

Qun Wang, Feng Li, Baohua Sun et al.

## Recent citations

- Effect of Vegetation on the Temperatures of TRAPPIST-1 Planets

Antonio Vecchio et al.
Journal of Physics Communications

<figure><img src="image_2.png" /></figure>

CrossMark

PAPER

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

T Alberti<sup>1,2</sup> ●, F Lepreti<sup>2</sup>, A Vecchio<sup>3,4</sup> and V Carbone<sup>2</sup>

¹ INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere 100, 00133 Roma, Italy

² Dipartimento di Fisica, Università della Calabria, Ponte P. Bucci Cubo 31C, 87036 Rende (CS), Italy

³ LESIA-Observatoire de Paris, PSL Research University, 5 place Jules Janssen, 92190 Meudon, France

⁴ Department of Astrophysics/IMAPP - Radboud University, P.O. Box 9010, 6500GL Nijmegen, The Netherlands

E-mail: tommaso.alberti@iaps.inaf.it

Keywords: climate, energy-balance model, fixed points and stability, daisyworld-like model

OPEN ACCESS

RECEIVED
16 January 2018

REVISED
6 June 2018

ACCEPTED FOR PUBLICATION
19 June 2018

PUBLISHED
29 June 2018

Original 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.

<figure><img src="image_3.png" /></figure>

## Abstract

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 $\gamma$ 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.

## 1. 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–6]. Energy Balance Models (EBMs) have been proposed to investigate the climate system in a global way by modeling the basic atmospheric mechanisms [7–9]. The mathematical formulation of EBMs, 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–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–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$

$$C_T \frac{dT}{dt} = E_{in}(S, \alpha) - E_{out}(T), \quad (1)$$
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 $\alpha(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 heterogeneity 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 self-regulates 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.

## 2. 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–30])

$$C_T \frac{dT}{dt} = [1 - \alpha(T, A)]S - R(T), \qquad (2)$$

$$\frac{dA}{dt} = A[\beta(T)(1 - A) - \gamma], \qquad (3)$$

where $\alpha(T, A)$ is the planetary albedo, $S$ is the mean incoming radiation flux, $R(T)$ is the flux outgoing from the planet, $\beta(T)$ and $\gamma$ 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_{\odot}a^{-2}$, where $Q_{\odot}$ is the mean solar radiation flux at the Earth's orbit modulated by the parameter $a = d/d_{\odot}$, where $d$ is the stellar distance and $d_{\odot}$ is the Sun-Earth distance. The albedo of the planet depends on the fraction of covered land

$$\alpha(T, A) = (1 - p)\alpha_o(T) + p[\alpha_v A + \alpha_g(1 - A)], \qquad (4)$$

where $p$ is the fraction of land on the planet, while $\alpha_o$, $\alpha_v$ and $\alpha_g$ represent the albedos for ocean, vegetation and bare ground, respectively, with $\alpha_v < \alpha_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
<table><caption>**Table 1.** Parameters used in the model (see also table 1 in [27]).</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>C<sub>T</sub></td><td>500</td><td>W yr K<sup>-1</sup> m<sup>-2</sup></td></tr><tr><td>Q<sub>☉</sub></td><td>342.5</td><td>W m<sup>-2</sup></td></tr><tr><td>p</td><td>0.3</td><td></td></tr><tr><td>α<sub>v</sub></td><td>0.1</td><td></td></tr><tr><td>α<sub>g</sub></td><td>0.4</td><td></td></tr><tr><td>α<sub>max</sub></td><td>0.85</td><td></td></tr><tr><td>α<sub>min</sub></td><td>0.25</td><td></td></tr><tr><td>T<sub>low</sub></td><td>263</td><td>K</td></tr><tr><td>T<sub>high</sub></td><td>300</td><td>K</td></tr><tr><td>T<sub>opt</sub></td><td>283</td><td>K</td></tr><tr><td>B<sub>0</sub></td><td>200</td><td>W m<sup>-2</sup></td></tr><tr><td>B<sub>1</sub></td><td>2.5</td><td>W K<sup>-1</sup> m<sup>-2</sup></td></tr><tr><td>k</td><td>0.004</td><td>yr<sup>-1</sup> K<sup>-2</sup></td></tr><tr><td>γ</td><td>0.2</td><td>yr<sup>-1</sup></td></tr></tbody></table>

$$ \alpha_o(T) = \begin{cases} \alpha_{\max} & \text{if } T \le T_{\text{low}} \\ \alpha_{\max} + \xi(T - T_{\text{low}}) & \text{if } T_{\text{low}} < T \le T_{\text{high}} \\ \alpha_{\min} & \text{if } T_{\text{high}} < T \end{cases} \quad (5) $$

where

$$ \xi = \frac{\alpha_{\min} - \alpha_{\max}}{T_{\text{high}} - T_{\text{low}}}. \quad (6) $$

Here $\alpha_{\max}$ refers to an ocean completely covered by ice, while $\alpha_{\min}$ refers to an ice-free ocean.

Generally, the outgoing energy from the planet is described by a black-body radiation process $R(T) = \sigma T^4$, $\sigma$ 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

$$ R(T) = B_0 + B_1(T - T_{\text{opt}}), \quad (7) $$

where the phenomenological constants $B_0$ and $B_1$, which are tuned on the optimal temperature $T_{\text{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,

$$ \beta(T) = \max\{0; 1 - k(T - T_{\text{opt}})^2\}, \quad (8) $$

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_{\text{opt}}$. The death-rate $\gamma$ 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].

### 3. Multiple steady-states

The fixed points $T_*$ and $A_*$ of the system described by equations (2)–(3), satisfy the following equations

$$ 0 = [1 - \alpha(T_*, A_*)]S - R(T_*), \quad (9a) $$
$$0 = A_*[\beta(T_*) (1 - A_*) - \gamma]. \qquad (9b)$$

Since $\beta(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.

### 3.1. 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 $R(T_*) = (1 - \alpha)S$, which reduces to the linear equation

$$T_* + \frac{S}{B_1}(1 - p)\alpha_o(T_*) = T_{opt} + \frac{S - B_0}{B_1} - p\alpha_g\frac{S}{B_1}, \qquad (10)$$

when we use the linear approximation (7). The equilibrium point then depends on the definition of $\alpha_o(T_*)$ 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 $\alpha_o$, being

$$T_* = T_{opt} + \frac{S}{B_1}(1 - \alpha_g) - \frac{B_0}{B_1}. \qquad (11)$$

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

$$a^{-2} > \frac{Q_{\odot}^{-1}}{(1 - \alpha_g)}(B_0 - T_{opt} B_1), \qquad (12)$$

a condition that is only satisfied when $(B_0/B_1) > T_{opt}$. 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 $0 \le p < 1$, the equilibrium solution depends on $\alpha_o$ and can be easily calculated from equation (10) as

$$T_* = \left( T_{opt} - \frac{B_0}{B_1} \right) + \frac{S}{B_1} \left[ 1 - \alpha_o(T_*) - p(\alpha_g - \alpha_o(T_*)) \right]. \qquad (13)$$

These solutions exist only for determined distances of the planet from the hosting star. The ice-free ocean planet (i.e., $\alpha_o = \alpha_{min}$ and $T_* \ge T_{high}$) exists only when the distance from the hosting star is small enough

$$a^{-2} \ge \frac{Q_{\odot}^{-1}[(T_{high} - T_{opt})B_1 + B_0]}{[1 - \alpha_{min} - p(\alpha_g - \alpha_{min})]}. \qquad (14)$$

In the opposite case, an ice-covered planet (i.e., $\alpha_o = \alpha_{max}$ and $T_* \le T_{low}$) can exist only when

$$a^{-2} \le \frac{Q_{\odot}^{-1}[(T_{low} - T_{opt})B_1 + B_0]}{[1 - \alpha_{max} - p(\alpha_g - \alpha_{max})]}. \qquad (15)$$

When only a fraction of ocean is covered by ice, say when $T_* \in [T_{low}, T_{high}]$, the fixed points are given by

$$T_* = \frac{(T_{opt} B_1 - B_0) + S[1 - p\alpha_g - (1 - p)(\alpha_{max} - \xi T_{low})]}{B_1 + S(1 - p)\xi}, \qquad (16)$$

which is satisfied when

$$a^{-2} \in [f(T_{low}), f(T_{high})], \qquad (17)$$

where

$$f(x) = \frac{Q_{\odot}^{-1}}{\Pi} [B_1(T_{opt} - x) - B_0], \qquad (18)$$

and

$$\Pi = T_{low}(1 - p)\xi - [1 - p\alpha_g - (1 - p)(\alpha_{max} - \xi T_{low})]. \qquad (19)$$

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.
<figure><img src="image_4.png" /><figcaption>Figure 1. Contours of the equilibrium solutions $T_*$ in the plane $(p, a)$, when the vegetation is set to zero, namely $A_* = 0$. Panel (a) refers to the maximum albedo for an ocean covered by ice, $\alpha_o = \alpha_{max}$; panel (b) refers to the minimum albedo $\alpha_o = \alpha_{min}$ for an ice-free ocean; panel (c) refers to the intermediate case when only a fraction of the ocean is covered by ice (depending on temperature as in equation (5)).</figcaption></figure>

## 3.2. 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

$$A_* = 1 - \frac{\gamma}{\beta(T_*)}, \qquad (20)$$

$$A_* = \frac{\mu}{p} \left\{ \left[ 1 - \frac{R(T_*)}{Q_\odot a^{-2}} - (1-p)\alpha_o(T_*) \right] - p\alpha_g \right\}, \qquad (21)$$

where $\mu = 1/(\alpha_v - \alpha_g)$.

First of all, it can be noted that when $\beta(T_*) = \gamma$, the only steady-state solution is the absence of vegetation, corresponding to
$$A_{\star} = 0, \qquad (22)$$

$$T_{\star} = T_{opt} \pm \sqrt{\frac{1-\gamma}{k}}, \qquad (23)$$

and derived from equation (20). This corresponds to the solution inside the growth range of the vegetation (i.e., $T_{\star} \in (T_{low}, T_{high}$]). Conversely, the steady-state solution corresponding to the absence of vegetation outside this range is obtained from equation (21), with a temperature $T_{\star} \sim 242$ K [27].

On the contrary, a solution with $A_{\star} \in (0, 1]$ corresponds to $\beta(T_{\star}) > \gamma$ and $T_{\star} \in (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 $\beta(T_{\star})$, while the second one, derived from equation (21), is a linear function of $T_{\star}$. This suggests that the number of fixed points for $T_{\star} \in (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 $\gamma$ 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 $\gamma$ correspond to two climate states, and one single climate state for distances $a < 1$ with a low equilibrium temperature. As $\gamma$ 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 \sim 1 - 1.2$). Higher values of $\gamma$ 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.

## 4. 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 $\lambda_k$ (in our case $k = 2$) of the Jacobian matrix $J(T, A)$,

$$J = \begin{bmatrix} -\frac{1}{C_T} \left( S \frac{\partial \alpha}{\partial T} + \frac{\partial R}{\partial T} \right) & -\frac{1}{C_T} \left( S \frac{\partial \alpha}{\partial A} \right) \\ \frac{\partial \beta(T)}{\partial T} A(1-A) & \beta(T)(1-2A) - \gamma \end{bmatrix} \qquad (24)$$

its determinant $Det(J)$ and the discriminant $\Delta(J) = Tr(J)^2 - 4Det(J)$. Table 2 shows the different nature of fixed points, according to the sign of $Det(J)$, $\Delta(J)$, and $Tr(J)$.

### 4.1. Stability of the climatic states in the absence of vegetation

In the absence of vegetation $A_{\star} = 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

$$\lambda_1 = -\frac{1}{C_T} \left[ S(1-p) \left( \frac{\partial \alpha_o}{\partial T} \right)_{T=T_{\star}} + B_1 \right], \\ \lambda_2 = \beta(T_{\star}) - \gamma. \qquad (25a)$$

For both extreme cases, namely ice-free and ice covered ocean, the ocean albedo does not depend on temperature, and $\beta(T_{\star}) = 0$ because either $T_{\star} > T_{high}$ or $T_{\star} < T_{low}$. In this case the two eigenvalues $\lambda_{1,2}$ of the matrix are $\lambda_1 = -B_1/C_T$ and $\lambda_2 = -\gamma$, 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

$$S(1 - p)\xi + B_1 > 0, \\ \beta(T_{\star}) - \gamma < 0, \qquad (26)$$

which corresponds to a planet not too close to its hosting star

$$a^2 < (1 - p) \left( \frac{Q_{\odot}}{B_1} \right) \left( \frac{\alpha_{\max} - \alpha_{\min}}{T_{high} - T_{low}} \right), \qquad (27)$$
<figure><img src="image_5.png" /><figcaption>Figure 2. The (T, A) phase plane with the null isoclines obtained from equation (20) (black curve) and equation (21) (colored lines), respectively. The fixed points are the intersections of the black curves with the colored ones. The colored lines refer to three different values of the parameter a. Left and right panels are obtained by fixing p = 0.3 and p = 0.8, respectively, and four different values of γ (γ = 0.02, 0.1, 0.2, 0.4).</figcaption></figure>

with an equilibrium temperature clearly different from $T_{opt}$, namely

$$ (T_{*} - T_{opt})^{2} > \left( \frac{1 - \gamma}{k} \right), \qquad (28) $$
<table><caption>**Table 2.** Nature and stability of the fixed points.</caption><thead><tr><th>Fixed point</th><th>Det(J)</th><th>Δ(J)</th><th>Tr(J)</th></tr></thead><tbody><tr><td>Node</td><td>+</td><td>+</td><td>+ (unstable), - (stable)</td></tr><tr><td>Focus</td><td>+</td><td>-</td><td>+ (unstable), - (stable)</td></tr><tr><td>Saddle</td><td>-</td><td>+</td><td>(unstable)</td></tr><tr><td>Center (marginal stability)</td><td>+</td><td>-</td><td>=0</td></tr></tbody></table>

<figure><img src="image_6.png" /></figure>

which corresponds to

$$ 
\begin{aligned}
T_{\text{low}} < T_{\star} < T_{\text{opt}} & - \sqrt{\frac{1-\gamma}{k}}, \\
T_{\text{opt}} + \sqrt{\frac{1-\gamma}{k}} < T_{\star} \le T_{\text{high}}.
\end{aligned} 
\qquad (29) 
$$

In the case of Earth-like parameters, the above conditions result in $a < 1.25$ and $T_{\star} \in \{[263, 268] \cup [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.

## 4.2. 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 \in \{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.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Phase space nullclines and invariant manifolds for the vegetated-state with $\gamma = 0.2$ and $a = 1$ to consider the Earth-like case. The stable manifold is represented by the vertical line, while the unstable one is recovered through the parabolic-like function. The diamond refers to the unstable fixed point corresponding to ($T_*$, $A_*$) = (282.9, 0.78). Blue and red lines refer to the nullclines (see equations (20)–(21)).</figcaption></figure>

For the vegetated state in the Earth-like case (i.e., $p = 0.3$) we found that: (i) no fixed points are found when $a < 0.87$ and $a > 1.05$, (ii) only one fixed point is found when $a \in [0.87, 0.98]$, and (iii) two fixed points are found when $a \in [0.98, 1.05]$.

From figure 3 we note that, when only one fixed point is found for the vegetated steady-state, it is an unstable fixed point for each value of $\gamma$. More specifically, for $a \in [0.87, 0.98]$, the vegetated steady-state can be classified as an unstable saddle point. This means that the Jacobian matrix $J(T_{*}^{\text{saddle}}, A_{*}^{\text{saddle}})$ admits two distinct real eigenvalues ($\lambda_1, \lambda_2 \in \mathbf{R}$), opposite in sign, such that all trajectories in the phase-space initially flow towards the fixed point ($T_{*}^{\text{saddle}}, A_{*}^{\text{saddle}}$) 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 \in [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 $\gamma \sim 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 $J(T_{*}^{\text{node}}, A_{*}^{\text{node}})$ which has two distinct real negative eigenvalues ($\lambda_1, \lambda_2 \in \mathbf{R}^{-}$). This suggests that all phase-space trajectories flow into the stable fixed point ($T_{*}^{\text{node}}, A_{*}^{\text{node}}$), 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 \in [0.96, 1.05]$, and (iii) two fixed points are found when $a \in [1.05, 1.25]$. When only one fixed point is found, it is always a stable fixed point, classified as a focus for $\gamma < 0.34$ and as a node when $\gamma > 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 $\gamma$. In this case, the only found fixed point has a temperature in the range $T_* \in (290, 300)$ K, similar to the larger temperature fixed point obtained for $p = 0.3$ and $a \in [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 $\gamma$.

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_{*}^{\text{saddle}}, A_{*}^{\text{saddle}}$), and a stable node when temperature is larger ($T_{*}^{\text{node}}, A_{*}^{\text{node}}$). Figure 4 shows the phase space nullclines and invariant manifolds, for the case $a = 1$ and $\gamma = 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.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (Upper panels) Numerical results showing both an oscillatory and a damped-oscillatory behavior corresponding to the limit cycle (left panel) and to the supercritical Hopf bifurcation (right panel) when the system moves towards a stable focus. (Lower panels) Invariant manifolds for the two cases associated to the normal and supercritical Hopf bifurcations.</figcaption></figure>

An interesting situation can be investigated by looking at the fixed points where $\Delta(J) = 0$, indicating that a Hopf bifurcation occurs. This curve corresponds to the phase-plane points where a fixed point loses its stability and a periodic solution arises, with a pair of complex conjugate eigenvalues that cross the complex plane's imaginary axis. Interestingly, since we inserted a control parameter for solar irradiance (i.e., $a$) we found that a Hopf bifurcation occurs when $p = 0.3$, $a \in [0.98,1.02]$, and $\gamma \in (0.01,0.03)$; however, by fixing $a = 1$, a Hopf bifurcation occurs when $\gamma \approx 0.02$ as reported in previous studies [27]. Figure 5 shows both an oscillatory and a damped-oscillatory behavior observed through numerical solutions of equations (2)–(3) by setting $a = 1$ and for two different values of $p = 0.3, 0.8$. The first case corresponds to the Hopf bifurcation when amplitude temperature oscillations of $\Delta T \sim 2 - 3$ K on $\sim 800$-yr timescale are found. This behavior is in good agreement with that observed in [27] who found temperature oscillations with an amplitude of 2 K and larger amplitude oscillations for the vegetation (from no vegetation to about 70% of covered land). These results were also pointed out in several works [27, 33, 34] where the Daisyworld system exhibits self-sustained oscillations in a range of $C_T$ and $\gamma$. [30] also showed that, when small temperature perturbations are introduced in the initial condition, a striped pattern occurs in the range of low-to-intermediate luminosity values; conversely, for high luminosities (in this case $a < 1$), the greenhouse effect leads to the disappearance of the striped patterns. The main difference is that the albedo feedback of the Daisyworld model is due to two types of daisies, black and white, respectively. This bimodal behavior of the albedo is somehow similar to the role of sea ice in the model developed in [27].

When $p$ is set to 0.8 (right panels in figure 5), a different dynamical behavior is found, with oscillations that are damped, reaching a stable state (without oscillations) corresponding to $T_\star = 298.6$ K. From a dynamical system point of view, this corresponds to the case of a supercritical Hopf bifurcation for which the system approaches a stable (stationary) point which can be classified as a focus [35]. The main difference with the previous case, when $p = 0.3$, is that, although in both cases the fixed point is stable, when a Hopf bifurcation occurs it is marginally stable, meaning that the system oscillates around it, while, when a supercritical Hopf bifurcation is found the system tends to move towards a stationary solution (e.g., a focus point).

The observed oscillatory behavior needs to be investigated with more accuracy because it can reproduce several quasi-periodic behaviors observed in climatic time series [5, 6, 36–41]. Indeed, the climate system exhibits periodic or quasi-periodic oscillations, observed by different climatic proxies (e.g., $\delta^{18}$O, carbon dioxide), at different timescales (from years to millennia). These oscillations are related to several physical mechanisms of both internal and external origin [6] and present both sawtooth-like and smoothed behaviors, characterizing changes in the ice sheets' volume. This behavior can be seen as a nonlinear response of the climate
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 ice-oscillatory models to describe the observed nonlinear oscillations.

## 5. 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 $\gamma$ 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 \in [0.87, 1.05]$, and (ii) stable nodes for $a \in [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 $\gamma < 0.34$ and as a node when $\gamma > 0.32$, is found (being $\gamma$ 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 $\gamma$ and $a$ parameters (i.e., $p = 0.3, a \in [0.98,1.02]$, and $\gamma \in (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–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.

## Acknowledgments

We thank the anonymous reviewers for fruitful and helpful suggestions.

## ORCID iDs

T Alberti https://orcid.org/0000-0001-6096-0220

## References

[1] Tessier Y, Lovejoy S and Schertzer D 1993 *J. Appl. Meteorology* 32 223

[2] Schmitt F, Lovejoy S and Schertzer D 1995 *Geophys. Res. Lett.* 22 13

[3] Marsh N D and Ditlevsen P D 1997 *J. Geophys. Res.* 102 11219

[4] Lovejoy S and Schertzer D 2013 *The Weather and Climate: Emergent Laws and Multifractal Cascades* (Cambridge: Cambridge Univ. Press)

[5] Alberti T, Lepreti F, Vecchio A, Bevacqua E, Capparelli V and Carbone V 2014 *Clim. Past* 10 1751

[6] Shao Z G and Ditlevsen P D 2016 *Nature Commun.* 7 10951

[7] Budyko M I 1969 *Tellus* 21 611

[8] Sellers W D 1969 *J. Appl. Meteorol.* 8 392

[9] Watson A J and Lovelock J E 1983 *Tellus* B 35 284

[10] Beaulieu J-P et al 2006 *Nature* 439 437

[11] Anglada-Escudé G et al 2016 *Nature* 536 437

[12] Gillon M et al 2016 *Nature* 533 221
[13] Gillon M et al 2017 *Nature* 542 456

[14] Kasting J F, Whitmire D P and Reynolds R T 1993 *Icarus* 101 108

[15] Spiegel D S, Menou K and Scharf C A 2008 *Astrophys. J.* **681** 1609

[16] Merlis T M and Schneider T 2010 *J. Adv. Model. Earth Syst.* 2 13

[17] Pierrehumbert R T 2011 *Astrophys. J.* 726 L8

[18] Heng K, Menou K and Phillipps P J 2011 *Mon. Not. R. Astron. Soc.* 413 2380

[19] Heng K, Frierson D M W and Phillipps P J 2011 *Mon. Not. R. Astron. Soc.* 418 2669

[20] Menou K 2012 *Astrophys. J. Lett.* 744 L16

[21] Showman A P, Wordsworth R D, Merlis T M and Kaspi Y 2013 *Comparative Climatology of Terrestrial Planets* (Tuscon, AZ: Univ. Arizona Press) pp 277–326

[22] Wolf E T 2017 *Astrophys. J. Lett.* 839 L1

[23] Turbet M, Bolmont E, Leconte J, Forget F, Selsis F, Tobie G, Caldas A, Naar J and Gillon M 2017 *Astron. Astrophys.* 612 A86

[24] Alberti T, Carbone V, Lepreti F and Vecchio A 2017 *Astrophys. J.* 844 19

[25] Ghil M 1976 *J. Atmos. Sci.* 33 3

[26] Wood A J, Ackland G J, Dyke J G, Williams H T P and Lenton T M 2008 *Rev. Geophys.* 46 RG1001

[27] Rombouts J and Ghil M 2015 *Nonlin. Proc. Geophys.* 22 275

[28] Adams B, Carr J, Lenton T M and White A 2003 *J. Theor. Biol.* 223 505

[29] Adams B and Carr J 2003 *Nonlinearity* 16 1339

[30] Alberti T, Primavera L, Vecchio A, Lepreti F and Carbone V 2015 *Phys. Rev. E* 92 052717

[31] Lenton T M and Lovelock J E 2001 *Tellus* B 53 288

[32] Lovelock J E and Watson A J 1982 *Planet. Space Sci.* 30 793

[33] De Gregorio S, Pielke R A and Dalu G A 1992 *J. Nonlin. Sci.* 2 293

[34] Nevison C, Gupta V and Klinger L 1999 *Tellus* B 51 806

[35] Benhabib J and Nishimura K 2012 The hopf bifurcation and existence and stability of closed orbits in multisector models of optimal economic growth *Nonlinear Dynamics in Equilibrium Models* (Berlin, Heidelberg: Springer)

[36] Denton G H and Karlén W 1973 *Quaternary Res.* 3 155

[37] Ditlevsen P D, Kristensen M S and Andersen K K 2005 *J. Climate* 18 2594

[38] Ditlevsen P D, Andersen K K and Svensson A 2007 *Clim. Past.* 3 129

[39] Rahmstorf S 1994 *Nature* 372 82

[40] Rahmstorf S 2003 *Geophys. Res. Lett.* 30 1510

[41] Schulz M 2002 *Paleoceanography* 17 4

[42] Ganopolski A and Calov R 2011 *Clim. Past.* 7 1415

[43] Stocker T F 1998 *Science* 282 61