Variational techniques for a one-dimensional energy balance model

A one-dimensional 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 latitude-dependent albedo. Sufficient conditions are provided for the existence of at least three steady-state 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 time-dependent 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 non-decreasing 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.

1 Introduction

1.1 Low-dimensional 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 ice-albedo feedback and the Stefan–Boltzmann law. The positive ice-albedo 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 multi-stable 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₂), 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 noise-induced transitions resulting from unresolved processes in climate models or the representation of short-timescale 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 non-linear 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 zero-dimensional (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

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{C}_T{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}& ={\stackrel{\mathrm{‾}}{Q}}_{\mathrm{0}}\mathit{\beta }\left(T\right)+q-{\mathit{\sigma }}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{0}}{T}^{\mathrm{4}},t>\mathrm{0},\\ T_{\mid t=\mathrm{0}}& =T_{\mathrm{0}}.\end{array}\end{array}$$

In this equation, C_T>0 represents the heat capacity, ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{0}}>\mathrm{0}$ is the globally averaged solar radiation and the co-albedo β 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₂ on the energy budget (Bastiaansen et al., 2022). The term ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{e}}\left(T\right)={\mathit{\sigma }}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{0}}{T}^{\mathrm{4}}$ on the right-hand side of Eq. (1) accounts for the outgoing solar radiation, following the Stefan–Boltzmann law (where σ₀ denotes the Stefan–Boltzmann constant, and ε₀ is the globally averaged emissivity). The fixed points of the model are the solutions of the equation:

$${\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\mathrm{0},$$

corresponding to points in Fig. 1 where the absorbed radiation ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{a}}\left(T\right)={\stackrel{\mathrm{‾}}{Q}}_{\mathrm{0}}\mathit{\beta }\left(T\right)+q$ and the emitted radiation ${\stackrel{\mathrm{‾}}{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{‾}}{F}$ for the negative radiation budget $\stackrel{\mathrm{‾}}{R}$. In other words, $\stackrel{\mathrm{‾}}{F}$ is any regular function such that

$${\stackrel{\mathrm{‾}}{F}}^{\prime }\left(T\right)=\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\left(T\right)-\stackrel{\mathrm{‾}}{R_{\mathrm{a}}}\left(T\right)=-\stackrel{\mathrm{‾}}{R}\left(T\right).$$

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

$$\begin{array}{}\text{(2)}& \mathrm{d}T=\stackrel{\mathrm{‾}}{R}\left(T\right)\mathrm{d}t+\mathit{\epsilon }\mathrm{d}{W}_t,\end{array}$$

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{‾}}{\mathit{\nu }}$ (Lelièvre and Stoltz, 2016). An invariant measure is a probability distribution $\stackrel{\mathrm{‾}}{\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{‾}}{\mathit{\nu }}$ at some time t, then it remains so for all later times. It is given by

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{\mathit{\nu }}\left(\mathrm{d}T\right)={\displaystyle \frac{\mathrm{1}}{Z}}\exp \left(-{\displaystyle \frac{\mathrm{2}}{{\mathit{\epsilon }}^{\mathrm{2}}}}\stackrel{\mathrm{‾}}{F}\left(T\right)\right)\mathrm{d}T,\end{array}$$

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{‾}}{R}$ should be extended to negative values for the Kelvin temperature T in a way such that $\stackrel{\mathrm{‾}}{F}\to +\mathrm{\infty }$ as $T\to -\mathrm{\infty }$). The key observation from the explicit formula (Eq. 3) is that $\stackrel{\mathrm{‾}}{\mathit{\nu }}$ is concentrated around the minimum points of the function $\stackrel{\mathrm{‾}}{F}$. Indeed, if T₀ is a strict minimum point and T₁≠T₀ is a point close to T₀ such that $\stackrel{\mathrm{‾}}{F}\left(T_{\mathrm{1}}\right)>\stackrel{\mathrm{‾}}{F}\left(T_{\mathrm{0}}\right)$, then the mass given by the measure ν in a small neighbourhood of T₁ is exponentially lower than the mass around T₀; more specifically, the ratio between the two masses is given by $\exp \left(-\frac{\mathrm{2}}{{\mathit{\epsilon }}^{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{F}\left(T_{\mathrm{1}}\right)-\stackrel{\mathrm{‾}}{F}\left(T_{\mathrm{0}}\right)\right)\right)$.

Figure 1 (a) Absorbed radiation ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{a}}$ and emitted radiation ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{e}}$ for a 0D EBM. The graphs intersect in the three fixed points of the model $T_{\mathrm{S}}<T_{\mathrm{M}}<T_{\mathrm{W}}$; T_S and T_W are stable, and T_M is unstable. (b) Double-well potential $\stackrel{\mathrm{‾}}{F}$ associated with the 0D EBM. The function $\stackrel{\mathrm{‾}}{F}$ satisfies ${\stackrel{\mathrm{‾}}{F}}^{\prime }={\stackrel{\mathrm{‾}}{R}}_{\mathrm{e}}-{\stackrel{\mathrm{‾}}{R}}_{\mathrm{a}}$. The minimum points T_S and T_W of $\stackrel{\mathrm{‾}}{F}$ correspond to stable fixed points.

A one-dimensional (1D) EBM is given by a parabolic partial differential equation where the space variable is one-dimensional (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 non-linear 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 non-degenerate reaction–diffusion equation:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\partial }_{t}u& =\mathit{\kappa }\mathrm{\Delta }u+R(x,u;q),t>\mathrm{0},x\in (-\mathrm{1},\mathrm{1}),\\ u_x(t,-\mathrm{1})& =u_x(t,\mathrm{1})=\mathrm{0},t\ge \mathrm{0}\\ u(\mathrm{0},x)& =\stackrel{\mathrm{̃}}{u}\left(x\right),x\in [-\mathrm{1},\mathrm{1}],\end{array}\end{array}$$

where $\mathrm{\Delta }={\partial }_{xx}$ denotes the Laplace operator in dimension one, the Neumann boundary conditions impose no-heat flux at the poles and $\stackrel{\mathrm{̃}}{u}$ is an initial condition. The steady-state solutions of this model, representing the asymptotic solutions for the time-evolving dynamics, correspond to the non-negative solutions u=u(x) of the following elliptic problem:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{0}& =\mathit{\kappa }{u}^{\prime \prime }+R(x,u;q),x\in (-\mathrm{1},\mathrm{1}),\\ u^{\prime }(-\mathrm{1})& =u^{\prime }\left(\mathrm{1}\right)=\mathrm{0},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(6)}& F_q\left(u\right)=\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\mathcal{R}(x,u;q)\mathrm{d}x+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}\mathrm{|}\mathrm{|}{u}^{\prime }\mathrm{|}{\mathrm{|}}_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

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

$$\begin{array}{}\text{(7)}& \mathrm{d}u=\left(\mathit{\kappa }\mathrm{\Delta }u+R(x,u;q)\right)\mathrm{d}t+\mathit{\epsilon }\mathrm{d}{W}_t,\end{array}$$

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{‾}}{F}$ (see Sects. 2.1 and 3.1). Therefore, as in the zero-dimensional 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 steady-state 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” steady-state solutions for 1D EBMs. Firstly, in Theorem 1 we prove the existence of a steady-state solution for the 1D EBM by solving the associated variational problem

$$inf\left\{F_q\left(u\right)\mid u\in \mathbb{X}\right\},$$

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 steady-state 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 space-averaged global radiation balance ℛ of the 1D EBM should present a double-well potential with sufficiently deep minimum values attained at the two minimum points.

In essence, the conditions require that steady-state solutions of the spatially inhomogeneous 1D EBM (Eq. 5) are sufficiently well approximated by steady-state 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 space-averaged model. Then, the mountain pass theorem, a classical result from the calculus of variations, enables us to deduce the existence of a third steady-state 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

$$V\left(q\right)=inf\left\{F_q\left(u\right)\mid u\in \mathbb{X}\right\},$$

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 zero-measure 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 non-increasing.

• 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₀ for F_q. This establishes a one-to-one correspondence between the graph of V and the branch of the bifurcation diagram corresponding to u₀.

A byproduct of our analysis is that the global mean temperature is non-increasing 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 time-dependent 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 steady-state solutions. Secondly, we recall the properties of the steady-state 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 steady-state 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 non-increasing. 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 Background and methodology

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 Seller-type 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 non-degenerate parabolic PDE. Given an initial condition $\stackrel{\mathrm{̃}}{u}$, the non-linear, parabolic, reaction–diffusion PDE governing the model is given by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{C}_T{\displaystyle \frac{\partial u}{\partial t}}& ={\partial }_x\left[\mathit{\kappa }\left(x\right){\partial }_{x}u\right]+R_{\mathrm{a}}(x,u)-R_{\mathrm{e}}(u;q),\\ & t>\mathrm{0},x\in (-\mathrm{1},\mathrm{1}),\\ & \mathit{\kappa }(-\mathrm{1})u_x(t,-\mathrm{1})=\mathit{\kappa }\left(\mathrm{1}\right)u_x(t,\mathrm{1})=\mathrm{0},t\ge \mathrm{0}\\ & u(\mathrm{0},x)=\stackrel{\mathrm{̃}}{u}\left(x\right),x\in [-\mathrm{1},\mathrm{1}],\end{array}\end{array}$$

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:

$$R_{\mathrm{a}}(x,u)=Q_{\mathrm{0}}\left(x\right)(\mathrm{1}-\mathit{\alpha }(u\left)\right),$$

where Q₀ is the solar radiation per unit area, and α is the albedo. The solar radiation is assumed to be

$$Q_{\mathrm{0}}\left(x\right)={\widehat{Q}}_{\mathrm{0}}\left(c_{\mathrm{1}}-c_{\mathrm{2}}{x}^{\mathrm{2}}\right),c_i>\mathrm{0}$$

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,

$$\mathit{\alpha }\left(u\right)={\mathit{\alpha }}_{\mathrm{1}}+({\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{1}})\left[{\displaystyle \frac{\mathrm{1}+\tanh \left(K\right(u-u_{\mathrm{ref}}\left)\right)}{\mathrm{2}}}\right],$$

where K>0 is a rate parameter and α₁>α₂ 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

$$R_{\mathrm{e}}(u;q)={\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}{u}^{\mathrm{4}}-q,$$

where ε₀ and σ₀ are respectively the emissivity and Boltzmann's constant. The additive parameter q describes, in a simplified way, the radiative forcing by CO₂, i.e. the effect of atmospheric CO₂ 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₂ 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₂ concentration. If we linearize ${\widehat{R}}_{\mathrm{e}}$ with respect to temperature, we get

$${\widehat{R}}_{\mathrm{e}}(u,C)\approx A\left(C\right)+B(u-\widehat{u}),$$

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₂ is given by

$$A\left(C\right)=A_{\mathrm{1}}-A_{\mathrm{2}}\cdot \ln \left({\displaystyle \frac{C}{C_{\mathrm{0}}}}\right),$$

where C₀ is a reference CO₂ concentration, and $A_{\mathrm{1}},A_{\mathrm{2}}>\mathrm{0}$ are explicit constants. In conclusion,

$${\widehat{R}}_{\mathrm{e}}\approx A_{\mathrm{1}}+B(u-\widehat{u})-q,q=A_{\mathrm{2}}\cdot \ln \left({\displaystyle \frac{C}{C_{\mathrm{0}}}}\right),$$

and thus the radiative forcing of CO₂ has an additive structure. About point (ii), we adopt the well-mixed hypothesis for CO₂. In other words, we assume that atmospheric CO₂ is globally homogeneous, thereby inducing a radiative forcing q independent of latitude (IPCC, 2001). This assumption overlooks the spatial pattern of CO₂ concentration, which affects many aspects of the climate system, such as the poleward heat transport (Huang et al., 2017). It was the state-of-the-art 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 non-degenerate and given by

$$\mathit{\kappa }\left(x\right)=D(\mathrm{1}-x^{\mathrm{2}})+\mathit{\delta },D,\mathit{\delta }>\mathrm{0}.$$

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 non-negative, the solution remains non-negative 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{̃}}{u}\in H$, we consider the following SPDE:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\partial }_{t}u& =\mathit{\kappa }\mathrm{\Delta }u+Q_{\mathrm{0}}\left(x\right)\mathit{\beta }\left(u\right)-R_{\mathrm{e}}(u;q)+\mathit{\epsilon }\mathrm{d}{W}_t\\ u_{\mid t=\mathrm{0}}& =\stackrel{\mathrm{̃}}{u},\end{array}\end{array}$$

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 H-valued 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

$$\begin{array}{}\text{(10)}& \mathit{\nu }\left(\mathrm{d}u\right)={\displaystyle \frac{\mathrm{1}}{Z}}\exp \left(-{\displaystyle \frac{\mathrm{2}}{{\mathit{\epsilon }}^{\mathrm{2}}}}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\mathcal{R}(x,u;q)\mathrm{d}x\right)\mathit{\mu }\left(\mathrm{d}u\right),\end{array}$$

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 positive-definite, and its eigenvalues (λ_k)_k∈ℤ satisfy ${\sum }_{k\in \mathbb{Z}}{\mathit{\lambda }}_k<\mathrm{\infty }$. In the following lines, given a symmetric, positive-definite operator 𝒬 such that ${\sum }_{k\in \mathbb{Z}}{\mathit{\lambda }}_k<+\mathrm{\infty }$, we are going to explain how to construct an H-valued 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 H-valued random variable

$$X=\sum _{k\in \mathbb{Z}}\sqrt{{\mathit{\lambda }}_k}{R}_k{e}_k$$

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):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\left(\mathit{\kappa }\left(x\right)u^{\prime }\right)}^{\prime }& +Q_{\mathrm{0}}\left(x\right)\mathit{\beta }\left(u\right)+q-{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}{u}^{\mathrm{4}}=\mathrm{0},\\ & x\in (-\mathrm{1},\mathrm{1}),u^{\prime }(-\mathrm{1})=u^{\prime }\left(\mathrm{1}\right)=\mathrm{0},u\left(x\right)\ge \mathrm{0}.\end{array}\end{array}$$

These solutions can be either stable or unstable, depending on the long-term behaviour of their infinitesimal perturbations. As pointed out in Bastiaansen et al. (2022), if the reaction–diffusion equation for $u=u(t,x)$ was space-homogeneous, i.e. of the form

$$\begin{array}{}\text{(12)}& {\partial }_{t}u=\mathit{\kappa }\mathrm{\Delta }u+R\left(u\right),\end{array}$$

then the stable steady-state solutions would correspond to functions that are constant in space and time, with values given as the roots of

$$R\left(y\right)=\mathrm{0}.$$

A rigorous result in this direction has been shown in Gaspar and Guaraco (2018). Indeed, for a fixed double-well symmetric potential, it has been proved that (i) if κ is large enough, the only steady-state solutions of Eq. (12) are the constants where the potential is critical, and (ii) the number of unstable steady-state solutions to Eq. (12) can be made arbitrarily large as κ→0. Introducing a spatial dependence in $R=R(x,u)$ leads to a space-heterogeneous model. Depending on the space heterogeneity, it can exhibit any number of both stable and unstable steady-state solutions (Bastiaansen et al., 2022). The variational approach to the study of steady-state 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 steady-state solutions. That is, there exists q₁<q₂ s.t. Eq. (8) has one steady-state solution if q<q₁ or q>q₂, and the steady-state 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 steady-state solution for q<q₁ and by u_W the unique one for q>q₂. Figure 2a shows the bifurcation diagram of the model in the $(q,\stackrel{\mathrm{‾}}{u})$ plane, where $\stackrel{\mathrm{‾}}{u}={\int }_{-\mathrm{1}}^{\mathrm{1}}u\left(x\right)\mathrm{d}x$ denotes the average temperature. Figure 2b depicts the three steady-state solutions for $q=\mathrm{25}\in (q_{\mathrm{1}},q_{\mathrm{2}})$.

Figure 2 (a) Bifurcation diagram of the steady-state solutions in the $(q,\stackrel{\mathrm{‾}}{u})$ plane, with $\stackrel{\mathrm{‾}}{u}={\int }_{-\mathrm{1}}^{\mathrm{1}}u\left(x\right)\mathrm{d}x$. Solid lines denote stable solutions u_S and u_W, while dashed lines denote the unstable solution u_M. (b) Steady-state solutions of the EBM for q=25. In every point x of the space domain, the three steady-state solutions satisfy $u_{\mathrm{S}}\left(x\right)<u_{\mathrm{M}}\left(x\right)<u_{\mathrm{W}}\left(x\right)$, with maximum temperature attained at the Equator and minimum temperature attained at the poles.

A stability analysis can be conducted to determine the stability of the steady-state 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 saddle-node type, where the unstable solution u_M collides with either u_W (for q=q₁) or u_S (for q=q₂) and then disappears. These numerical findings regarding the number and stability of the steady-state solutions will be supported and validated using rigorous arguments, as in the next section.

Table 1 Parameters and constants appearing in the Seller EBM (8).

3 Results

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 space-averaged potential $\stackrel{\mathrm{‾}}{\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 steady-state 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

$$R_{\mathrm{e}}\left(u\right)=\left\{\begin{array}{ll}{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}{u}^{\mathrm{4}}& \text{ if }u\ge \mathrm{0},\\ \mathrm{0}& \text{ if }u<\mathrm{0}.\end{array}\right.$$

And β is smoothly extended to $\stackrel{\mathrm{̃}}{\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

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathit{\nu }\left(\mathrm{d}u\right)\propto \exp & \left(-{\displaystyle \frac{\mathrm{2}}{{\mathit{\epsilon }}^{\mathrm{2}}}}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}{\displaystyle \frac{{\left(u^{\mathrm{5}}\right)}_{+}}{\mathrm{5}}}-Q_{\mathrm{0}}\left(x\right)B\left(u\right)-qu\mathrm{d}x\right)\\ & \mathit{\mu }\left(\mathrm{d}u\right),\mathit{\mu }\sim \mathcal{N}\left(\mathrm{0},-{\displaystyle \frac{{\mathit{\epsilon }}^{\mathrm{2}}}{\mathrm{2}\mathit{\kappa }}}{\mathrm{\Delta }}^{-\mathrm{1}}\right),\end{array}\end{array}$$

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

$$\mathit{\mu }\left(\mathrm{d}u\right)={\displaystyle \frac{\mathrm{1}}{Z_{\mathrm{1}}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle {\mathcal{Q}}^{-\mathrm{1}}u,u\rangle \right)\mathrm{d}u,$$

where ${\mathcal{Q}}^{-\mathrm{1}}=-\frac{\mathrm{2}\mathit{\kappa }}{{\mathit{\epsilon }}^{\mathrm{2}}}\mathrm{\Delta }.$ Here, Z₁ 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

$$\mathit{\mu }\left(\mathrm{d}u\right)={\displaystyle \frac{\mathrm{1}}{Z_{\mathrm{1}}}}\exp \left({\displaystyle \frac{\mathit{\kappa }}{{\mathit{\epsilon }}^{\mathrm{2}}}}\langle u^{\prime \prime },u\rangle \right)\mathrm{d}u={\displaystyle \frac{\mathrm{1}}{Z_{\mathrm{1}}}}\exp \left(-{\displaystyle \frac{\mathit{\kappa }}{{\mathit{\epsilon }}^{\mathrm{2}}}}\mathrm{|}\mathrm{|}{u}^{\prime }\mathrm{|}{\mathrm{|}}_{\mathrm{2}}^{\mathrm{2}}\right)\mathrm{d}u.$$

Plugging the previous identity into Eq. (10), we obtain

$$\begin{array}{rl}\mathit{\nu }\left(\mathrm{d}u\right)& \propto \exp \left(-{\displaystyle \frac{\mathrm{2}}{{\mathit{\epsilon }}^{\mathrm{2}}}}\left(\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}{\displaystyle \frac{{\left(u^{\mathrm{5}}\right)}_{+}}{\mathrm{5}}}-Q_{\mathrm{0}}\left(x\right)B\left(u\right)\right.\right.\\ & \left.\left.-qu\mathrm{d}x+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}\mathrm{|}\mathrm{|}{u}^{\prime }\mathrm{|}{\mathrm{|}}_{\mathrm{2}}^{\mathrm{2}}\right)\right)\mathrm{d}u\\ & \propto \exp \left(-{\displaystyle \frac{\mathrm{2}}{{\mathit{\epsilon }}^{\mathrm{2}}}}{F}_q\left(u\right)\right)\mathrm{d}u.\end{array}$$

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₁ is a global minimum point and u≠u₁, 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₁; in particular, the ratio between the two masses is given by $\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 infinite-dimensional 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

$$\begin{array}{rl}{F}_q\left(u\right)& =\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{u^{\mathrm{5}}}{\mathrm{5}}}{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}-Q_{\mathrm{0}}\left(x\right)B\left(u\right)-qu\mathrm{d}x\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\mathit{\kappa }\left(x\right)\left(u^{\prime }\right(x){)}^{\mathrm{2}}\mathrm{d}x,\end{array}$$

where B is a primitive of the co-albedo $\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 steady-state 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

$$\begin{array}{rl}\mathit{\delta }& F_q(u,h)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}s}}{F}_q(u+sh{)}_{\mid s=\mathrm{0}}=\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\left(u^{\mathrm{4}}{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}-Q_{\mathrm{0}}\left(x\right)\mathit{\beta }\left(u\right)-q\right)\\ & h\mathrm{d}x+\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\mathit{\kappa }\left(x\right)u^{\prime }\left(x\right)h^{\prime }\left(x\right)\mathrm{d}x\\ & =\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\left[u^{\mathrm{4}}{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}-Q_{\mathrm{0}}\left(x\right)\mathit{\beta }\left(u\right)-q-\left(\mathit{\kappa }\right(x\left)u^{\prime }\right(x){)}^{\prime }\right]h\mathrm{d}x\end{array}$$

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 steady-state solution for the EBM. In particular, local extremum points for F_q correspond to steady-state 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.

Theorem 1.

If q>0, then there exists a global regular non-negative minimizer for F_q. In other words, if we consider the variational problem

$$\begin{array}{}\text{(14)}& inf\left\{F_q\left(u\right)\mid u\in H^{\mathrm{1}},u\ge \mathrm{0}\right\},\end{array}$$

then there exists u₀∈C^∞ such that u₀ is a solution of the EBM and

$$F_q\left(u_{\mathrm{0}}\right)=inf\left\{F_q\left(u\right)\mid u\in H^{\mathrm{1}},u\ge \mathrm{0}\right\}.$$

In addition to this, if q belongs to a bounded interval, then u can be bounded uniformly with respect to q:

$$\begin{array}{}\text{(15)}& \exists M>\mathrm{0}\mathit{\text{ s.t. }}{u}_{\mathrm{0}}\left(x\right)\le M,\forall x\in [-\mathrm{1},\mathrm{1}].\end{array}$$

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{‾}}{\mathcal{R}}:\mathbb{R}\to \mathbb{R}$ coming from the space-averaged model

$$\stackrel{\mathrm{‾}}{\mathcal{R}}\left(u\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\mathcal{R}(x,u)\mathrm{d}x.$$

If the viscosity κ>0 is sufficiently large and the function $\stackrel{\mathrm{‾}}{\mathcal{R}}$ has a double-well 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 steady-state solution. Next, we characterize a situation in which there are three steady-state solutions, two of which are local minimizers (Jabri, 2003). This is summarized in the following result.

Theorem 2.

Denote by $B_{H^{\mathrm{1}}}\mathrm{(}v\mathrm{,}\mathit{\rho }\mathrm{)}\mathrm{=}\mathit{\{}u\mathrm{\in }{H}^{\mathrm{1}}\mathrm{\mid }\mathrm{|}\mathrm{|}u\mathrm{-}v\mathrm{|}{\mathrm{|}}_{H^{\mathrm{1}}}\mathrm{<}\mathit{\rho }\mathit{\}}$ the open ball in H¹ with centre v and radius ρ>0. Assume $\stackrel{\mathrm{‾}}{\mathrm{R}}$ has two non-negative minimum points u₁≠u₂, with F_q(u₁)≥F_q(u₂). 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{‾}}{\mathit{\epsilon }}\mathrm{>}\mathrm{0}$ satisfies

• i.

  ${\stackrel{\mathrm{‾}}{\mathcal{R}}}^{\prime \prime }\left(u_i\right)>f\left(\stackrel{\mathrm{‾}}{\mathit{\epsilon }}\right)$, for $i=\mathrm{1},\mathrm{2}$;

• ii.

  $\mathit{\kappa }>g\left(\stackrel{\mathrm{‾}}{\mathit{\epsilon }}\right)$;

• iii.

  $\stackrel{\mathrm{‾}}{\mathit{\epsilon }}\le \mathit{\omega }$;

then ${\stackrel{\mathrm{̃}}{F}}_q$ has two local minimum points, ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{u}}_{\mathrm{2}}$, such that

• a.

  $B_{H^{\mathrm{1}}}(u_{\mathrm{1}},\stackrel{\mathrm{‾}}{\mathit{\epsilon }})\cap B_{H^{\mathrm{1}}}(u_{\mathrm{2}},\stackrel{\mathrm{‾}}{\mathit{\epsilon }})=\mathrm{\varnothing }$;

• b.

  ${\stackrel{\mathrm{̃}}{u}}_i\in B_{H^{\mathrm{1}}}(u_i,\stackrel{\mathrm{‾}}{\mathit{\epsilon }}),$ for $i=\mathrm{1},\mathrm{2}$.

• c.

  If $\mathrm{|}\mathrm{|}u-u_{\mathrm{1}}\mathrm{|}{\mathrm{|}}_{H^{\mathrm{1}}}=\stackrel{\mathrm{‾}}{\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{‾}}{\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 space-inhomogeneous EBM to the stable solution of the corresponding space-averaged 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

$$V\left(q\right)=inf\left\{F_q\left(u\right)\mid u\in H^{\mathrm{1}},u\ge \mathrm{0}\right\}.$$

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 steady-state 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₃ such that u_S is the global minimum point of F_q for q<q₃, while u_W is the global minimum point for q>q₃. Further, for q=q₃ the function F_q has two different global minimum points u_S and u_W, and q=q₃ corresponds a non-differentiability point for V. In addition to this, the value function appears to be concave, thus with a decreasing derivative, where it exists.

Figure 3 Potential functional F_q evaluated in the three steady-state solutions u_S,u_M and u_W. For q<q₃, u_S is the global minimum point, while u_W is a local minimum point. On the other hand, for q>q₃ the opposite happens. Solid lines correspond to values of the functional attained on stable solutions; dashed lines are for values corresponding to unstable ones.

Summarizing, the numerical evaluations of V(q) suggest the following result, a rigorous proof of which is included in the Supplement.

Theorem 3.

Assume q belongs to a bounded interval. Then

• i.

  V is Lipschitz continuous.

• ii.

  q is a non-differentiable point for V if and only if there is more than one minimizer for F_q.

• iii.

  V is concave and V^′ is non-increasing.

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₀ 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}{g}_i$ is Lipschitz continuous if the constants L_i can be bounded uniformly. In our case, given u∈H¹ is non-negative, we have

$$\mathrm{|}{F}_{{\mathit{\mu }}_{\mathrm{1}}}\left(u\right)-F_{{\mathit{\mu }}_{\mathrm{2}}}\left(u\right)\mathrm{|}\le \mathrm{|}{\mathit{\mu }}_{\mathrm{1}}-{\mathit{\mu }}_{\mathrm{2}}\mathrm{|}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}\mathrm{|}u\left(x\right)\mathrm{|}\mathrm{d}x\le \mathrm{2}M\mathrm{|}{\mathit{\mu }}_{\mathrm{1}}-{\mathit{\mu }}_{\mathrm{2}}\mathrm{|}$$

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

Figure 4 Comparison between the value function graph (left) and bifurcation diagram (right) for the 1D EBM. The magenta-shaded area highlights the parts of the plots which are in one-to-one correspondence. (a) Functional F_q evaluated on steady-state solutions, as in Fig. 3. (b) Bifurcation diagram, as in Fig. 2a.

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 non-decreasing with respect to q. In other words, as the concentration of CO₂ 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{̃}}{u}$, this model is given by

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\partial }_{t}u& ={\displaystyle \frac{\mathrm{1}}{\mathrm{100}}}\mathrm{\Delta }u+u(\mathrm{1}-u^{\mathrm{2}})+q+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cos \left(\mathit{\pi }x\right),\\ & x\in (-\mathrm{1},\mathrm{1}),t>\mathrm{0},\\ u_x(t,-\mathrm{1})& =u_x(t,\mathrm{1})=\mathrm{0},t\ge \mathrm{0},\\ u_{\mid t=\mathrm{0}}& =\stackrel{\mathrm{̃}}{u}.\end{array}\end{array}$$

The associated elliptic problem for u=u(x) is

$$\begin{array}{}\text{(17)}& \begin{array}{rl}\mathrm{0}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{100}}}{u}^{\prime \prime }+u(\mathrm{1}-u^{\mathrm{2}})+q+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cos \left(\mathit{\pi }x\right),\\ & x\in (-\mathrm{1},\mathrm{1}),\\ u^{\prime }(-\mathrm{1})& =u^{\prime }\left(\mathrm{1}\right)=\mathrm{0}.\end{array}\end{array}$$

In this case, the potential functional takes the form

$$J_q\left(u\right)=\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{u^{\mathrm{4}}\left(x\right)}{\mathrm{4}}}-{\displaystyle \frac{u^{\mathrm{2}}\left(x\right)}{\mathrm{2}}}-u\left(x\right)(q+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cos (\mathit{\pi }x\left)\right)\mathrm{d}x,$$

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₅ or |q|<q₄, there exists a single steady-state solution, which is stable, (b) for $q_{\mathrm{4}}<\mathrm{|}q\mathrm{|}<q_{\mathrm{5}}$ there are three steady-state solutions, two of which are stable while the third is unstable. Further, $q=q_{\mathrm{4}},q_{\mathrm{5}}$ are bifurcation points of the saddle-node type. We denote by u₁ the steady-state solution for $q<-q_{\mathrm{5}}$, by u₂ and u₃ the steady-state solutions appearing at the bifurcation point $q=-q_{\mathrm{5}}$ and existing for $-q_{\mathrm{5}}<q<-q_{\mathrm{4}}$ in addition to u₁, and by u₄ and u₅ the steady-state solutions appearing at q=q₄ and existing for $q_{\mathrm{4}}<q<q_{\mathrm{5}}$ in addition to u₃. Regarding the potential functional J_q, in this case there exists $q_{\mathrm{6}}\in (q_{\mathrm{4}},q_{\mathrm{5}})$ such that u₁ is the global minimum point for the functional for $q<-q_{\mathrm{6}}$, and u₃ is the global minimum point for $-q_{\mathrm{6}}<q<q_{\mathrm{6}}$, while u₅ becomes the global minimum point for q>q₆. 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.

Figure 5 Comparison between the value function and the bifurcation diagram for the non-homogeneous ACE. The magenta-shaded area highlights the parts of the plots which are in one-to-one correspondence. (a) Potential functional evaluated on the steady-state solutions: u₁ is the global minimum point for $q<-q_{\mathrm{6}}$, u₃ is the global minimum point for $-q_{\mathrm{6}}<q<q_{\mathrm{6}}$, u₅ is the global minimum point for q>q₆. Note that $q=\pm q_{\mathrm{6}}$ are the non-differentiability point for the value function, corresponding to non-uniqueness of the minimizer. (b) Bifurcation diagram.

4 Conclusions

In this paper, we have considered a one-dimensional energy balance model depending on a bifurcation parameter q, describing the effect of CO₂ 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 steady-state solutions. The functional F_q has significant implications, as it is closely linked to both the stability of steady-state 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 steady-state 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 steady-state solutions u_S,u_M and u_W, we observed that the function exhibits Lipschitz continuity and concavity. Furthermore, non-differentiability points of V(q) coincide with points where multiple global minimizers exist for F_q. Lastly, when V is differentiable, its derivative is non-increasing 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)\mathrm{d}x$, where u₀ is the minimizer for F_q. Moreover, as a consequence of the explicit expression for V^′, the global mean temperature is non-decreasing with respect to q, and it is continuous, except for a Lebesgue zero-measure set of upwards jumps. These are the non-differentiability 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 space-inhomogeneous reaction–diffusion equations.

The diffusion function κ=κ(x) that we have examined is non-degenerate 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₂ 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₂. Finally, we have outlined simple conditions, adaptable to other multi-stable 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.

Appendix A: Numerical methods

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

$$f(x,u)=R_{\mathrm{e}}(x,u)-R_{\mathrm{a}}\left(u\right)$$

as the non-linear reaction term. We consider a uniform mesh for $[-\mathrm{1},\mathrm{1}]$ made of n+1 points

$$\begin{array}{rl}{x}_{\mathrm{0}}& =-\mathrm{1}<x_{\mathrm{1}}<\mathrm{\cdots }<x_n=\mathrm{1},x_i=-\mathrm{1}+i\mathrm{\Delta }x,\\ & i=\mathrm{0},\mathrm{\dots },n,\mathrm{\Delta }x={\displaystyle \frac{\mathrm{2}}{n}}.\end{array}$$

Then, the solution to the problem can be approximated by considering the following system:

$$\begin{array}{rl}{\displaystyle \frac{u_{i-\mathrm{1}}{\mathit{\kappa }}_{i-\frac{\mathrm{1}}{\mathrm{2}}}-u_i({\mathit{\kappa }}_{i-\frac{\mathrm{1}}{\mathrm{2}}}+{\mathit{\kappa }}_{i+\frac{\mathrm{1}}{\mathrm{2}}})+u_{i+\mathrm{1}}{\mathit{\kappa }}_{i+\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}}& +f(x_i,u_i)=\mathrm{0},\\ & i=\mathrm{0},\mathrm{\dots },n,\\ {\displaystyle \frac{u_{\mathrm{1}}-u_{-\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }x}}={\displaystyle \frac{u_{n+\mathrm{1}}-u_{n-\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }x}}& =\mathrm{0},\end{array}$$

where u₋₁and u_n+1 are ghost points, and

$$u_i=u\left(x_i\right),{\mathit{\kappa }}_{i\pm \frac{\mathrm{1}}{\mathrm{2}}}=\mathit{\kappa }\left(x_{i\pm \frac{\mathrm{1}}{\mathrm{2}}}\right),x_{i\pm \frac{\mathrm{1}}{\mathrm{2}}}:=x_i\pm \mathrm{\Delta }x/\mathrm{2}.$$

The system of equations can be written in vector form as

$$\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{x}^{\mathrm{2}}}}& \left[\begin{array}{cccccc}-\mathrm{\Delta }x/\mathrm{2}& \mathrm{0}& \mathrm{\Delta }x/\mathrm{2}& & & \\ {\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}& -({\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}})& {\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}}& & & \\ & \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }& & & \\ & & & {\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}& -({\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}})& {\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}}\\ & & & -\mathrm{\Delta }x/\mathrm{2}& \mathrm{0}& \mathrm{\Delta }x/\mathrm{2}\end{array}\right]\\ & \mathit{u}+\left[\begin{array}{c}\mathrm{0}\\ f_{\mathrm{0}}\\ \mathrm{⋮}\\ f_n\\ \end{array}\right]=\mathbf{0},\end{array}$$

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

$$-{\displaystyle \frac{{\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}{u}_{-\mathrm{1}}+{\displaystyle \frac{{\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}{u}_{\mathrm{0}}-{\displaystyle \frac{f_{\mathrm{0}}}{\mathrm{2}}}=\mathrm{0}.$$

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

$${\displaystyle \frac{{\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}{u}_n-{\displaystyle \frac{{\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}{u}_{n+\mathrm{1}}-{\displaystyle \frac{f_n}{\mathrm{2}}}=\mathrm{0}.$$

In this way, the Neumann version of the elliptic problem has the following form:

$$\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{x}^{\mathrm{2}}}}& \left[\begin{array}{cccccc}-{\displaystyle \frac{{\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}& \frac{{\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}}}{\mathrm{2}}& & & & \\ {\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}& -({\mathit{\kappa }}_{-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}})& {\mathit{\kappa }}_{\mathrm{1}/\mathrm{2}}& & & \\ & \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }& & & \\ & & & {\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}& -({\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}})& {\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}}\\ & & & & \frac{{\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}}}{\mathrm{2}}& -\frac{{\mathit{\kappa }}_{n-\mathrm{1}/\mathrm{2}}+{\mathit{\kappa }}_{n+\mathrm{1}/\mathrm{2}}}{\mathrm{2}}\end{array}\right]\\ & \mathit{u}+\left[\begin{array}{c}-f_{\mathrm{0}}/\mathrm{2}\\ f_{\mathrm{0}}\\ \mathrm{⋮}\\ f_n\\ -f_n/\mathrm{2}\end{array}\right]=\mathbf{0}\end{array}$$

and consists of a set of (n+3) non-linear 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.