## Abstract

We propose a generalization of the classical Goody model by taking into account greenhouse gas emission effects. We develop an asymptotic approach that allows us to obtain an expression for the greenhouse gas flux via the temperature and fluid fields. We show that there is a possible tipping point in atmospheric dynamics resulting from greenhouse gas emissions, where the climate system becomes bistable under sufficiently intensive greenhouse gas emissions.

## 1. Introduction

Over the past few decades, Earth’s warming climate has received great attention. In particular, it is interesting to estimate the effects connected with greenhouse gas emissions and propagation in the atmosphere. For example, methane is a dangerous greenhouse gas, and the problem of methane emissions from the Siberian tundra has been a subject of significant interest [1–5]. An interesting idea, proposed recently, is connected with the so-called methane hydrate gun in the Arctic shelf [6,7]. It was suggested that the famous Permian catastrophe was a result of fast methane emissions that perhaps occurred many years ago (‘methane hydrate gun’ hypothesis) [6–9].

To describe extreme events, the concept of tipping elements and tipping points has been proposed [10,11]. When Earth’s climate system crosses some threshold, it can trigger a transition to a new state at a rate determined by the climate system itself. This makes tipping points identical to the dangerous bifurcations in nonlinear dynamics. In addition to tipping points, the concept of tipping elements has been introduced, where well-defined subsystems of the climate can function fairly independently, and cause abrupt changes. The work of Lenton *et al.* [10] considers permafrost as a tipping element. Arctic permafrost is a large source of methane, and it connects to the atmospheric state through positive feedback. Positive feedback loops give mechanisms through which small perturbations can have large impacts. Thus, this feedback can lead the climate system to a tipping (bifurcation) point. The work in Shakhova *et al*. [6] provided experimental research into this problem and originally formulated the methane hydrate gun hypothesis. In its original form, the hypothesis proposed that the methane hydrate gun could cause abrupt runaway warming over a time scale less than a human lifetime [8], and may be responsible for warming events. However, this runaway methane hydrate breakdown may have caused a drastic alteration of the ocean and atmospheric environments in the past. Possibly, the Permian extinction event, when 96 per cent of all marine species vanished 250 Ma, is a result of such an extreme event. Discussions by Kerr [7] in *Science* can be briefly summarized as such ‘Arctic armageddon, needs more science, less hype’.

Here, we propose a new approach that makes this problem mathematically tractable and allows us to describe catastrophic bifurcations in the atmosphere induced by soil greenhouse gas sources. As a result of some asymptotic estimates, this approach allows us to obtain an explicit, simple relation for a critical emission level involving fundamental physical parameters. These parameters can be estimated theoretically and by experimental data.

To obtain these estimates, we use the classical Goody model [12–14] that serves as a paradigm of fluid mechanical stability and turbulence. Meteorologists, however, have not paid much attention to the Goody model because it appears to be too simple a model of atmospheric convection. Note that our approach is applicable to more complicated and realistic models. We consider the Goody model as a basis to simplify technical mathematical details in our statement.

We extend the Goody model, taking into account effects connected with greenhouse gas production, a chemical degradation of greenhouse gas molecules, and greenhouse gas diffusion and convection. This extended model also describes the influence of greenhouse gases on thermal radiative transfer. Although the model proposed is a mathematical idealization, such a model can be useful as a first approximation, and it is close to classical models for Rayleigh–Bénard convection, where many linear and nonlinear stability results have been established [15–18].

### (a) Main results and organization of the paper

In §2, we introduce the Goody model, which takes into account greenhouse gas production. In §3, we show that, under some conditions, the extended Goody model with greenhouse gas diffusion and convection is well posed, and solutions are bounded; in particular, the atmosphere kinetic energy is bounded. The model contains a phenomenological parameter and, in §4, we obtain some integral relations that allow us to estimate this parameter. In §§4 and 5, we discuss possible bifurcations. A part of bifurcations is well studied earlier for the Goody model, and the corresponding patterns are quite similar to usual Rayleigh–Bénard cells. In this case (§4), we estimate how the greenhouse gas emissions affect the dynamics close to the bifurcation point. We show that this influence can either accelerate or slow down the dynamics.

In §§5 and 6, we describe new bifurcations caused only by the greenhouse gas emissions. The bifurcation parameter is an intensity of emissions. In §5, we obtain a nonlinear equation that defines the critical emission level. Using this equation in §6, we find an asymptotic relation for the critical emission level that involves different physical parameters (diffusion and degradation coefficients for the greenhouse gas, the coefficient of heat transfer, radiative absorption coefficient and others).

In §7, we obtain a simple nonlinear equation for the surface temperature (here, we use a rough approximation that the planet is a plate, and the greenhouse gas emissions from the soil are homogeneous processes). This equation has a single solution for small emission levels and three solutions for a large greenhouse gas flux. Finally, we see that the climate system can become bistable under a great greenhouse gas outburst.

## 2. Systems under consideration

We start with the Goody model in a dimensionless form,
2.1
2.2
2.3
where **v**=(*u*(*x*,*y*,*z*,*t*),*v*(*x*,*y*,*z*,*t*),*w*(*x*,*y*,*z*,*t*)) is the fluid velocity with the vertical component *w*, *θ*(*x*,*y*,*z*,*t*) is the temperature field, *P* is the pressure, *Γ* is a parameter (dimensionless adiabatic lapse rate), *σ* is the Prandtl number, *σ*_{1} is the buoyancy parameter, *F* is a radiative flux and *Θ*_{0} is the reference temperature. We consider, for simplicity, the motion in the flat layer *Ω* defined by 0<*z*<*h*, (*x*,*y*)∈*Π*, where *Π* is a rectangle (many relations hold in a general case where (*x*,*y*) lie on a platform *Π* of arbitrary form).

Actually, it is the standard Oberbeck–Boussinesq approximation taking into account the radiative flux. An equation governing the radiative flux *F* can be derived by making the Eddington (two-stream) approximation to Schwartzschild’s equation of radiative transfer, as in Goody & Young [19] or Goody [12]. The result is
2.4
where we assume that *α* is a constant proportional to the coefficient of absorption of radiation per unit volume. We suppose that *α* is small enough, whereas *θ*=*O*(1). Then, equation (2.4) can be simplified
2.5
and (2.2) takes the form
2.6
Bifurcations in this model have been studied extensively by Larson [20]. Equation (2.6) does not contain *F*.

Now, we extend this model by an equation for greenhouse gas convection, degradation and diffusion,
2.7
where *C*(*x*,*y*,*z*,*t*) is the greenhouse gas concentration, *d* is the diffusion coefficient and is a positive constant (the term describes a natural degradation of greenhouse gas molecules in the atmosphere owing to chemical reactions).

Note that one can use a more realistic equation for *θ* and *C*. The operator *Δ* can be replaced by a more general nonlinear diffusion operator,
2.8
(see Olbers [21]). Here, coefficients *K*_{h} and *K*_{v} of the horizontal and vertical diffusion may depend on *θ*. In the simplest form of (eddy) diffusion, these coefficients are taken to be constant [21], and all our mathematical analysis is valid, but formulae become more complicated.

Equation (2.7) can be resolved if we know the fluid flux *v* and formally, equations for *θ*,**v** and for *C* are independent. To describe the greenhouse gas influence on fluid circulation, we assume that the coefficient *α* depends on *C*: *α*=*α*(*C*). Because the greenhouse gas concentration *C* is small, it is natural to assume that this dependence is linear,
2.9
where *α*_{0} and *α*_{1} are constants. Here, *α*_{1} is a phenomenological coefficient that can be defined by an analysis of experimental data on the greenhouse gas mass in the atmosphere by some relations that will be obtained below.

Then, equations (2.1), (2.3), (2.6) and (2.7) become interconnected, and it is possible to compute feedback between the fluid and temperature fields as well, as the greenhouse gas concentration. These equations form our main system that will be under consideration.

Let us consider the boundary conditions. For **v,** we set the standard no-slip conditions
2.10
The boundary conditions for the temperature *θ* at the upper boundary *z*=*h* have the form
2.11
The boundary conditions for the temperature *θ* at the surface *z*=0 have the form
2.12
*κ*_{0} is a positive parameter defined by [21]
2.13
where *ρ* is the density, *K*_{v} is the thermal conductivity along the vertical *z*-direction, *c*_{p} is the specific heat capacity and
2.14
Here, *Q*_{SW} is the incident energy flux of shortwave radiation (which can be computed from a radiation model that is an essential part of a full climate model, see [21], §4.2), *α*_{s} is the surface albedo, the term follows from the Stefan–Boltzmann law, where *ϵ** is the emissivity of the surface, *σ*_{SB} is the Stefan–Boltzmann constant and *T*_{s} is the surface temperature. We can set *T*_{s}=*θ*(*x*,*y*,0,*t*).

The boundary condition for *C* describes the absence of the greenhouse gas flux at *z*=*h* and the greenhouse gas production at the bottom,
2.15
and
2.16
where the function *μ* describes the greenhouse gas flux intensity at the bottom. We assume that this flux is inhomogeneous in space and depends on the temperature. To simplify our mathematical analysis, we set periodic boundary conditions with respect to *x*,*y*. If *Π*={(*x*,*y*):0<*x*<*l*_{1}, 0<*y*<*l*_{2}}, then these conditions read
2.17
2.18
2.19

Depending on the function *μ* form, there are different possible situations. We can assume, for example, following Komarov [3], that there holds a relation of Arrhenius type
2.20
where *V* _{0} is a potential barrier and *k*_{B} is the Boltzmann constant. For example, if we consider lakes as a main greenhouse gas (basically, methane) source, we can assume that
2.21
where *c*_{g}>0 is a constant, *Ω*_{i} is a two-dimensional domain occupied by the *i*th lake and *χ*_{V} denotes the characteristic function of the set *V* : *χ*_{V}=1 if (*x*,*y*)∈*V* and *χ*_{V}=0 otherwise.

Note that *Ω*_{i} can extend; therefore, we should consider the case when *g* depends on time. However, this time evolution is slow, *Ω*_{i}=*Ω*_{i}(*τ*), where *τ* is a slow time *τ*=*ϵ*_{1}*t*, *ϵ*_{1}≪1, and we can assume below that *Ω*_{i} are fixed domains.

Equation (2.20) is consistent with some phenomenological models. We mention here the following model. The first is a result of observations under a large peatland ‘Bachkarevskoe’ within 1994–2002 [4]. The following empirical formula was proposed (for a similar relation, see [22]):
2.22
where *F* is the flux measured in Mg m^{−2} per hour, *T*_{12} the temperature at 12 cm depth, *W* is the water level. This formula can be obtained by a Taylor series from the Arrhenius relation (2.20) that allows us to obtain an approximative value of *V* _{0}/*k*_{B}.

## 3. Problem is well posed: solutions are bounded

Let us describe briefly how to prove global existence and uniqueness of solutions for the system of equations (2.1), (2.3), (2.6), (2.7) under the boundary conditions formulated earlier. Different approaches to this initial boundary problem are well known [23–25]. The local (in time) existence and uniqueness follow from these standard results. Global existence (for all times *t*) can be obtained by local existence and an estimate of a weak norm (associated with an energy, for example). From the physical point of view, it is important to show that the energy of the fluid flow is bounded for all times as well as the amplitudes (the norm ) of *θ*. We can obtain *a priori* estimates for the energy and the amplitudes only for the case *Γ*=0.

### (a) *A priori* estimates

To show that the initial-valued problem is well posed, and solutions exist for all *t*, we have to obtain *a priori* estimates of the solutions. We suppose the following.

### Assumption 3.1

*The function* *is positive for sufficiently large θ,
*
3.1
*and
*
3.2
*The function μ satisfies
*
3.3

Note that conditions (3.3) and (3.1) admit a transparent physical interpretation. The first one means that the greenhouse flux is directed from the soil (the surface: *z*=0) into the atmosphere. The second one holds in the case (2.14) because the Stefan–Boltzmann term is much greater than the albedo contribution for large *θ*. This means that the surface absorbs the heat for high temperatures. Conditions (3.1) and (3.2) are fulfilled for *q* defined by (2.14) if the albedo *α*_{s}∈(0,1).

Under this assumption, the estimates proceed as follows. First, we show that *C*(*x*,*y*,*z*,*t*)≥0, if initial concentration *C*(*x*,*y*,*z*,0)>0. Second, we prove that *θ*(*x*,*y*,*z*,*t*)≥0, if *θ*(*x*,*y*,*z*,0)>0 and is bounded. Then, it follows that ∥**v**(⋅,⋅,*t*)∥ is bounded uniformly in *t*. Here,
denote standard *L*_{2} and supremum norms.

The assertions about *C*,*θ* can be proved in an elementary way by the maximum principle (see Friedman [26]). Let us consider a time interval *t*∈[0,*T*]. Assume *C*(*x*,*y*,*z*,*t*)<0 holds at a point (*x*,*y*,*z*,*t*), 0<*t*≤*T*. Let (*x*_{0},*y*_{0},*z*_{0},*t*_{0}) be a point, where 0<*t*_{0}≤*T* and where *C* has a minimum. Denote **x**=(*x*,*y*,*z*) and **x**_{0}=(*x*_{0},*y*_{0},*z*_{0}). Then, *C*_{0}=*C*(**x**_{0},*t*_{0})<0. Because *μ*<0, *C*_{z}(0)<0, therefore *z*_{0}>0. Then, we have
and (2.7) leads to a contradiction, therefore, *C*(*x*,*y*,*z*,*t*)≥0 for all (*x*,*y*,*z*,*t*). This implies that *α*_{0}+*α*_{1}*C*>0. Then, in a similar way, we can show that *θ*≥0 (here, we use *q*(0)<0).

Let us show that *θ*(*x*,*y*,*z*,*t*)≤*M* for all (*x*,*y*,*z*)∈*Ω* and *t*>0 for some *M*>0. Let us take a number *M*>0 such that and *M*>*θ*_{0} (where *θ*_{0} is defined by (3.1)). Assume *θ*(*x*,*y*,*z*,*t*)>*M* is fulfilled at a point (*x*,*y*,*z*,*t*), 0<*t*≤*T*. Let (*x*_{0},*y*_{0},*z*_{0},*t*_{0}) be a point, such that 0<*t*_{0}≤*T* and *θ* has the maximum at **x**_{0},*t*_{0}, where **x**_{0}=(*x*_{0},*y*_{0},*z*_{0}). Then, *θ*(**x**_{0},*t*_{0})>*M*. Because *q*(*M*)>0 and *θ*_{z}(0)>0, one has *z*_{0}>0. One obtains, as above,
and heat equation (2.6) leads to a contradiction (under our assumption *Γ*=0); therefore, *θ*(*x*,*y*,*z*,*t*)≤*M* for all (*x*,*y*,*z*,*t*).

Note that *C* can increase as a linear function of time (thus, blow-up effects are impossible). In fact, let us denote by 〈〈*f*〉〉 the value of function *f* averaged over all the atmosphere: and |*Ω*|=*hL*_{1}*L*_{2} is the volume of the elementary periodical box. Then, integrating (2.7) over all *Ω*, one has
which, according to assumption (3.1), gives 〈〈*C*〉〉≤*C*_{0}+*M*_{1}*t*, where *C*_{0} is the averaged concentration at the initial moment.

Because |*θ*| is bounded, the fluid equation entails the standard energetic estimate [23,24]
where *M*_{3}>0 is a constant. Let us use now the Poincaré inequality ∥**v**∥≤*c*∥∇**v**∥ that gives
with *σ*_{1}>0. Therefore, by the Gronwall lemma,
This proves that ∥**v(t)**∥ is bounded in *t*.

Therefore, we have shown that our initial-value problem is well posed. To conclude this section, let us note that the greenhouse concentration can be increased as a linear function of time whereas the fluid energy and the temperature stay bounded. The physical reason of this effect is the presence of the Stefan–Boltzmann term in the boundary conditions at the soil.

### (b) Relations for mean temperature time evolution

Let us obtain an integral relation for the mean temperature 〈〈*θ*〉〉. Let us integrate the heat equation (2.6) over *Ω*. Integrating by parts and taking into account boundary conditions (2.7) and (2.11), one obtains
3.4
where *Q*_{s} is the heat flux from the surface, , *S* is the surface area and *q* is defined by (2.14).

The emission greenhouse gas contribution in warming is then given by 3.5

Using this relation and having experimental data on the total atmospherical greenhouse gas mass, one can find a rough estimate of the coefficient *α*_{1}. We can simplify relation (3.5). Owing to the atmosphere turbulent mixing, one can assume that
Then, the value 〈〈*θ*〉〉 can be obtained from (3.4): 3*α*_{0}〈〈*θ*〉〉=*Q*_{s}/*S*=*q*(*T*_{s}). The value *W*_{met} can be found from experimental data; therefore, relation (3.5) allows us to find *α*_{1}.

## 4. Greenhouse gas emission influence on bifurcations

The Goody model describes bifurcations, leading to patterns similar to Rayleigh–Bénard cells [27,16,17,20,28]. In this section, using the extended Goody model, we investigate the influence of greenhouse gas emissions on these bifurcations.

Fundamental variational relations, first obtained by Sorokin [15] and others (see also [17,29]), can be applied to this problem in order to estimate how the greenhouse gas can affect the fluid dynamics. Unfortunately, the possibilities of the analytical technique are restricted to a small neighbourhood of the bifurcation point; nonetheless, some results can be obtained and, possibly, extrapolated to general fluid flows.

We follow the standard scheme [17]. To start, we remove nonlinear terms from the Navier–Stokes equations (2.1)–(2.3).

We postulate a basic state in which **v**_{0}=**0** and the temperature . In this section, we denote by *θ*′ the temperature deviation with respect to this basic state *θ*_{0}, and *P*′ denotes the pressure perturbation. We assume here that, when we analyse the atmosphere dynamics, this basic state can be considered as time independent. We then obtain the following linearized equations for small corrections **v**,*θ*′,*P*′ to the basic state
4.1
4.2
4.3
where we assume that is an almost stationary greenhouse gas concentration slowly evolving in time *t*. We assume that all unknown functions are *l*_{1}-periodic in *x* and *l*_{2}-periodic in *y*. From boundary conditions (2.10)–(2.12), we obtain the following boundary conditions for perturbations:
4.4
4.5
4.6
Here, the boundary value *θ*_{0}(0)=*a* satisfies *q*(*a*)=0.

For the case *α*_{1}=0 (thus, *α*(*C*)=*α*_{0} is a constant), the fluid velocity and thermal fields can be defined, up to small corrections, by the normal mode solutions
4.7
where **V**=(*V*,*U*,*W*),*T* are eigenfunctions of some linear operator acting on corrections to **v**,*θ*. Here, *X*(*t*) is an unknown time function determining the time evolution of the correction fields. If we consider the linearized equations (removing the influence of a small external force and small effects with *α*_{1}*C*), then , where *λ* is, in general, a complex eigenvalue *λ*=*β*+i*ω*, where . The functions *W*,*T* can be decomposed as follows:
4.8
For *f*, we have
4.9
where is the horizontal Laplacian.

Furthermore, one can obtain ordinary differential equations for *W*(*z*) and *T*(*z*) (see [30,17]). The bifurcation point is defined by the condition *β*=0.

In the general case, when *α*_{1}≠0, we can obtain some relations for *λ* following Gershuni & Zhukhovickij [17]. For **V** and *T*, we have
4.10
4.11
4.12
We multiply (4.10) by the complex conjugated **V*** and (4.11) by *T**. Then, we multiply equations conjugated to (4.10) and (4.11) by **V** and *T* and obtain the following relations:
4.13
and
4.14
where d*v*=d*x* d*y* d*z* and *γ*′=−*γ*−*Γ*, and a boundary contribution *I*_{S} is defined by
4.15

If *γ*′>0, then one can show again, following Gershuni & Guchovizkii [31], that *Im* (*λ*)=*ω*=0. Then, using (4.14) and (4.15), one obtains
4.16
where
4.17
This relation for the term *H*_{0} is the key result of this section. For small *α*_{1}, this term describes the shift of the eigenvalue *λ* owing to the greenhouse gas emissions.

The quantity *H*_{0} can be computed in some cases when we can obtain some expressions for the time-averaged greenhouse gas concentration *C*. Averaging the equation *C* over time and assuming (i.e. correct at the bifurcation point) that the convective term is small, one has
4.18
under the boundary conditions
4.19
and
4.20
where the function *μ* was defined earlier.

If *b*_{0} is small with respect to , or convection induces a complicated dynamics, then one can assume that the greenhouse gas concentration is almost uniform in space, and we obtain
4.21
where is the averaged greenhouse gas concentration in the atmosphere. For perturbations *δ*_{met}*λ* of the parameter *λ*, we obtain
4.22
where .

Note that the integral in (4.21) always positive. We observe, therefore that if the atmosphere state is close to a bifurcation point, the greenhouse gas emission influence on the atmosphere dynamics depends on the sign of the term *K*. Note that Re(*λ*) determines damping of the main mode. This greenhouse presence can accelerate the atmosphere dynamics, or, inversely, it can slow down this dynamic. Indeed, if *K*>0, then the perturbation *δ*_{met}*λ* of *λ* is negative, and we have a slow down. Otherwise, *δ*_{met}*λ*>0, and we obtain the acceleration of the main mode that evolves in time, roughly, as .

These effects depend on the sign *γ*′=−*γ*−*Γ* and the atmosphere state. Remember that *γ* defines a linear non-perturbed temperature profile *θ*_{0}(*z*). The temperature profile has been studied experimentally [13], and as it turned out, it is a more complicated function of *z*. It is natural to assume that *γ*<0 because we consider the lowest atmosphere layer where the temperature decreases in *z*. Thus, it is possible that *γ*′>0 and then . In this case, *K*>0 if, for example, an integral contribution ∥**V**∥ of the speed velocity perturbations is much more than the integral contribution ∥*T*∥ of the temperature perturbations. If, otherwise, ∥*T*∥≫∥**V**∥, we have *K*<0.

## 5. New bifurcations in extended Goody models

We describe here new bifurcations that are possible in our extended Goody models, and impossible in the standard Goody model. These bifurcations appear as a result of atmosphere warming by the greenhouse gas. A straightforward analytical procedure allows us to compute the critical level of the greenhouse gas emissions.

We linearize the main equations at zero state **v**=**0**,*θ*=*T*_{0}(*z*),*C*=*C*_{0} assuming, for the sake of simplicity, that *μ* in boundary condition (4.20) is a function of *θ* and independent of *x*,*y*. Below, we denote by *θ* the temperature deviation with respect to the base state *T*_{0}. Then, the linearized equations (2.6) and (2.7) for small perturbations *θ*,*C* take the form
5.1
and
5.2
where *α*_{0},*b*_{0}>0. We consider these equations in the layer *Ω*={0<*z*<*h*, *x*∈(−*L*_{1},*L*_{1}), *y*∈(−*L*_{2},*L*_{2})} and assume that *h*≪*L*_{1},*L*_{2}. Let us set the simplest boundary conditions describing a uniform greenhouse gas emission,
5.3
5.4
5.5
5.6

The last condition means that we consider a simple linear approximation: emission is proportional to the temperature deviation. Here, *β*=d*μ*(*θ*)/d*θ*|_{θ=T0(0)}. We have the greenhouse gas flux inside the atmosphere if *β*<0. The coefficient *r*_{0} is a result of linearization (2.12): *r*_{0}=(d*q*/d*θ*)(*T*_{0}(0)).

We seek solutions where, in general, *λ* is a complex spectral parameter. Then,
5.7
and
5.8
Because *T*_{0}=*T*_{0}(*z*), we can use the Fourier decomposition that gives
5.9
where *k*=(*k*_{1},*k*_{2}) is a wavenumber.

Let us set *X*=*T*_{k}(0). Given *X*≠0, let us find equations and boundary conditions for vertical distributions *T*_{k},*ϕ*_{k}. An easy computation then gives
5.10
and
5.11
where , . The boundary conditions take the form
5.12
5.13
5.14
5.15
The solution of boundary problems (5.11), (5.14) and (5.15) is given by
5.16
To find *T*_{k}, we have to resolve the boundary problems (5.10), (5.12) and (5.13) under the additional restriction *T*_{k}(0)=*X*. Owing to this restriction, this problem has no solutions for general *λ*. However, we can find such a non-trivial solution for some *λ*. These special *λ* satisfy a nonlinear equation that can be derived using a special trick.

Note that relation (5.16) for *Φ*_{k} gives
where
Multiplying equation (5.10) by , integrating by parts and taking into account boundary conditions (5.12) and (5.13) for *T*_{k}, one finally has
5.17
This nonlinear equation should be resolved for each *k* and have roots depending on *λ*=*λ*(*k*). Both the left- and right-hand sides of (5.17) depend on *k* only through the wavenumber module |*k*|, therefore, *λ*=*λ*(|*k*|). The case *k*=0 is admissible.

We consider *β*<0 as a bifurcation parameter. For each , we can find a *β*=*β*_{c} such that Re *λ* (*k*,*β*)=0 for *β*_{c} and Re *λ* (*k*,*β*)>0 for *β*<*β*_{c}<0.

Equation (5.17) is complicated and, to simplify computation of *β*_{c}, we apply an asymptotical approach that we will state in the following sections.

## 6. Dynamics and critical emission levels

### (a) Asymptotic calculation of *T*,*ϕ*

In the general case, the critical value *β*_{c} can be computed only numerically because the analytical expression for *T* has a complicated form. However, it is important to obtain an analytical expression for *β*_{c} in order to see how this value depends on the main physical parameters *d*,*β*,*α*_{0},*K*,*b*_{0} and *α*_{1}.

We thus consider a particular case when one can compute asymptotics for *β*_{c}. Let us set, to simplify calculations, *k*=0 (numerical simulations show that the bifurcation with a minimal critical level *β*_{c} occurs for the minimal |*k*|=0). Our main assumption is that . From the physical point of view, this means that the greenhouse degradation is a much weaker effect than the eddy diffusion (induced by convection and turbulence). Then, because *k*_{c}≈0 for *k*=0, and small , the greenhouse gas distribution of *Φ*=*Φ*_{k}(*z*) along *z*, given by (5.16), is almost uniform. We replace the basic temperature profile *T*_{0}(*z*) by a constant , the mean value of *T*_{0} along *z*∈(0,*h*). Then, the solution *T*_{k}(*z*) of (5.10) has the form *T*_{k}(*z*)≈1. In fact, without loss of generality, we can assume that *X*=*T*_{k}(0)=1. Function *f*_{k}(*z*) in (5.10) is close to a constant; thus, the solution *T*_{k}(*z*) is also a constant (such a solution satisfies boundary conditions (5.10), (5.12) and (5.13)). Consequently, we obtain *T*_{k}(*z*)≈1.

### (b) Critical value of *β*

Now we can find the critical value of *β*. We compute the integral on the right-hand side of (5.17) under the assumptions of §6*a* (i.e. *k*_{c}(*λ*)≪1). An easy asymptotic estimate of the integral gives
Remember that we suppose *β*<0. This equality allows us to find a critical *β*_{c}. Let us denote *k*_{T}=*k*_{T}(0) and *k*_{c}=*k*_{c}(0). We finally have
6.1

### (c) Taking into account convection effects

In the analysis carried out earlier, we removed the fluid velocity **v**, setting **v**=0. This assumption of course is not too realistic.

Let us consider equations
6.2
and
6.3
where we take into account convective contributions presented by terms with **v**^{0}, where **v**^{0} is a given velocity field. This means that we consider an influence on bifurcations of some time-stationary fluid flow. There are three key points.

— The flow

**v**^{0}does not influence energetic relations for real*λ*. In fact, if the eigenvalue*λ*is a real number, then (6.2) and (6.3) imply (we multiply these equations by*θ*and*ϕ*, respectively) 6.4 and 6.5 and we see that the flow is not involved in this relation.— Earlier, we found that, in the absence of the flow, the most dangerous bifurcations occur for

*λ*=0 and the wavenumber**k**=**0**. This means that the main eigenfunctions*θ*=*T*(*z*),*ϕ*=*Φ*(*z*) depend only on the vertical coordinate*z*. This yields that the flow**v**^{0}does not shift the bifurcation parameter*β*_{c}if this flow is purely horizontal and we neglect the vertical component:*W*=0. In fact, then**v**^{0}⋅∇*Φ*=0 and**v**^{0}⋅∇*T*=0, i.e. the corresponding contributions vanish.— One can eliminate the fluid velocity

**v**using the idea of eddy (turbulent) diffusion. This means that we replace the terms*dΔC*+(**v**⋅∇)*C*by an effective diffusion operator, for example, of the form (2.8). In this case, the bifurcation analysis can be carried out in a similar way.

### (d) Parameter estimates in an explicit relation for a critical emission level

We shall obtain here an explicit relation for the emission critical value *β*_{c}. In (5.17), we set *k*_{T}*h*,*k*_{c}*h*≫1 and *r*_{0}=0. We assume, to simplify calculations, that *T*_{0}(*z*)=*T*_{0}=const. Removing the small terms, we obtain
6.6
The parameters involved in this relation can be estimated as follows.

— The parameter

*K*is thermal diffusivity, we suppose that*K*∈0.15×10^{−8}m s^{−2}.— defines a rate of the natural greenhouse gas concentration

*C*decay. It is known that and*C*falls twice within 10–12 years, therefore –0.5×10^{−9}s^{−1}.—

*d*is the greenhouse gas diffusion coefficient that we take from 10^{−6}to 10^{−7}m^{2}s^{−1}.— A first approximation for averaged over space temperature is . We assume that the averaged heat flux in the atmosphere is given by the solar constant

*Q*_{sun}. The solar constant includes all types of solar radiation, not just the visible light. It is measured by a satellite to be roughly 1.366 kW m^{−2}. The value of*q*in our balance equation is*q*=*Q*_{sun}/*ρc*_{p}, where*ρc*_{p}is the volumetric heat capacity. We then obtain a rough estimate of*α*_{0}.

To obtain a rough estimate of sensitivity *α*_{1}, we integrate the heat transfer (2.6) over time *t* within *t*,*t*+*δt* and over all space, assuming, to simplify, that the concentration *C*=const. This gives that the corresponding decrement , induced by the greenhouse gas emission term, is given by *DT*=*α*_{1}*C*, and knowing *C*, we can evaluate *α*_{1}.

## 7. A nonlinear equation for surface temperature

In this section, we derive a nonlinear equation for *T*_{s} and show that, for a sufficiently large greenhouse gas emission from the soil, the climate system becomes bistable.

We consider a very simple, plate planet where the surface is homogeneous, and thus the greenhouse gas emissions are also homogeneous in *x*,*y*. We seek solutions to (2.1), (2.3), (2.6) and (2.7) in the following form.

—

*C*and*θ*depend only on*z*:*C*=*C*(*z*),*θ*=*θ*(*z*) and therefore*T*_{s}=*θ*(0) is a constant.—

*W*=0, and we neglect convective terms. Moreover,*μ*depends only on*θ*.— We replace the term

*α*_{1}*CT*by , where defines a vertical temperature profile, and therefore, .

Then, our initial boundary problem reduces to a nonlinear boundary problem, 7.1 7.2 7.3 7.4 7.5 7.6

Here, *μ*(*T*_{s}) and *q* are defined in §2 (see relation (2.14) and below (2.20)). The solution of equations (7.1)–(7.3) has the form
7.7

The solution *θ*(*z*) of boundary problems (7.4)–(7.6) has to satisfy the natural condition *θ*(0)=*T*_{s}. This allows us to find an equation for *T*_{s} without finding *θ* in an explicit form. Multiplying equation (7.4) by , where *k*_{0}=(3*α*_{0}/*K*)^{1/2}, integrating by parts and taking into account boundary conditions (7.5) and (7.6), one finally has the following equation:
7.8
where *S*=−*T*_{s}*μ*(*T*_{s})>0 and

Via *q* and *μ*, equation (7.8) involves the important physical parameters, namely the solar constant, the albedo (*ϵ**) and also the Stefan–Boltzmann constant. We consider here *ξ* as a parameter that determines the emission level. Note that function *S* has a sigmoidal form, typical for the neural network theory, and, thus, can be used to describe an excitable system.

A simple analysis shows that, for small emission levels, *ξ*<*ξ*_{0}, equation (7.8) has a unique solution, whereas for large *ξ*, (7.8) can have three solutions. In fact, for small *ξ*, the function *R*(*T*_{s})=*T*_{s}*k*_{0}−*ξS*(*T*_{s})+*q*(*T*_{s}) is increasing in *T*_{s}. This can be illustrated by figure 1, where for *S*(*T*_{s}), a sigmoidal approximation is taken, , that describes a trigger regime of excitation at *T*_{s}=300 K. Therefore, we conclude that there is possibly a bistable regime for our toy model of climate.

## 8. Concluding remarks

In this paper, we have considered a simplified approach that allows us to take into account the influence of greenhouse gases on atmospheric dynamics. This approach uses the Goody model coupled with an equation for greenhouse gas diffusion and convection. The coupling is achieved by introduction of a small nonlinear term proportional to the greenhouse gas concentration. The main goal of this approach is to verify the hypothesis that an outburst of a greenhouse gas can lead to a climate bifurcation. Indeed, we show that such a model has a large spectrum of bifurcations, including new ones generated by radiation effects and greenhouse gas emissions. A mathematical analysis leads to some nonlinear equations involving many physical parameters. However, they are not differential equations, and are quite tractable mathematically. Under some natural assumptions, it is possible to obtain a relation for the critical emission level leading to a bifurcation. Interestingly, the corresponding pattern is always the same, with our ‘toy planet’ sinking slowly into a homogeneous greenhouse gas fog.

## Acknowledgements

This research was presented at the IMA Conference on Mathematics of the Climate System (University of Reading, UK, September 2011). The authors are grateful to the Institute of Mathematics and its Applications for financial support for I. S.’s participation in this conference (IMA Small grant no. GS30/11). S.V. is grateful to the Russian Foundation for Basic Research for financial support. We thank Prof. Kenneth M. Golden (University of Utah, USA) for great help. We also are grateful to the referees for interesting and useful remarks.

## Footnotes

One contribution of 13 to a Theme Issue ‘Mathematics applied to the climate system’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.