Stability of a planetary climate system with the biosphere competing for resources

TL;DR

This paper analyzes the stability of planetary climate systems influenced by biosphere feedbacks, demonstrating that high biodiversity maintains stability, while reduced biodiversity can lead to climate change, using a coupled climate-biosphere model.

Contribution

It provides a mathematical framework linking climate stability to biosphere biodiversity, highlighting the importance of biodiversity in maintaining planetary habitability.

Findings

High biodiversity sustains climate stability.

Reduced biodiversity can cause climate instability.

The model exhibits slow-fast dynamical behavior.

Abstract

With the growing number of discovered exoplanets, the Gaia concept finds its second wind. The Gaia concept defines that the biosphere of an inhabited planet regulates a planetary climate through feedback loops such that the planet remains habitable. Crunching 'Gaia' puzzle has been a focus of intense empirical research. Much less attention has been paid to the mathematical realization of this concept. In this paper, we consider the stability of a planetary climate system with the dynamic biosphere by linking a conceptual climate model to a generic population dynamics model with random parameters. We first show that the dynamics of the corresponding coupled system possesses multiple timescales and hence falls into the class of slow-fast dynamics. We then investigate the properties of a general dynamical system to which our model belongs and prove that the feedbacks from the biosphere dynamics cannot break the system's stability as long as the biodiversity is sufficiently high. That may explain why the climate is apparently stable over long time intervals. Interestingly, our coupled climate-biosphere system can lose its stability if biodiversity decreases; in this case, the evolution of the biosphere under the effect of random factors can lead to a global climate change.

Figures

Equations

$$\frac{dT}{dt}={\lambda }^{-1}\left(-e\sigma T^4+\frac{{\mu }_0{I}_0}{4}(1-A)\right),$$

$$\frac{d{x}_i}{dt}=x_i(-{\mu }_i+{\varphi }_i(v)-{\gamma }_i{x}_i),i=1,\dots ,m,$$

$$\frac{d{v}_k}{dt}=D_k(S_k-v_k)-\sum _{i=1}^M{b}_{ki}{x}_i{\varphi }_i(v),k=1,\dots ,n,$$

$$0\le {\varphi }_j(v)\le C_{+},|{\varphi }_j(v)-{\varphi }_j(\tilde{v})|\le L_j|v-\tilde{v}|,$$

$${\varphi }_k(v)=0,forallk,v\in \partial {\mathbf{R}}_{+}^m$$

$$\frac{\partial {\varphi }_k(v)}{\partial v_j}\ge 0,forallk,j,v\in \partial {\mathbf{R}}_{+}^M.$$

$$\phi_{k}(v)=r_{k}\min\Big{\{}\frac{v_{1}}{K_{k1}+v_{1}},...,\frac{v_{m}}{K_{km}+v_{m}}\Big{\}}$$

$$A=A(x)=A_0+m^{-1}\sum _{j=1}^m{c}_j{x}_j.$$

$$\frac{d{x}_i}{dt}=x_i(-{\mu }_i+{\varphi }_i(v)-{\gamma }_i{x}_i),i=1,\dots ,m,$$

$$\frac{d{v}_k}{dt}=D_k(S_k(T)-v_k)-\sum _{i=1}^m{b}_{ki}{x}_i{\varphi }_i(v),k=1,\dots ,n.$$

$$\frac{dT}{dt}={\lambda }^{-1}\left(-e\sigma T^4+\frac{{\mu }_0{I}_0}{4}\left(1-A_0+m^{-1}\sum _{j=1}^m{c}_j{x}_j\right)\right).$$

$$A_0=A_{ice}{S}_{ice}{S}_{arc}^{-1},c_j\propto S_{free}=S_{arc}-S_{ice},$$

$$\frac{dT}{dt}={\lambda }^{-1}\left(-e\sigma T^4+\frac{{\mu }_0{I}_0}{4}\left(1-A_0+m^{-1}\sum _{j=1}^m{c}_j{X}_j(T)\right)\right).$$

$$\frac{d{y}_i}{dt}=\kappa g_i(y,x),$$

$$\frac{d{x}_j}{dt}=\sum _{l=1}^p{A}_{jl}{x}_l+{\kappa }_1{F}_j(y,x),$$

$$F_j(y,x)=\sum _{k=1}^m{b}_{jk}{f}_k(y,x).$$

$$Re\sigma (A)<-{\delta }_0<0.$$

$$x_{j}=\Phi_{j}(y,\kappa,\kappa_{1})=\kappa_{1}\Big{(}\sum_{k=1}^{m}c_{jk}f_{k}(y,0))+\tilde{X}_{j}(y,\kappa,\kappa_{1})\Big{)},$$

$$c_{ik}=-\sum _{j=1}^m(A^{-1}{)}_{ij}{b}_{jk}.$$

$$|{\tilde{X}}_j(\cdot ,\kappa ,{\kappa }_1){|}_{C^1(U_{\kappa ,{\kappa }_1})}\to 0(\kappa ,{\kappa }_1\to 0).$$

$$\frac{d{y}_i}{dt}=\kappa g_i(y,\Phi (y,\kappa ,{\kappa }_1)),$$

$$|\tilde{F}{|}_{C^1({\mathbb{D}}^n)}<ϵ.$$

$$\frac{d{y}_i}{dt}=g_i(y,\Phi (y)),$$

$$|g_i{|}_{C^2({\mathbb{B}}^n\times {\mathbb{R}}^p)}<C_g,$$

$${\Phi }_i(y)=m^{-1}\sum _{j=1}^m{c}_{ij}{f}_j(y),$$

$$|f_j{|}_{C^2({\mathbb{B}}^n)}<C_f,$$

$$y(0)=y^{(0)}.$$

$$c_{ij}\in (-R_0,R_0),$$

$$E(c_{ij}-\bar{c})(c_{kl}-\bar{c})=0ifi\ne korj\ne l.$$

$$\frac{d{\bar{y}}_i}{dt}={\bar{g}}_i(y),$$

Full text

††preprint: APS/PRE

I Introduction

Understanding of the mechanisms and scenarios of climate change as well its current and potential effects on ecosystems and biodiversity have been a focus of keen attention and intense research over the last few decades Chen et al. (2014); IPCC (2013, 2014). There is a general consensus that climate change will likely have an adverse impact on the ecological systems and population communities resulting in species extinction and a considerable biodiversity loss worldwide.

Whilst the top-down effect of climate on ecosystems is thus well established, relatively little attention has been paid to a possibility of an opposite, bottom-up effect that ecosystems may have on the climate. The mainstream of research often tends to consider the ecosystems and population communities as ‘biological actors on the physical stage’ Mortimer (1975) often disregarding possible feedback. Meanwhile, in planetary science, there is the concept of Gaia Lovelock (2000) that postulates the biosphere regulates its planetary climate to mitigate it for its own survival. While this hypothesis has been introduced quite long ago, current research in planetary and earth sciences inspires new applications of this hypothesis. The work Alcabes et al. (2020) is shown that even if a model exoplanet has significant climate perturbations then the Gaia concept is still acceptable (the original Gaia concept is based on a static planetary climate). Another work Lenton et al. (2018) supports the Gaia concept considering Earth’s biosphere stability over climate change through the existence of climate feedback loops and climate tipping points Lenton et al. (2008); Ashwin et al. (2012).

In this paper, we present a new mathematical realization of the Gaia hypothesis by a model of coupled climate-biosphere dynamics. We consider the effect that the biosphere of a planet may have on a planetary climate by changing the global energy balance through modifying the planetary albedo.

Modeling of physical processes in the climate system leads to difficult problems, involving complicated systems of partial differential equations for biological and chemical processes IPCC (2013). There exist climate models with different levels of realism; they can include thousands and even millions of equations, thousands of parameters to adjust. Usually, one investigates these models by computer simulations Hurrell et al. (2013). However, it is difficult to estimate the reliability of these computations, since it is connected with a difficult mathematical problem on the structural stability of attractors Katok and Hasselblatt (1995); Smale (1967). The theory of linear response of climate systems to perturbations Lucarini (2008) is based on the Ruelle theory of linear response for dynamical systems that holds on the formal hypothesis that the dynamical system is of the type axiom $A$ one. The last fact implies structural stability. However, S. Smale’s $A$-axiom systems Smale (1967) seldom appear in practical applications. The class of structurally stable systems is very narrow; this mainly includes systems with hyperbolic or almost hyperbolic behavior. One can expect, therefore, that the attractors of climate systems are not structurally stable: their topological structure can change under small perturbations. Therefore, one can expect that they can exhibit complicated bifurcations under small parameter perturbations. Possibly, an adequate approach is to take into account random fluctuations and study random dynamical systems. Indeed, the climate system as a complex system has a large number of interconnected and interacting subsystems, including the following: the atmosphere, the oceans, the biosphere, etc. Determination of how the dynamics of these subsystems change to reach equilibria of the entire system is the main problem of so-called conceptual climate models.

There are different types of conceptual climate models. Many of them are energy balance climate models and are defined by an ordinary differential equation describing energy conservation in the climate system. The most popular model is a zero-dimensional model North and Kim (2017) based on the theory of blackbody radiation determining global temperature changes due to the difference in incoming and outgoing solar radiation. This difference may be caused by changing of control parameters: surface albedo, greenhouse gas emission, and even the solar constant. The equilibria and the ideas how to find them by the bifurcation theory tools can be found here Fraedrich (1979).

In the context of Gaia hypothesis Ashwin et al. (2012); Lenton et al. (2008) here arises a key question: Why does climate stays stable over long time intervals (e.g. hundreds of thousands of years)? To answer this question, we consider conceptual climate models where the dynamical variables may be decomposed as slow and fast modes. Then for large times fast mode dynamics is captured by the slow dynamics on a stable slow manifold of a slow-fast system. The slow variables determine a long-term climate evolution under external factors whereas the fast modes may be associated with rapid factors.

The paper is organized as follows. In the next section, we introduce a planetary climate model with a biosphere component that arises from coupling between the conceptual zero-dimensional global energy balance model of climate dynamics and a generic ecosystem dynamics model (a multispecific population system living on multiple food sources). In Section III, we consider a general class of systems to which our model belongs and discuss the stability of those systems. We then show in Section IV that, in the case of our climate-biosphere model, a planetary climate remains stable with regard to a variation of the ecosystem model parameters as long as biodiversity is sufficiently large but it can lose stability (hence potentially resulting in regime shifts and a global climate change) if the number of species is small. A discussion and conclusions can be found in the last section.

II The model

The energy balance system is one of simplest climate models. It is defined by the following equation Fraedrich (1979) :

[TABLE]

where $\lambda$ is thermal inertia, $T$ is the averaged surface temperature, $t$ is time, and $A$ is the albedo of the surface. The left term characterizes the time-dependent behavior of the climate system. On the right hand side, the first term is the outgoing emission and the second term represents the incoming star’s radiation. Generally, incoming radiation to the planetary surface from a star is modified by a parameter, $\mu_{0}$, to allow for variations in the stellar irradiance per unit area, $I_{0}$ (the solar constant in case of the Earth), or for long-term variations of the planetary orbit van Woerkom (1953). On the other side, the outgoing emission depends on the fourth power of temperature, the effective emissivity $e$ and a Stefan-–Boltzmann constant $\sigma$.

This model can be coupled with the modeled biosphere’s dynamics as follows. The complete averaged albedo $A$ can depend on the biosphere state. For simplicity, we mostly focus our analysis on a single global ecosystem competing for several resources. We consider the following classical model:

[TABLE]

[TABLE]

cf. Huisman and Weissing (1999); Kozlov et al. (2017), where $x=(x_{1},x_{2},...,x_{n})$ are the species abundances, $m\gg 1$, and $v=(v_{1},...,v_{n})$ the resource concentrations. Here $\mu_{i}$ are the species mortalities, $D_{k}>0$ are resource turnover rates, and $S_{k}$ is the supply of the resource $v_{k}$, $\phi_{i}$ is the specific growth rate of species as a function of the availability of the resource (also known as Michaelis–Menten’s function). The coefficients $\gamma_{i}>0$ define self-limitation effects Roy and Chattopadhyay (2007). We assume that each of the resources $v_{k}$, $k=1,\ldots,n$, is consumed by all species so that the content of $k$-th resource in the $i$-th species is positive $b_{ik}>0$.

We consider general $\phi_{j}$ which are bounded, non-negative and Lipshitz continuous

[TABLE]

i.e., $\phi_{k}$ have a minimal smoothness, they are bounded and non-negative. The last restriction means that species consume resources.

Moreover, we suppose

[TABLE]

where $\partial{\bf R}^{m}_{+}$ denotes the boundary of the hyperoctant ${\bf R}_{+}^{m}=\{v:v_{j}\geq 0,\ \forall j\}$. Moreover, we suppose that

[TABLE]

This assumption means that as the amount of the $j$-th resource increases all the functions $\phi_{l}$ also increase.

Conditions (4) and (5) can be interpreted as a generalization of the well known von Liebig law, where

[TABLE]

(cf. Huisman and Weissing (1999)) where $r_{k}$ and $K_{kj}$ are positive coefficients, and $k=1,...,M$. The coefficient $r_{k}$ is the maximal level of the resource consumption rate by the $k$-th species and coefficients $K_{ki}$, $i=1,...,M$ define the sharpness of the consumption curve $\phi_{k}(v)$.

A simple way to couple climate subsystem (1) and the modeled biosphere defined by (2) and (3) is to suppose that the resource supply parameters $S_{k}$ depends on the surface temperature $T$. Moreover, we can suppose the albedo is a linear function of $x_{i}$:

[TABLE]

Finally, we obtain the following climate-biosphere system

[TABLE]

[TABLE]

[TABLE]

As an example, let us consider a model planet where the surface significantly covered by ice Shields et al. (2013) and the ice-albedo feedback is the main regulator of the planetary climate dynamics Cess (1978). Let the area of some region of the planet be $S_{arc}$, the area occupied by ice be $S_{ice}$ while the free ice area be $S_{free}$ Sudakov et al. (2015), where $S_{free}=S_{arc}-S_{ice}$. One can suppose that different species coexist in free ice domain and the averaged albedo of this domain is a linear combination of contributions of different species. Then we obtain

[TABLE]

where $A_{ice}$ is the albedo of the ice-covered area. This relation will be useful below.

Suppose that species populations $x_{i}$ and resources $v_{k}$ are fast variables, while the temperature $T$ evolves slowly. Such a situation arises if, for example, $\gamma_{i}>>1$ (see Kozlov, V. et al. (2019)). Then one can show that for large times $t$ $x_{i}(t)\approx X_{i}(T)$, where $X_{i}(T)$ are time averaged equilibrium species populations for fixed $T$ (see section IV). Then we obtain the following equation:

[TABLE]

Note that the equation (13) formally resembles the well-known ice-albedo feedback modification of the zero-dimensional energy balance model Fraedrich (1979).

When the system (9), (10) and (13) is regarded as a model of a particular biosphere, the choice of coefficients $c_{k}$ is determined by the environmental conditions at the given location and the corresponding species properties. Since we are aiming at building a global model, we want the eqs. (9), (10) and (13) to be applicable to any part of a modeled planet. Thus, we consider the coefficients unspecified. More precisely, we suppose that coefficients $c_{k}$ are random numbers described by certain probability distributions. We introduce these coefficients randomly assuming the randomness of the biological evolution.

In the coming section, we will consider a general class of slow-fast system with random coefficients, which includes the system (13) as a particular case.

III A general class of systems

III.1 A slow-fast system

In this section, we consider the following class of systems:

[TABLE]

[TABLE]

where $i=1,\ldots,n$, $j=1,\ldots,p$, and

[TABLE]

In these equations, the unknown vector-valued function $y(t)=(y_{1}(t),...,y_{n}(t))$ consists of slow components, the unknown function $x=(x_{1},...,x_{p})$ determines fast components, $\kappa,\kappa_{1}$ are small positive parameters, $g_{i},f_{k}$ are given smooth and uniformly bounded functions, $b_{jk}$ are bounded coefficients, and the square matrix $A_{jl}$ defines a linear operator $A$ with the spectrum $\sigma(A)$ such that

[TABLE]

Then for sufficiently small $\kappa,\kappa_{1}>0$ the system of equations (14) and (15) has a locally attracting smooth and locally invariant in an open neighborhood $U_{\kappa,\kappa_{1}}$ of $x=0$ manifold $\mathcal{M}$ defined by

[TABLE]

where

[TABLE]

Here $A^{-1}$ stands for a matrix inverse to $A$ and sufficiently smooth functions $\tilde{X}_{j}(y,\kappa,\kappa_{1})$ define small corrections such that

[TABLE]

Existence of $\mathcal{M}$ follows from the known results (for example, Henry (1981); Constantin et al. (1989a); Katok and Hasselblatt (1995)).

As a result, we obtain the following system for slow variables:

[TABLE]

where $\Phi(y,\kappa,\kappa_{1})=(\Phi_{1}(y,\kappa,\kappa_{1}),...,\Phi_{p}(y,\kappa,\kappa_{1}))$.

For consideration of the systems with random parameters we need to use arguments from dynamical system theory and the Hoeffding inequality, one of concentration inequalities.

Recall the basic concept of structural stability introduced by A. Andronov and S. Pontryagin in 1937 Smale (1980). Consider a smooth vector field $F$ on compact domain ${\mathbb{D}}^{n}$ of $\mathbb{R}^{n}$ with a smooth boundary (or on a compact smooth manifold $M$ of dimension $n$). Assume that $F\in C^{1}({\mathbb{D}}^{n})$ and consider all $\epsilon$-small perturbations $\tilde{F}$ such that

[TABLE]

Consider systems of differential equations $dx/dt=F(x)$ and $dx/dt=F(x)+\tilde{F}(x)$ and suppose that they define global semiflows $S_{F}^{t}$ and $S_{F+\tilde{F}}^{t}$ on ${\mathbb{D}}^{n}$. The system $dx/dt=F(x)$ is called structurally stable if there exists an $\epsilon_{0}$ such that for all positive $\epsilon<\epsilon_{0}$ trajectories of semiflows $S_{F}^{t}$ and $S_{F+\tilde{F}}^{t}$ are orbitally topologically conjugated (there exists a homeomorphism, which maps trajectories of the first system into trajectories of the second one). Roughly speaking, the original system is structurally stable if any sufficiently small $C^{1}$ perturbations of that system conserve the topological structure of its trajectories, for example, the equilibrium point stays an equilibrium (maybe, slightly shifted with respect to the equilibrium of non-perturbed system), or the perturbed cycle is again a cycle (maybe slightly deformed and shifted). We will refer the number $\epsilon_{0}(F)$ the structural stability constant of the system $dx/dt=F(x)$.

Note that structurally stable dynamics may be, in a sense, ”chaotic”. There is a rather wide variation in different definitions of ”chaos”. Chaotic (not periodic and no rest point) hyperbolic sets occur in some model systems Katok and Hasselblatt (1995); Smale (1980); Anosov (1995); Newhouse et al. (1978); Ruelle (1976, 1979, 2014).

III.2 Systems with random parameters

We consider systems (18), which arise, in a natural way, from systems decomposed in slow and fast variables. We will use the following notation. We denote by $EX$ the expectation of a random quantity $X$, and by $Var\ X$ its variance. Moreover, $\Pr[A]$ denotes the probability of a random event $A$. In this section, we formulate general principles on averaging with respect to the parameters that are applicable to fast-slow climate models.

Consider the following general system of differential equations:

[TABLE]

where $i=1,\ldots n$, $y(t)=(y_{1}(t),\ldots y_{n}(t))$ is a unknown vector-function, and $\Phi=(\Phi_{1},...,\Phi_{p})$, $\Phi_{l}(y)$ are functions, which will be defined below. Let ${\mathbb{B}}^{n}$ be a compact subdomain of $\mathbb{R}^{n}$ with a smooth boundary $\partial{\mathbb{B}}^{n}$. We suppose that $g_{i}(y,\Phi)$ are smooth functions uniformly bounded as are the first and second derivatives with respect to all variables $y,\Phi$:

[TABLE]

where $C_{g}$ is a positive constant.

We assume, moreover, that the functions $\Phi_{l}(y)$ are linear combinations of other functions $f_{j}(y)$ with random coefficients $c_{ij}$:

[TABLE]

We suppose that the $f_{j}$ are non-random, fixed functions and they have uniformly bounded derivatives

[TABLE]

where a positive constant $C_{f}$ is uniform in $i,j$.

For (20) we set the initial data

[TABLE]

Let the following assumptions hold:

Assumption 1

Let $c_{ij}$ be independent random quantities such that $Ec_{ij}=\bar{c}$, and, moreover, almost surely

[TABLE]

where $R_{0}>0$ is a constant uniform in $m$.**

As a consequence, one has that if

[TABLE]

Together with system (20) we consider the corresponding averaged system:

[TABLE]

where

[TABLE]

where $i=1,\ldots,n$, and $y(t)=(y_{1}(t),\ldots y_{n}(t))$ is a unknown vector-function, and $\bar{\Phi}_{i}(y)$ are averages of functions $\Phi_{i}(y)$ over the random parameters $c_{ij}$:

[TABLE]

We assume that there hold the following conditions:

[TABLE]

and

[TABLE]

where $e(y)$ is a normal vector to the boundary $\partial{\mathbb{B}}^{n}$ at the point $y$ directed inward on the domain ${\mathbb{B}}^{n}$. For the system (25) we set the same initial data (24). Condition (28) implies that the Cauchy problem (24) and (25) defines a global semiflow on the domain ${\mathbb{B}}^{n}$.

III.3 Main features of the systems with random parameters

For slow variables systems (18) we prove an averaging theorem assuming that $c_{ik}$ are random independent parameters (see the Appendix). This theorem asserts the attractor of the original system is close to the attractor of averaged one with a probability $Pr_{m}$, which is exponentially close to $1$ for large $m$. So, our main idea in the explanation the relative stability of climate is that a large number of independent factor can mutually cancel each other out. The probability $Pr_{m}$ satisfies an inequality that includes the number $\epsilon_{0}$ is a measure of stability under perturbations. If $\epsilon_{0}>0$ is small, i.e., the original system is weakly stable and conserves its dynamics only under very small perturbations, then estimate (1) makes a sense only for large $m>m_{0}(\epsilon_{0})$ (in fact, for bounded $m$ the right hand side of (1) is negative).

Moreover, the structurally stable system are seldom found in real applications (if we exclude the cases $n=1$ and $n=2$, where they are generic). According to basic result of S. Smale Katok and Hasselblatt (1995); Ruelle (2014), for dimensions $n>2$ structurally stable systems are not generic. To overcome this difficulty, we consider an approach, which allows us to show that solutions of the original system stay in a small neighborhood of a local attractor of the corresponding averaged system.

The stability of many dynamical regimes can be proved by using Lyapunov functions. Recall that $L(y)$ is a Lyapunov function of a system $dy/dt=g(y)$ in a domain ${\mathbb{V}}\subset{\mathbb{R}}^{n}$ if $L$ is at least $C^{1}$ smooth and $L(y(t))$ does not increases along trajectories $y(t)$ of the system:

[TABLE]

For example, if $y^{*}$ is a stable rest point of the system, then one can construct a $L(y)$ close to a quadratic form, which is Lyapunov function in a small neighborhood ${\mathbb{V}}$ of $y^{*}$ and

[TABLE]

for some $c>0$.

The next statement (see the Appendix) can be proven for the Lyapunov functions. If the averaged system defined by (25) has a Lyapunov function then the original system (20) has the same Lyapunov function a probability $Pr_{L}$, which is exponentially close to $1$ as $m$ large.

This theorem can be applied to the Budyko–Sellers energy balance system (13) as follows. Suppose that the averaged system is gradient-like (note that (13) enjoys this property). Let $\bar{A}$ be an attractor of the original system, which consists of stable equilibria. Suppose that all equilibria of the averaged system are hyperbolic. Then there exists a Lyapunov function $L_{(}y)$ such that

[TABLE]

for all $y\in{\mathbb{V}}(\bar{A})$ and some $\epsilon>0$, where ${\mathbb{V}}(\bar{A})$ is an open subset of the attraction basins of $\bar{A}$. This subset contains all points $y$ except for small $\delta$-neighborhoods of equilibria, where $\delta\to$ as $\epsilon\to 0$. Then with probability $Pr_{\delta,\epsilon,m}$ all original system also has the same Lyapunov function with analogous properties.

IV Stability of the coupled climate-biosphere system

Let us apply (1) and (2) to a system defined by (13). In the general case this system is complicated. To simplify the problem, we suppose that the $c_{i}$ are random independent quantities such that $Ec_{i}=\bar{c}$, and, moreover, we apply the approximation obtained in Kozlov et al. (2017); Sudakov et al. (2017); Vakulenko et al. (2018). We assume that the turnovers satisfy $D_{k}>>1$. Then

[TABLE]

Suppose that all species $X_{j}$ survive and have positive abundances. Then

[TABLE]

[TABLE]

where we use notation $f_{+}=max(f,0)$. Then eq. (13) take the form (we remove the terms the order $O(D^{-1})$)

[TABLE]

We apply (1) and (2), with $p=1$ and

[TABLE]

The averaged system takes the form

[TABLE]

where

[TABLE]

Let all $\phi_{i}(S)$ be uniformly bounded by a constant $a$, $\phi_{i}(S)<a$ for all $i=1,...,m$ and $S$. Then we find that, with a probability exponentially close to $1$, there exists a Lyapunov function defined by

[TABLE]

where

[TABLE]

Non-degenerate local minima of this function are steady states (local attractors) of the averaged system, and local extrema are saddle points or repellers of that system. If $\bar{c}$ is small enough, we have only a single local attractor $T=\bar{T}_{e}$. Our theorems (see the Appendix) imply that the original system then also has (with a probability close to $1$) a single local attractor $T=T_{e}(m)$ and $|T_{e}(m)-\bar{T}_{e}|\to 0$ as $m\to\infty$.

The situation dramatically changes if the condition $\phi_{i}<a$ is violated, say, one species dominates or if $m$ is small. Then it is impossible to guarantee that $|T_{e}(m)-\bar{T}_{e}|\to 0$. This means that biodiversity decreases can produce global climate changes.

To find possible bifurcations, we consider the simplest case when we are dealing with a single resource $v_{1}=v$ and the growth function are identical for all species, $\phi_{i}(v)=v(K_{i}+v)^{-1}$. We assume that $S(T)=S_{0}+S_{1}\Delta(T)$, where the coefficient $S_{1}$ defines an influence of temperature on the resource supply and

[TABLE]

This means that there exist an optimal temperature $T_{0}$ for species growth and a characteristic spread of this temperature $\sigma_{T}$. Then we obtain eq. (33) with

[TABLE]

and the equation for temperature steady state takes then the form

[TABLE]

where

[TABLE]

Note that for $C>0$ the species diminish the averaged planetary albedo. Moreover, simulations show that for small $m$ variations in the species parameters, for example, in $K_{i}$ can decrease albedo.

Depending on $C$ we have either a single root of (34), or three, when two roots give us local attractors and the third root is a saddle point, as it is shown on Fig. 1. So, we observe here a pitchfork bifurcation, which is essentially the same as for the ice-albedo feedback problem, see above.

It is interesting to understand as global warming affects the described bifurcation effect. Consider the cold planetary region and relation (12). We observe that a decrease of the area occupied by ice increase the coefficient $\bar{c}$ and decreases $A_{0}$ thus it reinforces the bifurcation effect and can lead to climate bifurcation.

V Discussion and conclusions

Understanding the planetary climate dynamics and identification of factors and processes that can affect its stability are problems of significant importance, in particular, because of their prominent effect on the biosphere functioning. There is growing evidence that the biosphere can have a variety of feedback loops to climate and a comprehensive understanding is only possible based on the analysis of the coupled climate-biosphere system. For example, in the Earth system, the perturbation of the carbon cycle Sudakov and Vakulenko (2013) or water-vapor Held and Soden (2000) cycle is one such feedback, but there are many more. Here we focus on the feedback induced by an interaction between the biosphere of an ice planet and climate. In this paper, we have endeavored to address this issue theoretically, in the framework of the Gaia hypothesis by considering a conceptual model of climate-biosphere dynamics arising from the coupling between a global energy balance model to a generic multi-specific model of population dynamics.

The climate-biosphere system is an extremely complex system and the corresponding mathematical model, even a relatively simple ’conceptual’ one, is usually too complicated for a comprehensive analytical study. A possibility of nontrivial model reduction lays in the observation that different processes often go with very different rates, i.e. take place on very different timescales. In particular, many complex systems, including climate models, have slow and fast components. According to classical results Constantin et al. (1989b), large time dynamics of fast modes are captured by a dynamics of slow modes on a slow invariant manifold. It is well known that even low dimensional systems exhibit complex bifurcations Engler et al. (2017); Kazi (2014); Sieber et al. (2018). Moreover, such models exhibit multistationarity, i.e., existence of many stationary states that, according to Emanuel (2002), provides the climate stability under variations of astronomical factors.

Referring back the Gaia concept, why, however, was the climate system stable over long periods of time in the past? To answer this question, we assumed that parameters of fast subsystems are random and mutually independent. Under such assumptions, we prove a general theorem on connection between attractors of averaged and original systems. If the attractor $\bar{A}$ of the averaged system has a low fractal dimension then, with a probability close to $1$, the attractor of the original system is close to $\bar{A}$. We think that this result may have applications for many different fields such as global network systems with unknown parameters, foodwebs, gene networks etc.

So, climate stability can be explained by the fact that many independent factors are canceled out. In our realization of the Gaia concept, the stability of the climate system is ensured by growing biodiversity. Our findings are summarized in Fig. 2. Interestingly, our analysis suggests a possibility of a positive feedback of the biosphere on the climate change. Consider a scenario of a slow change in the energy balance resulting, for instance, in a gradual increase of the mean temperature. Its is well known that such an increase will eventually results in species extinctions and biodiversity loss (see Arrow 1 in Fig. 2). Our results predict that, as long as the number of extinctions is not too large, the biodiversity loss will not have any notable feedback on the climate dynamics (Arrow 2 in Fig. 2). However, when the biodiversity loss becomes considerable, i.e. the number of surviving species becomes small, the failing biota will have a positive feedback on the climate resulting in its destabilization (Arrow 3 in Fig. 2). The global climate change resulting from this destabilization is likely to have a stronger negative effect on the biosphere, hence accelerating the extinction rate.

This model may be used to reconstruct and project climate change on the ice planets of the Solar System Schulze-Makuch et al. (2013); Russell et al. (2017) and some exoplanets Yang et al. (2017). Another possible application of our approach is paleoclimate modeling. For example, in the Cryogenian period, the planet has been transformed into so-called ‘snowball Earth’, where early life survived under the environmental stress William T. Hyde and Peltier (2000), became stable and even diverse van Maldegem et al. (2019). Our model may help to evaluate how biodiversity could contribute to global ice melting in this period.

Acknowledgments

This research is supported by the Swedish Foundation for International Cooperation in Research and Higher Education (STINT), Grant IB 2018-7517. This work was assisted by attendance as a Short-term Visitor at the National Institute for Mathematical and Biological Synthesis, an Institute supported by the National Science Foundation through NSF Award DBI-1300426, with additional support from The University of Tennessee, Knoxville. We also would like to thank the Mathematical Biosciences Institute (MBI) at Ohio State University, for helping initiate this research. MBI receives its funding through the NSF grant DMS-1440386. IS and SVP acknowledge the kind hospitality of Banff International Research Station (BIRS) for Mathematical Innovation and Discovery where they worked on this research. DL has been supported by the Government of the Russian Federation (the grant 08-08). IS also thanks University of Dayton Research Council Seed Grant, 2019.

Appendix

In this Appendix, constants $c$ and $C_{i}$ can depend on system parameters but are uniform in $m$ for large $m$. Note that we sometimes denote different constants by the same index if it does not lead to confusion. Our proving plan can be outlined as follows. To simplify our statement, we first prove three auxiliary lemmas, and then we state short demonstrations of theorems.

Probabilistic estimates. Let us fix some points $y^{(k)}\in{\mathbb{B}}^{n}$, where $k=1,2,...,M$ and $M$ is an positive integer, which will be adjusted later. Let us define the events ${\mathcal{A}}_{\epsilon,i}(k)$ by

[TABLE]

[TABLE]

where $Not\ B$ denotes the negation of $B$ and $\bar{g}_{i}(y)$ are defined by relation (26).

The next auxiliary lemma is elementary but useful.

Lemma 1

One has

[TABLE]

Proof. That lemma can be proved by de Morgan’s rule.

Furthermore, we use Chernoff bounds to estimate \Pr\big{[}{\mathcal{A}}_{out,\epsilon,i}(k))\big{]}. Let $C_{\bar{g},\Phi}$ be a Lipshitz constant of $\bar{g}$ with respect to the variables $\Phi_{1},...,\Phi_{p}$, i.e., for all $y\in{\mathbb{B}}^{n}$ and $i=1,...,n$

[TABLE]

where $|\Phi|=\max_{l}\ |\Phi_{l}|$. This constant $C_{\bar{g},\Phi}$ exists due to assumption (21) to $g$. Moreover, an analogous estimate holds for derivatives with respect to $y$: all $i,j$:

[TABLE]

Lemma 2

One has

[TABLE]

where

[TABLE]

Proof. Our the first step is to estimate differences $\Phi_{i}(y^{(k)})-E\Phi_{i}(y^{(k)})$. To this end, let us fix indices $i$ and $k$ and introduce $X_{j}$ by

[TABLE]

Then

[TABLE]

Our Assumption 1 on $c_{ij}$ implies that $X_{j}$ are independent random variables. Moreover, taking into account that $C^{1}$ - norms of $f_{j}$ are uniformly bounded we have

[TABLE]

Let us recall the Hoeffding inequality. Let $X_{j}$, $j=1,...,m$ be independent random variables strictly bounded in intervals $[a_{j},b_{j}]$, i.e. surely $X_{i}\in[a_{i},b_{i}]$. Let $\bar{X}=m^{-1}\sum_{j=1}^{m}X_{j}$ be the average of those quantities. Then (see Hoeffding (1963))

[TABLE]

Therefore, according to Hoeffding’s inequality for each $\epsilon>0$ we obtain

[TABLE]

where $l=1,...,p$.

The second step is as follows. Consider the events

[TABLE]

Let ${\mathcal{B}}_{\delta}=\prod_{l=1}^{p}{\mathcal{B}}_{\epsilon,l,k}$ Then, due to (1) and (42),

[TABLE]

Let us take $\delta=\epsilon/4C_{\bar{g},\Phi}^{2}$. Then, if the event ${\mathcal{B}}_{\delta}$ takes place, we have (because $g$ is a Lipshitz map with the Lipshitz constant $C_{\bar{g},\Phi})$ and by definition of $\bar{g}$) that

[TABLE]

i.e., the event opposite to $A_{out,\epsilon,i}(k)$ takes place. Now conditions (21) and estimate (43) lead to inequality (2) that completes the proof of the lemma.

Let us define now the events ${\mathcal{A}}_{out,\epsilon,i,j}(k)$ and ${\mathcal{A}}_{\epsilon,i,j}(k)$ by

[TABLE]

where

[TABLE]

and

[TABLE]

There holds the following Lemma:

Lemma 3

One has

[TABLE]

where $\tilde{C}$ is defined by (38).

The proof of (3) repeats the same arguments used in the proof of (2) so do not present it.

Demonstrations of Theorem 1 and Theorem 2

Theorem 1

Suppose condition (28) holds and that averaged system defined by (25) defines a global dissipative semiflow on the domain ${\mathbb{B}}^{n}$. Moreover, let us assume that averaged system (25) is structurally stable with a structural stability constant $\epsilon_{0}(\bar{g})$ and that system has an attractor $\bar{A}$. Then with probability $\Pr_{\bar{A}}$ the original system (20) also defines a global dissipative semiflow on ${\mathbb{B}}^{n}$, which has an attractor $A$ topologically equivalent to $\bar{A}$. The probability $\Pr_{\bar{A}}$ satisfies the inequality

[TABLE]

where $C_{1},C_{2}$ are positive constants uniform in $m$.

Proof. We use (1), (2) and (3) and the following auxiliary construction. The domain ${\mathbb{B}}^{n}$ has the dimension $n$ therefore we can cover it by $N(r\epsilon)\sim(r\epsilon)^{-n}$ balls $\Omega_{\epsilon,k}$ of the radius $\epsilon$ centered at some points $y^{(k)}\in{\mathbb{B}}^{n}$. Here $r$ is a positive constant uniform in $\epsilon$. We denote the union of all those balls by $U_{\epsilon}$, it is an open neighborhood of ${\mathbb{B}}^{n}$.

Let us consider the perturbation $\tilde{g}(y)=g(y,\Phi(y))-\bar{g}(y)$ and estimate the $C^{1}$ norm of $\tilde{g}$ on $U_{\epsilon}$. Suppose that all events ${\mathcal{A}}_{\epsilon,i}(k)$ and ${\mathcal{A}}_{\epsilon,i,j}(k)$ defined by (36) and (45), respectively, take place. Then

[TABLE]

Then, due to conditions (21) on $g$, and definition of $\bar{g}$ we have

[TABLE]

where a positive constant $C_{1}$ is independent of $m$. Therefore, for each $y\in{\mathbb{B}}^{n}$ one can find such point $y^{(k)}$ that there hold the estimates

[TABLE]

[TABLE]

Those last inequalities and (47) imply

[TABLE]

where $C_{2}$ is a positive constant. We set $r=1/2C_{2}$. Due to conditions (28) and (29) the vector fields $g$ and $\bar{g}$ are directed towards interior of ${\mathbb{B}}^{n}$ that allows us to apply now the definition of structural stability Ruelle (2014). Then for positive $\epsilon\leq\epsilon_{0}(\bar{g})$ the attractor of the original system is topologically equivalent to the attractor of the averaged system. Note that $\epsilon_{0}$ does not depend on $m$ and it is defined by the averaged system only.

Furthermore, we compute the probability that all the events defined by (47) take place by (1), (2) and (3). This finishes the proof.

Theorem 2

Suppose condition (28) holds and that the averaged system defined by (25) has a Lyapunov function such that

[TABLE]

where ${\mathbb{V}}$ is an open subdomain of ${\mathbb{R}}^{n}$ with a compact closure, and moreover,

[TABLE]

for a positive constant $C_{L}$. Then with the probability $\Pr_{L,\epsilon}$ the original system (20) has the same Lyapunov function such that

[TABLE]

The probability $\Pr_{L,\epsilon}$ satisfies the inequality

[TABLE]

where $\bar{C}_{1},\bar{C}_{2}$ are positive constants uniform in $m$.

Let us note that, similarly to the previous theorem, If $\epsilon>0$ is small, estimate (2) makes a sense only for large $m>m_{0}(\epsilon)$.

Proof. We apply the same idea used in the previous proof. The domain ${\mathbb{V}}$ can be covered by $N(r\epsilon)\sim(r\epsilon)^{-n}$ balls $\Omega_{\epsilon,k}$ of the radius $\epsilon$ centered at some points $y^{(k)}\in{\mathbb{B}}^{n}$. Here $r$ is a positive constant uniform in $\epsilon$. Let us introduce the functions

[TABLE]

Consider the events

[TABLE]

[TABLE]

Suppose that all events defined by (52) take place. Then

[TABLE]

Now we use the estimate

[TABLE]

where $Lip_{H}$ is a Lipshitz constant of $H$. Let us estimate that constant. By definition of $H$ one has

[TABLE]

Due to Assumption 1

[TABLE]

where

[TABLE]

The same estimate holds for the Lipshitz constant of $\bar{H}$. Therefore, (53) and (54) give

[TABLE]

where $C_{3}>0$ is a constant uniform in $m$. Let us set $r=1/4C_{3}$. Then condition (49) of (2) and (57) show that (50) is satisfied. Furthermore, to complete the proof, we compute the probability that all the events defined by (47) take place by estimates analogous to (1), (2) and (3).