the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Variational techniques for a onedimensional energy balance model
Gianmarco Del Sarto
Jochen Bröcker
Franco Flandoli
Tobias Kuna
A onedimensional climate energy balance model (1D EBM) is a simplified climate model for the zonally averaged global temperature profile, based on the Earth's energy budget. We examine a class of 1D EBMs which emerges as the parabolic equation corresponding to the Euler–Lagrange equations of an associated variational problem, covering spatially inhomogeneous models such as with latitudedependent albedo. Sufficient conditions are provided for the existence of at least three steadystate solutions in the form of two local minima and one saddle, that is, of coexisting “cold”, “warm” and unstable “intermediate” climates. We also give an interpretation of minimizers as “typical” or “likely” solutions of timedependent and stochastic 1D EBMs.
We then examine connections between the value function, which represents the minimum value (across all temperature profiles) of the objective functional, regarded as a function of greenhouse gas concentration, and the global mean temperature (also as a function of greenhouse gas concentration, i.e. the bifurcation diagram).
Specifically, the global mean temperature varies continuously as long as there is a unique minimizing temperature profile, but coexisting minimizers must have different global mean temperatures. Furthermore, global mean temperature is nondecreasing with respect to greenhouse gas concentration, and its jumps must necessarily be upward.
Applicability of our findings to more general spatially heterogeneous reaction–diffusion models is also discussed, as are physical interpretations of our results.
 Article
(1278 KB)  Fulltext XML

Supplement
(648 KB)  BibTeX
 EndNote
1.1 Lowdimensional energy balance models
Energy balance models are a fundamental tool used to understand the Earth's climate system and its energy dynamics. They represent the energy budget within the Earth's atmosphere, land, oceans and ice by quantifying the balance between incoming solar radiation and outgoing solar radiation. Although highly simplified compared to general circulation models, energy balance models (EBMs) are appreciated for their interpretability, mathematical tractability and ability to capture the essential dynamics of the Earth system (Budyko, 1969; Sellers, 1969; North, 1975; Ghil, 1976; Díaz, 1997; Cannarsa et al., 2022). Two important feedback mechanisms are typically present in such models: the icealbedo feedback and the Stefan–Boltzmann law. The positive icealbedo feedback occurs when the melting of ice and snow reduces the surface reflectivity (albedo), causing the planet to absorb more solar radiation. According to the Stefan–Boltzmann law, a warmer body emits more radiation, thereby providing a negative feedback which stabilizes the planet's temperature. Depending on the precise configuration, these mechanisms may endow EBMs with bistability, suggesting the existence of two stable climates commonly referred to as the snowball climate and the warm climate. The snowball climate, supported by palaeoclimatic evidence from the Cryogenian period around 650 million years ago, is characterized by the absence of vegetation and the presence of ice caps extending over the entire planet's surface. In contrast, the warm climate exhibits relatively low albedo, ice caps limited to the polar regions, and the presence of oceans and vegetation. Additionally, EBMs typically allow for a third possible climate, albeit unstable. Transitions between stable climates in an EBM, as well as in general multistable models, can occur in various ways. But two important mechanisms are the following. The first consists of changes in factors influencing the climate system, such as variations in greenhouse gas concentrations like carbon dioxide (CO_{2}), altering the balance of incoming and outgoing radiation and amplifying the greenhouse effect. Mathematically, this mechanism can be described by assuming that the model depends on one additional parameter, and changes in the parameter lead the model to undergo a bifurcation (Ashwin et al., 2012); the second consists of noiseinduced transitions resulting from unresolved processes in climate models or the representation of shorttimescale weather as stochastic forcing acting on slow variables, as observed in stochastic reduced models (Imkeller, 2001; Lucarini et al., 2022). These two types of transitions correspond to mechanisms recognized to induce climate tipping, that is, rapid nonlinear changes in the climate system with potentially irreversible and catastrophic consequences (Lenton et al., 2008, 2012; Scheffer et al., 2009; Lucarini and Bódai, 2019; Ghil and Lucarini, 2020).
A zerodimensional (0D) EBM is the simplest version of an EBM describing the evolution in time for the annual averaged global mean temperature T, without any space dependence (Berger, 1981; North, 1990; North and Kim, 2017; Ghil and Lucarini, 2020). This model is given by an ordinary differential equation (ODE) of the form
In this equation, C_{T}>0 represents the heat capacity, ${\stackrel{\mathrm{\u203e}}{Q}}_{\mathrm{0}}>\mathrm{0}$ is the globally averaged solar radiation and the coalbedo β is modelled by a continuous function (overbars typically denote globally averaged quantities). Further, q>0 is a positive parameter modelling the effect of the CO_{2} on the energy budget (Bastiaansen et al., 2022). The term ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{e}}\left(T\right)={\mathit{\sigma}}_{\mathrm{0}}{\mathit{\epsilon}}_{\mathrm{0}}{T}^{\mathrm{4}}$ on the righthand side of Eq. (1) accounts for the outgoing solar radiation, following the Stefan–Boltzmann law (where σ_{0} denotes the Stefan–Boltzmann constant, and ε_{0} is the globally averaged emissivity). The fixed points of the model are the solutions of the equation:
corresponding to points in Fig. 1 where the absorbed radiation ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{a}}\left(T\right)={\stackrel{\mathrm{\u203e}}{Q}}_{\mathrm{0}}\mathit{\beta}\left(T\right)+q$ and the emitted radiation ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{e}}\left(T\right)$ intersect. Figure 1a furthermore illustrates that this model is generally characterized by bistability, with two stable fixed points T_{S} and T_{W}. These points correspond to the snowball and warm climate states mentioned earlier and are separated by an unstable intermediate fixed point T_{M}. Furthermore, as highlighted by Fig. 1b, the stable points correspond to minimum points of a primitive function $\stackrel{\mathrm{\u203e}}{F}$ for the negative radiation budget $\stackrel{\mathrm{\u203e}}{R}$. In other words, $\stackrel{\mathrm{\u203e}}{F}$ is any regular function such that
To better capture the variability of global mean surface temperature, it has been proposed to add a stochastic forcing, such as white noise, to the radiation balance. This is interpreted as the effect of the fast components of the climate system, i.e. the weather, over slow components (Hasselmann, 1976; North and Cahalan, 1981; Imkeller, 2001; Díaz et al., 2009). For this reason, we are interested in considering the stochastic differential equation (SDE) given by
where ε>0 is the noise intensity and (W_{t})_{t≥0} is a Brownian motion (Baldi, 2017). This SDE is of gradient type and possesses a unique Gibbs invariant measure $\stackrel{\mathrm{\u203e}}{\mathit{\nu}}$ (Lelièvre and Stoltz, 2016). An invariant measure is a probability distribution $\stackrel{\mathrm{\u203e}}{\mathit{\nu}}$ in the state space of Eq. (2) (i.e. the real numbers in this case) with the property that if a solution T is distributed according to $\stackrel{\mathrm{\u203e}}{\mathit{\nu}}$ at some time t, then it remains so for all later times. It is given by
where Z is a normalization constant, and dT denotes the standard volume element on ℝ (we note the technical detail that to give meaning to Eqs. (2) and (3), the radiation budget $\stackrel{\mathrm{\u203e}}{R}$ should be extended to negative values for the Kelvin temperature T in a way such that $\stackrel{\mathrm{\u203e}}{F}\to +\mathrm{\infty}$ as $T\to \mathrm{\infty}$). The key observation from the explicit formula (Eq. 3) is that $\stackrel{\mathrm{\u203e}}{\mathit{\nu}}$ is concentrated around the minimum points of the function $\stackrel{\mathrm{\u203e}}{F}$. Indeed, if T_{0} is a strict minimum point and T_{1}≠T_{0} is a point close to T_{0} such that $\stackrel{\mathrm{\u203e}}{F}\left({T}_{\mathrm{1}}\right)>\stackrel{\mathrm{\u203e}}{F}\left({T}_{\mathrm{0}}\right)$, then the mass given by the measure ν in a small neighbourhood of T_{1} is exponentially lower than the mass around T_{0}; more specifically, the ratio between the two masses is given by $\mathrm{exp}\left(\frac{\mathrm{2}}{{\mathit{\epsilon}}^{\mathrm{2}}}\left(\stackrel{\mathrm{\u203e}}{F}\left({T}_{\mathrm{1}}\right)\stackrel{\mathrm{\u203e}}{F}\left({T}_{\mathrm{0}}\right)\right)\right)$.
A onedimensional (1D) EBM is given by a parabolic partial differential equation where the space variable is onedimensional (Budyko, 1969; Sellers, 1969; North and Kim, 2017). Denoting the temperature averaged in the zonal direction by $u=u(t,x)$, it extends the 0D EBM by introducing the sine of the latitude x=sin (ϕ), where $\mathit{\varphi}\in [\frac{\mathit{\pi}}{\mathrm{2}},\frac{\mathit{\pi}}{\mathrm{2}}]$ denotes the latitude and t≥0 represents time. We assume that the nonlinear radiation balance of the planet, denoted by $R(x,u;q)$, depends on the sine of the latitude and on an additive parameter q. This parameter models the effect of carbon dioxide concentration on the radiation budget (Bastiaansen et al., 2022). Atmospheric and ocean transport of heat between latitudes is modelled in a very simplified way by a diffusion term. Assuming spatially homogeneous diffusion in this introductory section and thus ignoring the dependence of κ on latitude and temperature, we obtain a nondegenerate reaction–diffusion equation:
where $\mathrm{\Delta}={\partial}_{xx}$ denotes the Laplace operator in dimension one, the Neumann boundary conditions impose noheat flux at the poles and $\stackrel{\mathrm{\u0303}}{u}$ is an initial condition. The steadystate solutions of this model, representing the asymptotic solutions for the timeevolving dynamics, correspond to the nonnegative solutions u=u(x) of the following elliptic problem:
where u=u(x) depends only on the space variable. This elliptic problem forms a necessary condition for u=u(x) to be our extremal (in particular a local minimizer) for the potential functional
where ${\partial}_{u}\mathcal{R}(x,u;q)=R(x,u;q)$, and $\mathrm{}\mathrm{}{u}^{\prime}\mathrm{}{\mathrm{}}_{\mathrm{2}}^{\mathrm{2}}={\int}_{\mathrm{1}}^{\mathrm{1}}{\left({u}^{\prime}\right)}^{\mathrm{2}}\mathrm{d}x$ is the square of the norm of u^{′} in ${L}^{\mathrm{2}}(\mathrm{1},\mathrm{1})$. Calculus of variations is a widely employed technique for studying the existence of a solution to the previous problem (North, 1975; North et al., 1979, 1981; Brezis, 2011). However, proving the existence of a local (but not global) minimum point is generally challenging, and this technique focuses on studying the existence of the global minimum point. The functional F_{q} in Eq. (6) has another interpretation though which renders it more important than being merely a characterization of solutions to the elliptic problem. Indeed, consider the stochastic partial differential equation (SPDE) on the Hilbert space $H={L}^{\mathrm{2}}(\mathrm{1},\mathrm{1})$, given by
obtained by adding a space–time white noise (W_{t})_{t≥0} modelled by a cylindrical Brownian motion on $H={L}^{\mathrm{2}}(\mathrm{1},\mathrm{1})$ to Eq. (7). R has a cutoff at negative temperature as in Sect. 3.1, and ε>0 is the noise intensity. We refer to Da Prato and Zabczyk (2014) for more details about SPDEs. It can be shown that this SPDE has a unique invariant Gibbs measure ν (Da Prato, 2004), given (broadly speaking) by an expression as in Eq. (3), with F_{q} replacing $\stackrel{\mathrm{\u203e}}{F}$ (see Sects. 2.1 and 3.1). Therefore, as in the zerodimensional case, ν concentrates on minimum points of the functional F_{q}. These minimizers satisfy the elliptic problem (Eq. 5), which therefore describes temperature profiles around which the solutions of the stochastic problem (Eq. 7) tend to cluster.
1.2 Main results and structure of the paper
This paper focuses on the study of the properties and the interpretation of the steadystate solutions of a 1D EBM depending on a bifurcation parameter. Motivated by 0D EBMs, there is a wide consensus in the literature, supported mainly by numerical simulations, regarding the existence of either one or three “interesting” steadystate solutions for 1D EBMs. Firstly, in Theorem 1 we prove the existence of a steadystate solution for the 1D EBM by solving the associated variational problem
i.e. showing the existence of a global minimum point for the functional F_{q} over a suitable function space 𝕏. Secondly, in Theorem 2 we provide sufficient conditions to have at least three steadystate solutions. These consist of two local minima and one saddle point. The conditions can be summarized as follows:
 i.
the viscosity κ should be sufficiently large,
 ii.
the spaceaveraged global radiation balance ℛ of the 1D EBM should present a doublewell potential with sufficiently deep minimum values attained at the two minimum points.
In essence, the conditions require that steadystate solutions of the spatially inhomogeneous 1D EBM (Eq. 5) are sufficiently well approximated by steadystate solutions of the spatially homogeneous model obtained from spatially averaging the terms in Eq. (5). These assumptions give us the possibility to prove the existence of two minimum points for F_{q}; further, these minimum points are also close to the minimum points of the spaceaveraged model. Then, the mountain pass theorem, a classical result from the calculus of variations, enables us to deduce the existence of a third steadystate solution (Ghil and Childress, 1987; Jabri, 2003). Thirdly, we investigate the uniqueness of the solution of the variational problem in terms of the value function
which is the minimum value attained by F_{q} as a function of q which relates to the greenhouse gas concentration.
 i.
In Theorem 3, we show that V is Lipschitz continuous, thus differentiable except for a Lebesgue zeromeasure set;
 ii.
Furthermore, the value function fails to be differentiable if and only if there are two or more co–existing global minimizers for F_{q}. Moreover, V is concave and hence its derivative V^{′} is nonincreasing.
 iii.
In Corollary 4, we demonstrate how the derivative of V is, up to the sign, the global mean temperature of the global minimum point u_{0} for F_{q}. This establishes a onetoone correspondence between the graph of V and the branch of the bifurcation diagram corresponding to u_{0}.
A byproduct of our analysis is that the global mean temperature is nonincreasing with respect to greenhouse gas concentration q. Moreover, it varies continuously with respect to q, as long as there is a unique minimizing temperature profile. However, for greenhouse gas concentration with coexisting global minimizers, the global mean temperature may exhibit as a discontinuous jump as coexisting minimizers cannot all have the same global mean temperature. Furthermore, if a jump occurs, it must necessarily be upward with increasing greenhouse gas concentration.
Our results have a number of interesting physical interpretations. The elliptic 1D EBM not only describes stationary solutions of the timedependent 1D EBM but moreover characterizes “likely” climates around which the solutions of the stochastic 1D EBM fluctuate. Global minimizers carry special importance as they are exponentially more likely than just local minimizers. Coexistence of global minimizers is just of special interest as these represent equally likely climate scenarios, and intuitively it seems plausible that rapid transitions between those climates are a dominant feature of the dynamics, although this point is not investigated further here. Furthermore coexistence of global minimizers implies a discontinuous change of global mean temperature which will jump upwards with increasing greenhouse gas concentration.
We expect that additional interesting physical conclusions can be drawn through identifying F_{q} with the (negative) entropy production rate (North and Kim, 2017, Sect. 7.4.2); this will be subject to future work.
This paper is organized as follows. In Sect. 2, we describe the methodology used throughout our work. Firstly, we review the 1D EBM proposed in Bastiaansen et al. (2022). This model serves as the reference for our paper, and it is characterized by the presence of an additive parameter in the radiation budget, which determines the number of steadystate solutions. Secondly, we recall the properties of the steadystate solutions of the 1D EBM that can be obtained from numerical simulations. Finally, we rigorously define the stochastic EBM by introducing space–time white noise. Specifically, we review the invariant measure formula for the resulting reaction–diffusion SPDE. In Sect. 3, we present our novel findings. In Sect. 3.1, we discuss the existence of a solution for the variational problem and outline the properties of the potential functional. Moreover, we explain why the invariant measure of the stochastic EBM concentrates around the global minimum points of the potential functional. Finally, we provide sufficient conditions to demonstrate the existence of at least three steadystate solutions. In Sect. 3.2, we characterize the uniqueness of the solution to the variational problem in terms of the value function. Additionally, we demonstrate that the value function is Lipschitz concave and that its derivative is nonincreasing. In Sect. 3.3, we illustrate how knowledge of the value function allows derivation of a portion of the bifurcation diagram and vice versa. In Sect. 4, we offer a comprehensive summary of our work. In Appendix A, we describe the finite difference method employed to conduct the numerical simulations presented in this study. Furthermore, the Supplement includes rigorous proofs of our main results.
2.1 A 1D energy balance model
The fundamental mechanism of 1D EBMs is that the temperature $u=u(t,x)$, averaged in the zonal direction, evolves in time due to (i) the diffusion of energy between adjacent regions, (ii) the energy absorbed by the planet and (iii) the energy emitted by the planet. The 1D EBM we consider in this paper is a Sellertype EBM where the absorbed radiation depends on an additive parameter (Bastiaansen et al., 2022). We only add a change in the diffusion term in order to get a nondegenerate parabolic PDE. Given an initial condition $\stackrel{\mathrm{\u0303}}{u}$, the nonlinear, parabolic, reaction–diffusion PDE governing the model is given by
where R_{a} and R_{e} represent the radiation absorbed and emitted by the planet per unit area, respectively. C_{T} is the heat capacity, and the differential term parameterizes the meridional heat transport. The boundary conditions impose no flux at the poles. We now provide further details regarding the parameterization of these terms. The values of the constants of the model can be found in Table 1.
Firstly, the absorbed radiation is assumed to have the following form:
where Q_{0} is the solar radiation per unit area, and α is the albedo. The solar radiation is assumed to be
where ${\widehat{Q}}_{\mathrm{0}}$ is the mean solar radiation, and c_{i} represents constants. The albedo, which is the proportion of the incident light or radiation that is reflected by a surface, is parameterized by a smooth monotonically decreasing function with a peak derivative in a reference temperature u_{ref} close to the melting point of ice. Specifically,
where K>0 is a rate parameter and α_{1}>α_{2} are respectively the ice albedo and the water albedo.
Second, the emitted radiation is modelled using the Stefan–Boltzmann law, in other words assuming that the Earth radiates as a black body. Under this assumption, the energy radiated is proportional to the fourth power of its temperature and it is given by
where ε_{0} and σ_{0} are respectively the emissivity and Boltzmann's constant. The additive parameter q describes, in a simplified way, the radiative forcing by CO_{2}, i.e. the effect of atmospheric CO_{2} on the energy budget (IPCC, 2014). It is worth explaining (i) the additive structure of q and (ii) its independence on the spatial variable x. About the first point, denote by C the global CO_{2} concentration in part per million (ppm) and assume, just for explanation purposes, a dependence of the outgoing radiation both on u and C. To avoid confusion, we denote this by ${\widehat{R}}_{\mathrm{e}}={\widehat{R}}_{\mathrm{e}}(u,C)$, the outgoing radiation depending on temperature and CO_{2} concentration. If we linearize ${\widehat{R}}_{\mathrm{e}}$ with respect to temperature, we get
where $\widehat{u}$ is a reference temperature. In Myhre et al. (1998), three radiative transfer models are used in order to get that the dependence of the outgoing radiation with respect to changes in CO_{2} is given by
where C_{0} is a reference CO_{2} concentration, and ${A}_{\mathrm{1}},{A}_{\mathrm{2}}>\mathrm{0}$ are explicit constants. In conclusion,
and thus the radiative forcing of CO_{2} has an additive structure. About point (ii), we adopt the wellmixed hypothesis for CO_{2}. In other words, we assume that atmospheric CO_{2} is globally homogeneous, thereby inducing a radiative forcing q independent of latitude (IPCC, 2001). This assumption overlooks the spatial pattern of CO_{2} concentration, which affects many aspects of the climate system, such as the poleward heat transport (Huang et al., 2017). It was the stateoftheart assumption 2 decades ago, although today it is common to keep in consideration the spatial distribution of radiative forcing (IPCC, 2001; Byrne and Goldblatt, 2014; Zhang et al., 2019).
The third component of the model is the term ∂_{x}(κ(x)u_{x}). It parameterizes the meridional heat transport, that is, the phenomenon resulting from the poleward transportation of heat by the Earth–atmosphere system due to the surplus of net radiation heating in the tropics and the deficit in the poleward regions. Usually, the diffusion function κ(x) is assumed null at the poles, i.e. with a form such as $\mathit{\kappa}\left(x\right)=D(\mathrm{1}{x}^{\mathrm{2}})$, where D is a diffusion constant. This choice is based on the paradigm of mimicking the conduction of heat on a sphere; see North and Kim (2017) for a derivation. On the other hand, it leads to mathematical difficulties in the treatment of the singular PDE arising, in particular in the study of the corresponding variational problem, from which all our results follow. To avoid these difficulties, which at the moment remain an open problem to solve, we add as a simplifying assumption that κ is nondegenerate and given by
We choose δ=0.003, but its value is not important for the results of this work, and different choices can be made.
For the parabolic problem (Eq. 8), the global existence and uniqueness of the solution can be demonstrated, given a regular initial condition (Temam, 1997). Furthermore, if the initial condition is nonnegative, the solution remains nonnegative for any time t>0. This can be shown proving that $[\mathrm{0},+\mathrm{\infty})$ is an invariant region for Eq. (8), exploiting the fact that there exist ${C}_{\mathrm{1}},{C}_{\mathrm{2}}>\mathrm{0}$ such that $R(x,u;q)>{C}_{\mathrm{1}}>\mathrm{0}$ for all $x\in [\mathrm{1},\mathrm{1}]$ and $u\in [\mathrm{0},{C}_{\mathrm{2}}]$ (Smoller, 2012).
We recall the formulation of stochastic EBMs using the theory of SPDEs (Da Prato and Zabczyk, 2014). Denote by Δ the Laplace operator with Neumann boundary conditions. Given an initial condition $\stackrel{\mathrm{\u0303}}{u}\in H$, we consider the following SPDE:
where ε>0, and (W_{t})_{t≥0} is a cylindrical Brownian motion on H. Under the minor cutoff modifications introduced in Sect. 3.1, it can be proved that the Hvalued stochastic process (u_{t})_{t} which solves in a mild sense (9) is unique and has continuous trajectories (Da Prato and Zabczyk, 2014). In addition to this there exists a unique Gibbs invariant measure
where ℛ is as in Eq. (6), and $\mathit{\mu}\sim \mathcal{N}(\mathrm{0},\frac{{\mathit{\epsilon}}^{\mathrm{2}}}{\mathrm{2}\mathit{\kappa}}{\mathrm{\Delta}}^{\mathrm{1}})$ is a symmetric Gaussian measure on H with covariance $\mathcal{Q}=\frac{{\mathit{\epsilon}}^{\mathrm{2}}}{\mathrm{2}\mathit{\kappa}}{\mathrm{\Delta}}^{\mathrm{1}}$ (Da Prato, 2004, 2006). The covariance operator $\mathcal{Q}:H\to H$ is the unique linear operator such that ${\int}_{H}\langle {h}_{\mathrm{1}},\mathit{\varphi}\rangle \langle {h}_{\mathrm{2}},\mathit{\varphi}\rangle \mathit{\mu}\left(\mathrm{d}\mathit{\varphi}\right)$ for each ${h}_{\mathrm{1}},{h}_{\mathrm{2}}\in H$, where $\langle \cdot ,\cdot \rangle $ denotes the scalar product in H. Further, it can be shown that 𝒬 is symmetric and positivedefinite, and its eigenvalues (λ_{k})_{k∈ℤ} satisfy ${\sum}_{k\in \mathbb{Z}}{\mathit{\lambda}}_{k}<\mathrm{\infty}$. In the following lines, given a symmetric, positivedefinite operator 𝒬 such that ${\sum}_{k\in \mathbb{Z}}{\mathit{\lambda}}_{k}<+\mathrm{\infty}$, we are going to explain how to construct an Hvalued random variable X with law 𝒩(0,𝒬). Indeed, consider a sequence (R_{k})_{k∈ℤ} of independent and identically distributed ℝ𝒩(0,1) random variables defined on a probability space $(\mathrm{\Omega},\mathcal{F},\mathbb{P})$. We can assume without loss of generality that the eigenvectors (e_{k})_{k∈ℤ} associated with the eigenvalues (λ_{k})_{k∈ℤ} form an orthonormal basis of H. Then, the Hvalued random variable
is well defined, i.e. the series defining X converges in ${L}^{\mathrm{2}}(\mathrm{\Omega},\mathcal{F},\mathbb{P};H)$ and has law 𝒩(0,𝒬) (Da Prato, 2006, Proposition 2.18). Further, the convergence also holds ℙ almost surely in H (Da Prato, 2006, Proposition 2.13).
As mentioned in the introduction, this measure is concentrated on minimum points of the functional F_{q}. A heuristic explanation of this fact can be found in Sect. 3.1.
The stationary problem associated with the 1D EBM is given by the elliptic equation for u=u(x):
These solutions can be either stable or unstable, depending on the longterm behaviour of their infinitesimal perturbations. As pointed out in Bastiaansen et al. (2022), if the reaction–diffusion equation for $u=u(t,x)$ was spacehomogeneous, i.e. of the form
then the stable steadystate solutions would correspond to functions that are constant in space and time, with values given as the roots of
A rigorous result in this direction has been shown in Gaspar and Guaraco (2018). Indeed, for a fixed doublewell symmetric potential, it has been proved that (i) if κ is large enough, the only steadystate solutions of Eq. (12) are the constants where the potential is critical, and (ii) the number of unstable steadystate solutions to Eq. (12) can be made arbitrarily large as κ→0. Introducing a spatial dependence in $R=R(x,u)$ leads to a spaceheterogeneous model. Depending on the space heterogeneity, it can exhibit any number of both stable and unstable steadystate solutions (Bastiaansen et al., 2022). The variational approach to the study of steadystate solutions provides a tool for characterizing the stable ones, which are the local minimum points of a functional.
In the following paragraph, we describe the properties of the solutions of Eq. (11). As the parameter q changes, numerical simulations for Eq. (11) suggest the existence of either one or three steadystate solutions. That is, there exists q_{1}<q_{2} s.t. Eq. (8) has one steadystate solution if q<q_{1} or q>q_{2}, and the steadystate solutions are three if ${q}_{\mathrm{1}}\le q\le {q}_{\mathrm{2}}$. In the latter case, we denote the solutions by ${u}_{\mathrm{S}}\le {u}_{\mathrm{M}}\le {u}_{\mathrm{W}}$, corresponding respectively to the snowball climate, a middle (or intermediate) climate and the warm climate. As an analogy, we denote by u_{S} the unique steadystate solution for q<q_{1} and by u_{W} the unique one for q>q_{2}. Figure 2a shows the bifurcation diagram of the model in the $(q,\stackrel{\mathrm{\u203e}}{u})$ plane, where $\stackrel{\mathrm{\u203e}}{u}={\int}_{\mathrm{1}}^{\mathrm{1}}u\left(x\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}x$ denotes the average temperature. Figure 2b depicts the three steadystate solutions for $q=\mathrm{25}\in ({q}_{\mathrm{1}},{q}_{\mathrm{2}})$.
A stability analysis can be conducted to determine the stability of the steadystate solutions. The results show that u_{S} and u_{W} are stable, while the middle climate, u_{M}, is unstable. Furthermore, it is worth noting that special values $q={q}_{\mathrm{1}},{q}_{\mathrm{2}}$ correspond to bifurcation points of saddlenode type, where the unstable solution u_{M} collides with either u_{W} (for q=q_{1}) or u_{S} (for q=q_{2}) and then disappears. These numerical findings regarding the number and stability of the steadystate solutions will be supported and validated using rigorous arguments, as in the next section.
3.1 Potential functional and its minimizer
In this section, we (i) provide an intuitive motivation for why the invariant measure for the stochastic EBM concentrates on minimum points of the functional F_{q}, (ii) prove the existence of global minimum points for F_{q} using the direct method, and (iii) present sufficient conditions on the viscosity κ and the spaceaveraged potential $\stackrel{\mathrm{\u203e}}{\mathcal{R}}\left(u\right)=\frac{\mathrm{1}}{\mathrm{2}}{\int}_{\mathrm{1}}^{\mathrm{1}}\mathcal{R}(x,u)\mathrm{d}x$, with ${\partial}_{u}\mathcal{R}=R={R}_{\mathrm{e}}{R}_{\mathrm{a}}$, to ensure that the 1D EBM has at least three steadystate solutions.
Firstly, consider the stochastic EBM (9). Assume that for a negative value of u, where the model has no physical meaning, the Stefan–Boltzmann law is extended as
And β is smoothly extended to $\stackrel{\mathrm{\u0303}}{\mathit{\beta}}$ by setting it to zero outside the physically relevant range, as described in the Supplement. Then, Eq. (9) possesses a unique Gibbs invariant probability measure given by
where $(u{)}_{+}=max\mathit{\{}u,\mathrm{0}\mathit{\}}$ is the positive part, $\mathcal{N}(\mathrm{0},\frac{{\mathit{\epsilon}}^{\mathrm{2}}}{\mathrm{2}}{\mathrm{\Delta}}^{\mathrm{1}})$ denotes a symmetric Gaussian measure with covariance operator $\mathcal{Q}=\frac{{\mathit{\epsilon}}^{\mathrm{2}}}{\mathrm{2}\mathit{\kappa}}{\mathrm{\Delta}}^{\mathrm{1}}$ over the Hilbert space $H={L}^{\mathrm{2}}(\mathrm{1},\mathrm{1})$ and Z is the normalization constant. See Da Prato (2004) for a rigorous derivation of the invariant measure for a reaction–diffusion model with a polynomial homogeneous reaction term. We move to explain in what sense ν is concentrated around minimum points of F_{q}. In fact, for u∈H the Gaussian measure μ is formally given by
where ${\mathcal{Q}}^{\mathrm{1}}=\frac{\mathrm{2}\mathit{\kappa}}{{\mathit{\epsilon}}^{\mathrm{2}}}\mathrm{\Delta}.$ Here, Z_{1} is a normalization constant, $\langle \cdot ,\cdot \rangle $ denotes the scalar product in H and du is a formal notation for the Lebesgue measure on H. If we perform an integration by parts, we get
Plugging the previous identity into Eq. (10), we obtain
From this heuristic formula, we see that points u such that F_{q}(u) is not a global minimum have exponentially smaller density than the minimum points. Indeed, if u_{1} is a global minimum point and u≠u_{1}, then the mass given by ν in a small neighbourhood around u is exponentially smaller than the mass given to a neighbourhood of the same size around u_{1}; in particular, the ratio between the two masses is given by $\mathrm{exp}\left(\frac{\mathrm{2}}{{\mathit{\epsilon}}^{\mathrm{2}}}\left({F}_{q}\left(u\right){F}_{q}\left({u}_{\mathrm{1}}\right)\right)\right)$. The previous derivation is formal, because the Lebesgue measure cannot be defined on an infinitedimensional Hilbert space. For a more rigorous explanation, see Sect. S2 in the Supplement.
Next, we discuss the properties of the functional ${F}_{q}:{H}^{\mathrm{1},\mathrm{2}}(\mathrm{1},\mathrm{1})\cap \left\{u\ge \mathrm{0}\right\}\to \mathbb{R}$ given by
where B is a primitive of the coalbedo $\mathit{\beta}\left(u\right)=\mathrm{1}\mathit{\alpha}\left(u\right)$, and ${H}^{\mathrm{1}}={H}^{\mathrm{1},\mathrm{2}}(\mathrm{1},\mathrm{1})$ denotes the Sobolev space of order 1 and exponent 2, i.e. the function space where a function u and its derivative u^{′} (in a weak sense) are both square integrable over $[\mathrm{1},\mathrm{1}]$. See Brezis (2011) for more details about Sobolev spaces. The functional F_{q}, depending on the parameter q, is known in the literature as a potential functional or Lyapunov function (North et al., 1979; North and Kim, 2017). The study of the functional F_{q} gives useful information thanks to its links with the invariant measure for the stochastic 1D EBM, as we have seen, and the stable steadystate solutions for the deterministic 1D EBM which emerge as necessary conditions for the stationarity of F_{q}. Going deeper with the former point, the first variation of F_{q} in the point u in direction h is given by
where in the last identity we have used integration by parts. Since h is arbitrary, u is a stationary point for the functional F_{q} if and only if it is a steadystate solution for the EBM. In particular, local extremum points for F_{q} correspond to steadystate solutions of the EBM. Any local minimizer of F_{q} represents a locally attractive solution of the deterministic 1D EBM. In view of our interpretation of F_{q} in terms of the invariant measure, however, global minimizers play a special role since if present and unique they are exponentially more likely than any other state (including minimizers that are just local). The following result establishes the existence of a global minimum point for F_{q}.
If q>0, then there exists a global regular nonnegative minimizer for F_{q}. In other words, if we consider the variational problem
then there exists u_{0}∈C^{∞} such that u_{0} is a solution of the EBM and
In addition to this, if q belongs to a bounded interval, then u can be bounded uniformly with respect to q:
A rigorous proof of the previous result can be found in Sect. S3. The proof relies on standard arguments from the direct method of calculus of variation, exploiting the fact that the outgoing radiation in the EBM model prevents the temperature from being too high.
Concerning the existence of two local minimum points, let us describe a sufficient condition. Consider the potential function $\stackrel{\mathrm{\u203e}}{\mathcal{R}}:\mathbb{R}\to \mathbb{R}$ coming from the spaceaveraged model
If the viscosity κ>0 is sufficiently large and the function $\stackrel{\mathrm{\u203e}}{\mathcal{R}}$ has a doublewell shape with sufficiently deep minimum values attained at the minimum points, then we are able to prove the existence of two minimum points for F_{q}. Further, it is possible to prove that the functional F_{q} satisfies a compactness condition known as the Palais–Smale condition. This property and the mountain pass theorem give the possibility to deduce the existence of a third steadystate solution. Next, we characterize a situation in which there are three steadystate solutions, two of which are local minimizers (Jabri, 2003). This is summarized in the following result.
Denote by ${B}_{{H}^{\mathrm{1}}}\mathrm{(}v\mathrm{,}\mathit{\rho}\mathrm{)}\mathrm{=}\mathit{\{}u\mathrm{\in}{H}^{\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\mathrm{\mid}\phantom{\rule{0.25em}{0ex}}\mathrm{}\mathrm{}u\mathrm{}v\mathrm{}{\mathrm{}}_{{H}^{\mathrm{1}}}\mathrm{<}\mathit{\rho}\mathit{\}}$ the open ball in H^{1} with centre v and radius ρ>0. Assume $\stackrel{\mathrm{\u203e}}{\mathrm{R}}$ has two nonnegative minimum points u_{1}≠u_{2}, with F_{q}(u_{1})≥F_{q}(u_{2}). Then, there exist ω>0 and $f\mathrm{,}g\mathrm{\in}O\mathrm{\left(}{\mathit{\epsilon}}^{\mathrm{}\mathrm{1}}\mathrm{\right)}$ as $\mathit{\epsilon}\mathrm{\to}{\mathrm{0}}^{\mathrm{+}}$ such that if $\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}\mathrm{>}\mathrm{0}$ satisfies
 i.
${\stackrel{\mathrm{\u203e}}{\mathcal{R}}}^{\prime \prime}\left({u}_{i}\right)>f\left(\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}\right)$, for $i=\mathrm{1},\mathrm{2}$;
 ii.
$\mathit{\kappa}>g\left(\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}\right)$;
 iii.
$\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}\le \mathit{\omega}$;
then ${\stackrel{\mathrm{\u0303}}{F}}_{q}$ has two local minimum points, ${\stackrel{\mathrm{\u0303}}{u}}_{\mathrm{1}},{\stackrel{\mathrm{\u0303}}{u}}_{\mathrm{2}}$, such that
 a.
${B}_{{H}^{\mathrm{1}}}({u}_{\mathrm{1}},\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}})\cap {B}_{{H}^{\mathrm{1}}}({u}_{\mathrm{2}},\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}})=\mathrm{\varnothing}$;
 b.
${\stackrel{\mathrm{\u0303}}{u}}_{i}\in {B}_{{H}^{\mathrm{1}}}({u}_{i},\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}),$ for $i=\mathrm{1},\mathrm{2}$.
 c.
If $\mathrm{}\mathrm{}u{u}_{\mathrm{1}}\mathrm{}{\mathrm{}}_{{H}^{\mathrm{1}}}=\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}$, then ${F}_{q}\left(u\right)\ge {F}_{q}\left({u}_{\mathrm{1}}\right)+\mathit{\delta}$, with $\mathit{\delta}=\mathit{\delta}\left(\stackrel{\mathrm{\u203e}}{\mathit{\epsilon}}\right)>\mathrm{0}.$
Note how the previous result can be also interpreted as giving sufficient conditions for the convergence of the stable solutions of a spaceinhomogeneous EBM to the stable solution of the corresponding spaceaveraged model, as the diffusion becomes large.
3.2 Value function and uniqueness for the functional minimizer
The key element of this section is the value function, which is given by
From Sect. 3.1, we know that the previous infimum is indeed a minimum, and so V(q) can be interpreted as the minimum possible value attained by the potential functional over the possible temperature profiles u. Since a minimum point for F_{q} is also a stationary point for the functional, the value function can be evaluated numerically by computing the minimum of the three steadystate solutions u_{S},u_{M} and u_{W}. Following this strategy, Fig. 3 shows $q\mapsto {F}_{q}\left({u}_{*}\right)$, with ${u}_{*}\in \left\{{u}_{\mathrm{S}},{u}_{\mathrm{M}},{u}_{\mathrm{W}}\right\}$. Particularly, there exists a point q_{3} such that u_{S} is the global minimum point of F_{q} for q<q_{3}, while u_{W} is the global minimum point for q>q_{3}. Further, for q=q_{3} the function F_{q} has two different global minimum points u_{S} and u_{W}, and q=q_{3} corresponds a nondifferentiability point for V. In addition to this, the value function appears to be concave, thus with a decreasing derivative, where it exists.
Summarizing, the numerical evaluations of V(q) suggest the following result, a rigorous proof of which is included in the Supplement.
Assume q belongs to a bounded interval. Then
 i.
V is Lipschitz continuous.
 ii.
q is a nondifferentiable point for V if and only if there is more than one minimizer for F_{q}.
 iii.
V is concave and V^{′} is nonincreasing.
We also see numerically that u_{M} is actually never a global minimizer for the specific functional F_{q} considered here, but we do no have rigorous proof of this fact. Let us briefly discuss the proofs of the previous points. The proof of (i) follows from the facts that the sup norm of the minimizer u_{0} can be bounded uniformly in q and that, given a family {g_{i}}_{i∈I} of L_{i} Lipschitz functions g_{i}, the $\underset{i\in I}{inf}\phantom{\rule{0.25em}{0ex}}{g}_{i}$ is Lipschitz continuous if the constants L_{i} can be bounded uniformly. In our case, given u∈H^{1} is nonnegative, we have
where M>0 is the constant appearing in Eq. (15). On the other hand, the proof of point (ii) is less straightforward, although being very similar to the one for the existence of a solution for the variational problem. The proof of point (iii) makes use of the concept of semiconcavity, a generalization of that of concavity, which is fundamental in optimal control (Cannarsa and Sinestrari, 2004). The main reason for the concavity of V though is that V is an infimum over functions that are affine in q. Hence, the fact that q is additive is essential for this result. More details can be found in Sects. S4 and S5 in the Supplement.
3.3 Value function graph and bifurcation diagram
An additional property of the value function can be observed when comparing the bifurcation diagram (Fig. 4a) and the graph of the value function (Fig. 3).
Corollary 4. If V is differentiable, then
${V}^{\prime}\left(q\right)={\int}_{\mathrm{1}}^{\mathrm{1}}{u}_{\mathrm{0}}\left(x\right)\mathrm{d}x$, where u_{0} is the only minimizer for F_{q}.
In other words, the part of the bifurcation diagram that corresponds to the global minimizer, represented by the subgraph $(q,{\int}_{\mathrm{1}}^{\mathrm{1}}{u}_{\mathrm{0}}(x),\mathrm{d}x)$, can be determined based on the knowledge of V^{′}, and vice versa. Figure 4 compares Figs. 2a and 3, highlighting in magenta the corresponding parts of the two graphs. From the mathematical point of view, the previous result is a consequence of the proof of Theorem 3.
It is worth pointing out that by combining Theorem 3 and Corollary 4, a valuable property emerges, i.e. the global mean temperature of the functional minimizer is nondecreasing with respect to q. In other words, as the concentration of CO_{2} rises, the global mean temperature increases. Additionally, through this monotonicity and Froda's theorem (Rudin, 1976, Theorem 4.30), we also establish that the global mean temperature is continuous, except for, at most, a countable number of upward jumps.
In the second part of this section, we demonstrate the applicability of Corollary 4 to other reaction–diffusion equations. We use as an example a spatially heterogeneous Allen–Cahn equation (ACE), already considered in Bastiaansen et al. (2022). For an initial condition $\stackrel{\mathrm{\u0303}}{u}$, this model is given by
The associated elliptic problem for u=u(x) is
In this case, the potential functional takes the form
and all the properties discussed in Sect. 3.1 and 3.2 can be extended to this equation. Specifically, Theorems 1, 2 and 3 hold. But in this case, the structure of the bifurcation diagram is more complex, even if symmetric with respect to q=0. Indeed, through numerical experiments, it is possible to deduce the existence of $\mathrm{0}<{q}_{\mathrm{4}}<{q}_{\mathrm{5}}$ such that (a) for q>q_{5} or q<q_{4}, there exists a single steadystate solution, which is stable, (b) for ${q}_{\mathrm{4}}<\mathrm{}q\mathrm{}<{q}_{\mathrm{5}}$ there are three steadystate solutions, two of which are stable while the third is unstable. Further, $q={q}_{\mathrm{4}},{q}_{\mathrm{5}}$ are bifurcation points of the saddlenode type. We denote by u_{1} the steadystate solution for $q<{q}_{\mathrm{5}}$, by u_{2} and u_{3} the steadystate solutions appearing at the bifurcation point $q={q}_{\mathrm{5}}$ and existing for ${q}_{\mathrm{5}}<q<{q}_{\mathrm{4}}$ in addition to u_{1}, and by u_{4} and u_{5} the steadystate solutions appearing at q=q_{4} and existing for ${q}_{\mathrm{4}}<q<{q}_{\mathrm{5}}$ in addition to u_{3}. Regarding the potential functional J_{q}, in this case there exists ${q}_{\mathrm{6}}\in ({q}_{\mathrm{4}},{q}_{\mathrm{5}})$ such that u_{1} is the global minimum point for the functional for $q<{q}_{\mathrm{6}}$, and u_{3} is the global minimum point for ${q}_{\mathrm{6}}<q<{q}_{\mathrm{6}}$, while u_{5} becomes the global minimum point for q>q_{6}. A picture for the bifurcation diagram just described and the value function is shown in Fig. 5. Note that $q=\pm {q}_{\mathrm{6}}$ represents the only values of the parameter q for which the value function is not differentiable and also the only points in which the global minimizer of the variational problem is not unique.
In this paper, we have considered a onedimensional energy balance model depending on a bifurcation parameter q, describing the effect of CO_{2} concentration in the atmosphere and affecting the energy absorbed by the planet. Numerical simulations show that this model can exhibit either one or three asymptotic solutions, depending on the values of q. We began our analysis by introducing the potential functional F_{q} associated with the steadystate solutions. The functional F_{q} has significant implications, as it is closely linked to both the stability of steadystate solutions of the EBM and the invariant measure for the stochastic EBM obtained by perturbing the model with an additive Gaussian white noise. In particular, the invariant measure of the system concentrates on global minimizers of F_{q}, giving them exponentially larger weight than local minimizers. By analysing the first variation of F_{q} and applying standard arguments from the direct method of calculus of variations, we established that F_{q} possesses a global regular minimizer for all values of the parameter q. Furthermore, we provide sufficient conditions to prove the existence of at least three steadystate solutions for the 1D EBM.
We then introduced the value function V(q), which represents the minimum value attained by the potential functional among all possible temperature profiles. By evaluating V(q) numerically using the steadystate solutions u_{S},u_{M} and u_{W}, we observed that the function exhibits Lipschitz continuity and concavity. Furthermore, nondifferentiability points of V(q) coincide with points where multiple global minimizers exist for F_{q}. Lastly, when V is differentiable, its derivative is nonincreasing and equal to the negative global mean temperature, i.e. ${V}^{\prime}\left(q\right)={\int}_{\mathrm{1}}^{\mathrm{1}}{u}_{\mathrm{0}}\left(x\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}x$, where u_{0} is the minimizer for F_{q}. Moreover, as a consequence of the explicit expression for V^{′}, the global mean temperature is nondecreasing with respect to q, and it is continuous, except for a Lebesgue zeromeasure set of upwards jumps. These are the nondifferentiability points of V, corresponding to the case where two or more global minimizers, hence multiple climates equally probable, exist for the stochastic EBM. These findings, which we are able to prove rigorously, allow us to establish a correspondence between the bifurcation diagram and the graph of the value function. Additionally, we applied our results to a spatially inhomogeneous Allen–Cahn equation to show how our results still hold for more general spaceinhomogeneous reaction–diffusion equations.
The diffusion function κ=κ(x) that we have examined is nondegenerate at the boundary of the spatial domain. As noted in Sect. 2.1, this is an assumption to simplify the study of the variation problem. At present, there remains a problem with how to extend our results to the case where κ is degenerate at the boundary.
Regarding the impact of our work on current climate change, we have characterized climate as an invariant measure within a stochastic equation that describes temperature. The emission of CO_{2} is considered a parameter influencing the shape of this invariant measure, particularly in relation to the points around which the measure is concentrated. From our perspective, the climate we are currently witnessing reflects changes in the invariant measure, representing a realization of a random variable with that invariant measure as its distribution. Moreover, we have demonstrated the monotonic relationship between global mean temperature and CO_{2}. Finally, we have outlined simple conditions, adaptable to other multistable reaction–diffusion models, to establish the existence of three asymptotic climate states.
Concerning future development of this work, one interesting aspect we are working on is to understand how the invariant measure for the stochastic EBM changes close to bifurcation points. This points in the direction of using statistical indicators to detect the approach of tipping points, which in our model correspond to points of discontinuity of the global mean temperature with respect to the parameter q.
In this section, we describe the numerical method adopted to approximate the solutions of the elliptic problem (11) numerically. We used a classical finite difference scheme, which we are going to illustrate (Quarteroni and Valli, 2008; Thomas, 2013). To simplify the notation, let's define
as the nonlinear reaction term. We consider a uniform mesh for $[\mathrm{1},\mathrm{1}]$ made of n+1 points
Then, the solution to the problem can be approximated by considering the following system:
where u_{−1}and u_{n+1} are ghost points, and
The system of equations can be written in vector form as
with $\mathit{u}={\left[\begin{array}{c}{u}_{\mathrm{1}},\mathrm{\cdots}{u}_{n+\mathrm{1}}\end{array}\right]}^{T}$ and ${f}_{i}=f({x}_{i},u({x}_{i}\left)\right)$. At this point, multiplying the first equation by $\mathrm{2}\frac{{\mathit{\kappa}}_{\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta}x}$, subtracting the second one and dividing by 2, we get
In a symmetric way, multiplying the last equation by $\mathrm{2}\frac{{\mathit{\kappa}}_{n\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta}x}$, subtracting the second last equation and dividing by 2, we get
In this way, the Neumann version of the elliptic problem has the following form:
and consists of a set of (n+3) nonlinear equations, whose solution u can be approximated using the Newton–Raphson method (NRM). The initial guess used to start the iteration in NRM is obtained via a shooting method, thus reducing the boundary value problem given by the elliptic PDE in Eq. (11) to an initial value problem (IVP). A linear search is applied to find the shooting parameter, i.e. the initial condition of the IVP. Lastly, the solution of the IVP is approximated using the classical Euler's method for ODEs.
This work does not include any externally supplied code, data, or other material. All material in the text and figures was produced by the authors using standard mathematical and numerical analysis by the authors. The code is available at Zenodo (https://doi.org/10.5281/zenodo.10469451, Del Sarto et al., 2024).
The supplement related to this article is available online at: https://doi.org/10.5194/npg311372024supplement.
GDS conceptualized the paper, performed the numerical simulations, and took the lead role in writing and revising the paper. JB, FF and TK conceptualized the paper and supervised the writing. All authors provided critical feedback and helped shape the research.
The contact author has declared that none of the authors has any competing interests.
Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We are very grateful to the referee for carefully reading the paper and for their valuable comments. We acknowledge fruitful discussions with Valerio Lucarini and Robbin Bastiaansen. Gianmarco Del Sarto would like to thank the Department of Mathematics and Statistics, University of Reading, for its hospitality. Tobias Kuna and Jochen Bröcker would like to thank the Scuola Normale Superiore for its hospitality. Gianmarco Del Sarto is supported by the Italian national interuniversity PhD course in sustainable development and climate change.
Franco Flandoli's research is funded by the European Union (ERC, NoisyFluid, no. 101053472).
This paper was edited by Natale Alberto Carrassi and reviewed by two anonymous referees.
Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noiseinduced and ratedependent examples in the climate system, Philos. T. Roy. Soc. A, 370, 1166–1184, https://doi.org/10.1098/rsta.2011.0306, 2012. a
Baldi, P.: Stochastic Calculus, Springer International Publishing, https://doi.org/10.1007/9783319622262, 2017. a
Bastiaansen, R., Dijkstra, H. A., and von der Heydt, A. S.: Fragmented tipping in a spatially heterogeneous world, Environ. Res. Lett., 17, 045006, https://doi.org/10.1088/17489326/ac59a8, 2022. a, b, c, d, e, f, g
Berger, A. (Ed.): Climatic Variations and Variability: Facts and Theories, Springer, Netherlands, https://doi.org/10.1007/9789400985148, 1981. a
Brezis, H.: Functional analysis, Sobolev spaces and partial differential equations, vol. 2, Springer, https://doi.org/10.1007/9780387709147, 2011.Please provide persistent identifier (DOI preferred). a, b
Budyko, M. I.: The effect of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619, 1969. a, b
Byrne, B. and Goldblatt, C.: Radiative forcing at high concentrations of well‐mixed greenhouse gases, Geophys. Res. Lett., 41, 152–160, https://doi.org/10.1002/2013gl058456, 2014. a
Cannarsa, P. and Sinestrari, C.: Semiconcave Functions, Hamilton–Jacobi Equations, and Optimal Control, Birkhäuser Boston, ISBN 9780817644130, https://doi.org/10.1007/b138356, 2004. a
Cannarsa, P., Lucarini, V., Martinez, P., Urbani, C., and Vancostenoble, J.: Analysis of a twolayer energy balance model: long time behaviour and greenhouse effect, arXiv [preprint], https://doi.org/10.48550/arXiv.2211.15430, 2022. a
Da Prato, G.: Kolmogorov Equations for Stochastic PDEs, Birkhäuser, Basel, https://doi.org/10.1007/9783034879095, 2004. a, b, c
Da Prato, G.: An Introduction to InfiniteDimensional Analysis, Springer, Berlin, Heidelberg, https://doi.org/10.1007/3540290214, 2006. a, b, c
Da Prato, G. and Zabczyk, J.: Stochastic Equations in Infinite Dimensions, Cambridge University Press, https://doi.org/10.1017/cbo9781107295513, 2014. a, b, c
Del Sarto, G., Flandoli, F., Kuna, T., and Bröcker, J.: Variational Techniques for a OneDimensional Energy Balance Model (Version MatlabR2023b), Zenodo [code], https://doi.org/10.5281/zenodo.10469451, 2024. a
Díaz, J. I.: On the mathematical treatment of energy balance climate models, in: The Mathematics of Models for Climatology and Environment, edited by: Díaz, J. I., Springer, Berlin Heidelberg, Berlin, Heidelberg, 217–251, ISBN 9783642606038, 1997. a
Díaz, J., Langa, J., and Valero, J.: On the asymptotic behaviour of solutions of a stochastic energy balance climate model, Phys. D, 238, 880–887, https://doi.org/10.1016/j.physd.2009.02.010, 2009. a
Gaspar, P. and Guaraco, M. A.: The Allen–Cahn equation on closed manifolds, Calc. Var. Partial Dif., 57, 1–42, 2018. a
Ghil, M.: Climate Stability for a SellersType Model, J. Atmos. Sci., 33, 3–20, https://doi.org/10.1175/15200469(1976)033<0003:CSFAST>2.0.CO;2, 1976. a
Ghil, M. and Childress, S.: Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics, Springer, New York, https://doi.org/10.1007/9781461210528, 1987. a
Ghil, M. and Lucarini, V.: The physics of climate variability and climate change, Rev. Mod. Phys., 92, 035002, https://doi.org/10.1103/RevModPhys.92.035002, 2020. a, b
Hasselmann, K.: Stochastic climate models Part I. Theory, Tellus, 28, 473–485, https://doi.org/10.1111/j.21533490.1976.tb00696.x, 1976. a
Huang, Y., Xia, Y., and Tan, X.: On the pattern of CO_{2} radiative forcing and poleward energy transport, J. Geophys. Res., 122, 10578–10593, 2017. a
Imkeller, P.: Energy balance models – viewed from stochastic dynamics, in: Stochastic Climate Models, edited by: Imkeller, P. and von Storch, J.S., Birkhäuser Basel, Basel, 213–240, ISBN 9783034882873, 2001. a, b
IPCC: Climate Change 2001: The Scientific Basis, Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Houghton, J. T., Ding, Y., Griggs, D. J., Noguer, M., van der Linden, P. J., Dai, X., Maskell, K., and Johnson, C. A., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 881 pp., ISBN 9780521014953, 2001. a, b
Intergovernmental Panel on Climate Change (IPCC): Climate Change 2013 – The Physical Science Basis: Working Group I25 Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, https://doi.org/10.1017/CBO9781107415324, 2014. a
Jabri, Y.: The Mountain Pass Theorem, Cambridge University Press, https://doi.org/10.1017/cbo9780511546655, 2003. a, b
Lelièvre, T. and Stoltz, G.: Partial differential equations and stochastic methods in molecular dynamics, Acta Numer., 25, 681–880, https://doi.org/10.1017/s0962492916000039, 2016. a
Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping elements in the Earth's climate system, P. Natl. Acad. Sci. USA, 105, 1786–1793, 2008. a
Lenton, T. M., Livina, V. N., , V., van Nes, E. H., and Scheffer, M.: Early warning of climate tipping points from critical slowing down: comparing methods to improve robustness, Philos. T. Roy. Soc. A, 370, 1185–1204, https://doi.org/10.1098/rsta.2011.0304, 2012. a
Lucarini, V. and Bódai, T.: Transitions across melancholia states in a climate model: Reconciling the deterministic and stochastic points of view, Phys. Rev. Lett., 122, 158701, https://doi.org/10.1103/PhysRevLett.122.158701, 2019. a
Lucarini, V., Serdukova, L., and Margazoglou, G.: Lévy noise versus Gaussiannoiseinduced transitions in the Ghil–Sellers energy balance model, Nonlin. Processes Geophys., 29, 183–205, https://doi.org/10.5194/npg291832022, 2022. a
Myhre, G., Highwood, E. J., Shine, K. P., and Stordal, F.: New estimates of radiative forcing due to well mixed greenhouse gases, Geophys. Res. Lett., 25, 2715–2718, https://doi.org/10.1029/98gl01908, 1998. a
North, G. R.: Theory of EnergyBalance Climate Models, J. Atmos. Sci., 32, 2033–2043, https://doi.org/10.1175/15200469(1975)032<2033:TOEBCM>2.0.CO;2, 1975. a, b
North, G. R.: Multiple solutions in energy balance climate models, Global Planet. Change, 2, 225–235, https://doi.org/10.1016/09218181(90)90003U, 1990. a
North, G. R. and Cahalan, R. F.: Predictability in a Solvable Stochastic Climate Model., J. Atmos. Sci., 38, 504–513, https://doi.org/10.1175/15200469(1981)038<0504:PIASSC>2.0.CO;2, 1981. a
North, G. R. and Kim, K.Y.: Energy Balance Climate Models, Wiley, https://doi.org/10.1002/9783527698844, 2017. a, b, c, d, e
North, G. R., Howard, L., Pollard, D., and Wielicki, B.: Variational Formulation of BudykoSellers Climate Models, J. Atmos. Sci., 36, 255–259, https://doi.org/10.1175/15200469(1979)036<0255:VFOBSC>2.0.CO;2, 1979. a, b
North, G. R., Cahalan, R. F., and Coakley, J. A.: Energy balance climate models, Rev. Geophys., 19, 91–121, https://doi.org/10.1029/rg019i001p00091, 1981. a
Quarteroni, A. and Valli, A.: Numerical approximation of partial differential equations, vol. 23, Springer Science & Business Media, https://doi.org/10.1007/9783540852681, 2008. a
Rudin, W.: Principles of Mathematical Analysis, 3rd Edn., ISBN 9780070542358, 1976. a
Scheffer, M., Bascompte, J., Brock, W., Brovkin, V., Carpenter, S., Dakos, V., Held, H., Nes, E., Rietkerk, M., and Sugihara, G.: EarlyWarning Signals for Critical Transitions, Nature, 461, 53–9, https://doi.org/10.1038/nature08227, 2009. a
Sellers, W. D.: A Global Climatic Model Based on the Energy Balance of the EarthAtmosphere System, J. Appl. Meteorol. Clim., 8, 392–400, https://doi.org/10.1175/15200450(1969)008<0392:AGCMBO>2.0.CO;2, 1969. a, b
Smoller, J.: Shock waves and reaction–diffusion equations, vol. 258, Springer Science & Business Media, https://doi.org/10.1007/9781461208730, 2012. a
Temam, R.: InfiniteDimensional Dynamical Systems in Mechanics and Physics, Springer, New York, https://doi.org/10.1007/9781461206453, 1997. a
Thomas, J. W.: Numerical partial differential equations: finite difference methods, vol. 22, Springer Science & Business Media, https://doi.org/10.1007/9781489972781, 2013. a
Zhang, X., Li, X., Chen, D., Cui, H., and Ge, Q.: Overestimated climate warming and climate variability due to spatially homogeneous CO_{2} in climate modeling over the Northern Hemisphere since the mid19th century, Sci. Rep.UK, 9, 17426, https://doi.org/10.1038/s41598019535137, 2019. a