Solving a North-type energy balance model using boundary integral methods

Simplified climate models, such as energy balance models (EBMs), are useful conceptual tools, in part because their reduced complexity often allows for studies using analytical methods. In this paper, we solve a North-type EBM using a boundary integral method (BIM). The North-type EBM is a diffusive, one-dimensional EBM with a nonlinear albedo feedback mechanism. We discuss this approach in light of existing analytical techniques for this type of equation. Subsequently, we test the proposed method by solving multiple North-type EBMs with a zonally symmetric continent featuring an altered ice-albedo feedback dynamic. We demonstrate that the introduction of a continent results in new equilibrium states characterized by multiple ice edges and ice belts. Furthermore, we show that the BIM serves as an efficient framework for handling unconventional ice distributions and model configurations for North-type EBMs.

1 Introduction

Despite the advancement in computational power, conceptual climate models remain valuable tools for understanding the Earth's climate system. The complexity of realistic models has highlighted the need for a hierarchical model structure where conceptual models provide a solid theoretical foundation as model complexity increases (Schneider and Dickinson, 1974; Claussen et al., 2002; McGuffie and Henderson-Sellers, 2014). Energy balance models (EBMs) stand out as some of the simplest climate models. Their simplicity allows for both analytical and numerical studies of climate responses to forcings. So-called zero-dimensional EBMs describe the Earth's global mean temperature and include no spatial variables. Zero-dimensional EBMs are readily examined using analytical tools (North, 1990; Ghil and Lucarini, 2020; Lohmann, 2020). One-dimensional EBMs with a latitude dependence, coupled with a linear transport term and a nonlinear albedo feedback, often called Budyko-type models, also lend themselves to analytical investigations (Budyko, 1969; Held and Suarez, 1974; Widiasih, 2013; Walsh and Widiasih, 2014). A more physically motivated transport mechanism (Rose and Marshall, 2009) may be included in the model by instead adding a diffusion term. One-dimensional EBMs with meridional heat transport by diffusion, hereafter called North-type models, have attracted considerable interest (North et al., 1981; Ghil, 1976; Bódai et al., 2015; Del Sarto et al., 2024). While the inclusion of the diffusion term complicates the models, mathematically speaking, in certain configurations these models may be studied using analytical methods. Pioneering analytical investigations into these models have been conducted by North (North, 1975a, b; North et al., 1981). A general solution expressed through a Fourier–Legendre series for the equilibrium temperature field was found using spectral methods for a step function albedo. Although this solution is rapidly converging for standard model configurations such as the idealized aquaplanet, a geographical input with parameter discontinuities across land–sea boundaries causes the spectral solution to converge slowly (Mengel et al., 1988; North and Kim, 2017).

In this paper, we solve the stationary form of an EBM with a meridional heat transport and a nonlinear albedo feedback using an analytical method: the boundary integral method (BIM). This North-type EBM describes the zonal mean surface temperature with its key features being a linear heat diffusion across latitudes, an ice-albedo feedback mechanism and the stabilizing effect of the outgoing longwave radiation. The model may be formulated as

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle }C{\displaystyle \frac{\partial T}{\partial t}}-D{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}& {\displaystyle }\left(\sin \mathit{\theta }{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\right)+BT=\\ {\displaystyle }& {\displaystyle }Qs\left(\mathit{\theta }\right)(\mathrm{1}-a(T\left)\right)-A\end{array}\end{array}$$

using spherical coordinates, where the polar angle, θ, is the latitude. The latitude ranges from θ=0, the north pole, to θ=π, the south pole. Here C represents the heat capacity of the lower atmosphere and hydrosphere with an assigned average constant value, and D is the diffusion rate determining the strength of the meridional heat transport. The parameter Q is defined as one-fourth of the mean total solar irradiance (TSI), as the disk silhouette capturing solar radiation is one-fourth of the Earth's total area. To address the solar radiation distribution across latitudes, the model incorporates an average annual latitudinal energy distribution function, denoted s(θ). Additionally, the model assumes a constant lapse rate, establishing a linear relationship between surface temperature and outgoing energy, expressed as $E_{\text{out}}=A+BT$, where A and B are constants (Budyko, 1969). The ice-albedo feedback mechanism is included by allowing for a temperature-dependent albedo,

$$\begin{array}{}\text{(2)}& a\left(T\right)=\left\{\begin{array}{ll}{a}_{\mathrm{1}},& T>-T_{\mathrm{s}},\\ a_{\mathrm{2}},& T<-T_{\mathrm{s}},\end{array}\right.\end{array}$$

where −T_s is the critical temperature for ice formation at the surface. Latitudes with an annual mean temperature below −T_s are deemed to have an ice cover. Consequently, there will be a critical latitude at which the ice cover ends and begins. A major challenge arises in determining the location of the critical latitude under the given constraints. We show that the BIM offers a convenient way to address this, even for scenarios with several critical latitudes.

The proposed method is tested on a model configuration where the idealized aquaplanet is given a zonally symmetric continent with an altered ice-albedo feedback mechanism. Equilibrium solutions to Eq. (1) are found, and a bifurcation diagram is drawn for three different systems with a zonally symmetric continent. Lin and North (1990) previously studied similar zonal band continent configurations, and a circular cap of land centered at one pole was studied by Mengel et al. (1988). In these studies, land and sea were differentiated by a change in heat capacity, with no change in the stationary equation. Here, we demonstrate that the BIM represents an analytical method that efficiently handles arbitrary parameter discontinuities at the land–sea interface in North-type EBMs. Additionally, we show that the introduction of a continent with altered equilibrium parameters gives rise to new equilibrium states, characterized by unconventional ice distributions featuring multiple ice edges and ice belts.

2 Results

BIMs are a general approach for boundary value problems, where the problem is reduced to boundary integral equations involving an associated Green's function (Hsiao and Wendland, 2008; Morino and Piva, 2012). To showcase the application of the BIM in the context of EBMs, we employ it to find the equilibrium solutions to the classical idealized aquaplanet. Using T_s as a scale for temperature, $T=T_{\mathrm{s}}\stackrel{\mathrm{̃}}{T}$, where $\stackrel{\mathrm{̃}}{T}=\stackrel{\mathrm{̃}}{T}\left(\mathit{\theta }\right)$ is the non-dimensional temperature field at latitude θ, we may write the stationary form of the energy balance equation (Eq. 1) in non-dimensional form:

$$\begin{array}{}\text{(3)}& -{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{T}}{\partial \mathit{\theta }}}\right)+\mathit{\beta }\stackrel{\mathrm{̃}}{T}=\mathit{\eta }s\left(\mathit{\theta }\right)(\mathrm{1}-a(\stackrel{\mathrm{̃}}{T}\left)\right)-\mathit{\alpha },\end{array}$$

where $\mathit{\beta }=\frac{B}{D}$, $\mathit{\alpha }=\frac{A}{T_{\mathrm{s}}D}$ and $\mathit{\eta }=\frac{Q}{T_{\mathrm{s}}D}$. Hereafter, the tilde notation on the non-dimensional temperatures will be omitted, as the subsequent analysis will focus exclusively on non-dimensional temperatures. Defining

$$\begin{array}{}\text{(4)}& \mathcal{L}(\cdot )=-{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}(\cdot )\right)+\mathit{\beta }(\cdot ),\end{array}$$

it can be shown that, for two functions v and u on the domain [θ₁,θ₂], we have

$$\begin{array}{}\text{(5)}& \underset{{\mathit{\theta }}_{\mathrm{1}}}{\overset{{\mathit{\theta }}_{\mathrm{2}}}{\int }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }\mathit{\{}v\mathcal{L}u-u\mathcal{L}v\mathit{\}}=\left\{u\sin \mathit{\theta }{\displaystyle \frac{\partial v}{\partial \mathit{\theta }}}-v\sin \mathit{\theta }{\displaystyle \frac{\partial u}{\partial \mathit{\theta }}}\right\}{|}_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}.\end{array}$$

Defining

$$\begin{array}{}\text{(6)}& h(T,\mathit{\theta })=\mathit{\eta }s\left(\mathit{\theta }\right)(\mathrm{1}-a(T\left)\right)-\mathit{\alpha },\end{array}$$

we may write Eq. (3) in the compact form:

$$\begin{array}{}\text{(7)}& \mathcal{L}T=h.\end{array}$$

Let K be a Green's function for operator (4). That is,

$$\begin{array}{}\text{(8)}& \mathcal{L}K(\mathit{\theta },\mathit{\xi })={\mathit{\delta }}_{\mathit{\xi }}\left(\mathit{\theta }\right),\end{array}$$

where δ_ξ(θ) is a Dirac delta function along the curve $\mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$ (see Appendix A, Eq. A4). Inserting v=T and u=K into the identity in Eq. (5), we have

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\int }_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }& {\displaystyle }\mathit{\{}T\mathcal{L}K-K\mathcal{L}T\mathit{\}}=\\ {\displaystyle }& {\displaystyle }\left\{K\sin \mathit{\theta }{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}-T\sin \mathit{\theta }{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}\right\}{|}_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}.\end{array}\end{array}$$

Inserting Eqs. (7) and (8) into Eq. (9) and rearranging for T, any exact solution T(θ) to Eq. (3) must satisfy the identity

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathit{\xi }\right)& {\displaystyle }={\int }_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h(T,\mathit{\theta })\\ {\displaystyle }& {\displaystyle }+\left\{K\sin \mathit{\theta }{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}-T\sin \mathit{\theta }{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}\right\}{|}_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}.\end{array}\end{array}$$

A suitable Green's function is found in Appendix A.

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }K(\mathit{\theta },\mathit{\xi })=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\frac{P_{\mathit{\lambda }}\left(\cos \mathit{\xi }\right)\left(\mathit{\pi }\cot \right(\mathit{\pi }\mathit{\lambda }\left)P_{\mathit{\lambda }}\right(\cos \mathit{\theta })-\mathrm{2}{Q}_{\mathit{\lambda }}(\cos \mathit{\theta }\left)\right)}{\mathrm{2}(\mathrm{1}+\mathit{\lambda })\left(P_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }+\mathrm{1}}\right(\cos \mathit{\xi })-P_{\mathit{\lambda }+\mathrm{1}}(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)\right)},& \mathit{\theta }>\mathit{\xi },\\ \frac{P_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)\left(\mathit{\pi }\cot \right(\mathit{\pi }\mathit{\lambda }\left)P_{\mathit{\lambda }}\right(\cos \mathit{\xi })-\mathrm{2}{Q}_{\mathit{\lambda }}(\cos \mathit{\xi }\left)\right)}{\mathrm{2}(\mathrm{1}+\mathit{\lambda })\left(P_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }+\mathrm{1}}\right(\cos \mathit{\xi })-P_{\mathit{\lambda }+\mathrm{1}}(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)\right)},& \mathit{\theta }<\mathit{\xi }.\end{array}\right.\end{array}\end{array}$$

Here P_λ and Q_λ are Legendre functions of order λ, where

$$\begin{array}{}\text{(12)}& \mathit{\lambda }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\sqrt{\mathrm{1}-\mathrm{4}\mathit{\beta }}-\mathrm{1}\right).\end{array}$$

This Green's function is continuous and bounded on the domain $\mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$, and its derivative is also bounded at the boundaries θ=0 and θ=π for a given $\mathit{\xi }\in [\mathrm{0},\mathit{\pi }]$.

2.1 No partial ice cover

The BIM relies on initially positing an ansatz regarding the distribution of ice and water, and then the domain is partitioned into regions where the albedo function remains invariant with respect to temperature. Subsequently, the identity in Eq. (10) is applied within these regions to obtain explicit expressions for the solution to Eq. (3). We start by analyzing solutions where the surface has no partial ice cover. This is the linear problem where the ice-albedo feedback is inactive due to extreme temperatures. The step function albedo (Eq. 2) leads to

$$\begin{array}{}\text{(13)}& h(T,\mathit{\theta })=\left\{\begin{array}{ll}\mathit{\eta }s\left(\mathit{\theta }\right)(\mathrm{1}-a_{\mathrm{1}})-\mathit{\alpha },& T>-\mathrm{1},\\ \mathit{\eta }s\left(\mathit{\theta }\right)(\mathrm{1}-a_{\mathrm{2}})-\mathit{\alpha },& T<-\mathrm{1}.\end{array}\right.\end{array}$$

Note that h is a function of non-dimensional temperature; hence, the critical temperature for the presence of surface ice is $T=-\mathrm{1}$. For theses extreme cases, where $T>-\mathrm{1}\forall \mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$ and $T<-\mathrm{1}\forall \mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$, the surface is either (1) devoid of ice entirely or (2) entirely covered by ice. Consequently, the function h is constant in the domain [0,π], and we may apply the relation in Eq. (10) to the full domain: letting ${\mathit{\theta }}_{\mathrm{1}}\to {\mathrm{0}}^{+}$ and ${\mathit{\theta }}_{\mathrm{2}}\to {\mathit{\pi }}^{-}$, we get

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathit{\xi }\right)& {\displaystyle }={\int }_{\mathrm{0}}^{\mathit{\pi }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h(T,\mathit{\theta })\\ {\displaystyle }& {\displaystyle }+\underset{{\mathit{\theta }}_{\mathrm{2}}\to {\mathit{\pi }}^{-}}{lim}K({\mathit{\theta }}_{\mathrm{2}},\mathit{\xi })\sin {\mathit{\theta }}_{\mathrm{2}}{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }-\underset{{\mathit{\theta }}_{\mathrm{2}}\to {\mathit{\pi }}^{-}}{lim}T\left({\mathit{\theta }}_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{2}}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}({\mathit{\theta }}_{\mathrm{2}},\mathit{\xi })\\ {\displaystyle }& {\displaystyle }-\underset{{\mathit{\theta }}_{\mathrm{1}}\to {\mathrm{0}}^{+}}{lim}K({\mathit{\theta }}_{\mathrm{1}},\mathit{\xi })\sin {\mathit{\theta }}_{\mathrm{1}}{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+\underset{{\mathit{\theta }}_{\mathrm{1}}\to {\mathrm{0}}^{+}}{lim}T\left({\mathit{\theta }}_{\mathrm{1}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}({\mathit{\theta }}_{\mathrm{1}},\mathit{\xi }).\end{array}\end{array}$$

Green's function, K, and its derivative are bounded at the boundary (see Appendix A). The gradient must vanish at the boundaries (North, 1975a) as we do not allow for heat transport at the poles, leading to the boundary conditions

$$\begin{array}{}\text{(15)}& \underset{\mathit{\theta }\to \mathrm{0}}{lim}\sin \mathit{\theta }{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left(\mathit{\theta }\right)=\mathrm{0}\end{array}$$

and

$$\begin{array}{}\text{(16)}& \underset{\mathit{\theta }\to \mathit{\pi }}{lim}\sin \mathit{\theta }{\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left(\mathit{\theta }\right)=\mathrm{0}.\end{array}$$

This ensures that Eq. (14) takes the simpler form

$$\begin{array}{}\text{(17)}& {\displaystyle }T\left(\mathit{\xi }\right){\displaystyle }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h(T,\mathit{\theta }).\end{array}$$

Using this, we may express the solution to Eq. (3) as

$$\begin{array}{}\text{(18)}& T\left(\mathit{\xi }\right)=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h_{\mathrm{1}}\left(\mathit{\theta }\right)\end{array}$$

for case (1) and

$$\begin{array}{}\text{(19)}& T\left(\mathit{\xi }\right)=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h_{\mathrm{2}}\left(\mathit{\theta }\right)\end{array}$$

for case (2), where $h_{\mathrm{1}}=\mathit{\eta }s\left(\mathit{\theta }\right)(\mathrm{1}-a_{\mathrm{1}})-\mathit{\alpha }$ and $h_{\mathrm{2}}=\mathit{\eta }s\left(\mathit{\theta }\right)(\mathrm{1}-a_{\mathrm{2}})-\mathit{\alpha }$.

2.2 Partial ice cover

For solutions to Eq. (1) where the zonal mean temperature profile is not strictly above or below the critical temperature, T_s, the surface will have a partial ice cover analogous to the Earth's current climate state. Critical latitudes, denoted as ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}$ and ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}$, mark the transitions between ice and water coverage. Assuming that T remains continuous across the critical latitudes, the non-dimensional temperature at these latitudes must necessarily be $T\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)=T\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)=-\mathrm{1}$. Given that the radiation distribution s(θ) prescribes an incoming radiation maximum at the equator, the ice cover must be centered at the poles (as illustrated in Fig. 1 but without the continent), and it is sensible to partition the domain into the following subdomains:

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }\mathit{\theta }\in (\mathrm{0},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}),\text{(21)}& {\displaystyle }& {\displaystyle }\mathit{\theta }\in ({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}})\end{array}$$

and

$$\begin{array}{}\text{(22)}& \mathit{\theta }\in ({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},\mathit{\pi }),\end{array}$$

such that

$$\begin{array}{}\text{(23)}& \mathcal{L}T=h_{\mathrm{1}}\end{array}$$

in the region in Eq. (21) and

$$\begin{array}{}\text{(24)}& \mathcal{L}T=h_{\mathrm{2}}\end{array}$$

in the regions in Eqs. (20) and (22).

Figure 1 Schematic of the domain $\mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$ of Eq. (1) for a planet with a zonally symmetric continent and partial ice cover. The extent of the continent, i.e., ${\mathit{\theta }}_{l_{\mathrm{1}}}$ and ${\mathit{\theta }}_{l_{\mathrm{2}}}$, is determined by Eq. (34). The ice caps extend from θ=0 to the critical latitude $\mathit{\theta }={\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}$ and from $\mathit{\theta }={\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}$ to θ=π.

The relation in Eq. (10) is applied in these regions, and we get

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathit{\xi }\right)& {\displaystyle }={\int }_{\mathrm{0}}^{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h_{\mathrm{2}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }+K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},\mathit{\xi })\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\\ {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{-}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }),\end{array}\text{(26)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathit{\xi }\right)& {\displaystyle }={\int }_{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}}^{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h_{\mathrm{1}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }+K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},\mathit{\xi })\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{-}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi })\\ {\displaystyle }& {\displaystyle }-K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},\mathit{\xi })\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\\ {\displaystyle }& {\displaystyle }-\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{+}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi })\end{array}\end{array}$$

and

$$\begin{array}{}\text{(27)}& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathit{\xi }\right)& {\displaystyle }={\int }_{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}}^{\mathit{\pi }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })h_{\mathrm{2}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }-K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},\mathit{\xi })\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }-\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{+}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi })\end{array}\end{array}$$

in the regions in Eqs. (20), (21) and (22), respectively. The solution to Eq. (3) within the three subdomains may be expressed through the relations in Eqs. (25)–(27), given the points ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}$ and ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}$, as well as the spatial derivative of T at these points. These unknown boundary values are determined by solving the following system of boundary integral equations, obtained by letting ξ approach the boundaries of the three subdomains in Eqs. (25)–(27):

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathrm{0}\right)& {\displaystyle }={\int }_{\mathrm{0}}^{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathrm{0})h_{\mathrm{2}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }+K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},\mathrm{0})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\\ {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\underset{\mathit{\xi }\to {\mathrm{0}}^{+}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{-}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }),\end{array}\text{(29)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }-\mathrm{1}& {\displaystyle }={\int }_{\mathrm{0}}^{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}})h_{\mathrm{2}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }+K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\\ {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\underset{\mathit{\xi }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{-}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{-}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }),\end{array}\text{(30)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }-\mathrm{1}& {\displaystyle }={\int }_{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}}^{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}})h_{\mathrm{1}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }+K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\underset{\mathit{\xi }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{+}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{-}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi })\\ {\displaystyle }& {\displaystyle }-K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\\ {\displaystyle }& {\displaystyle }-\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\underset{\mathit{\xi }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{+}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{+}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }),\end{array}\text{(31)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }-\mathrm{1}& {\displaystyle }={\int }_{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}}^{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}})h_{\mathrm{1}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }+K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\underset{\mathit{\xi }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{-}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{-}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi })\\ {\displaystyle }& {\displaystyle }-K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\\ {\displaystyle }& {\displaystyle }-\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)\underset{\mathit{\xi }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{-}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}^{+}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }),\end{array}\text{(32)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }-\mathrm{1}& {\displaystyle }={\int }_{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}}^{\mathit{\pi }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}})h_{\mathrm{2}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }-K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}})\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }-\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\underset{\mathit{\xi }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{+}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{+}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }),\end{array}\text{(33)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }T\left(\mathit{\pi }\right)& {\displaystyle }={\int }_{{\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}}^{\mathit{\pi }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\pi })h_{\mathrm{2}}\left(\mathit{\theta }\right)\\ {\displaystyle }& {\displaystyle }-K({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}},\mathit{\pi })\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right){\displaystyle \frac{\partial T}{\partial \mathit{\theta }}}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }-\sin \left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)\underset{\mathit{\xi }\to {\mathit{\pi }}^{-}}{lim}\underset{\mathit{\theta }\to {\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}^{+}}{lim}{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}(\mathit{\theta },\mathit{\xi }).\end{array}\end{array}$$

The system of equations (Eqs. 28–33) can be solved for T(0), T(π), $\frac{\partial T}{\partial \mathit{\theta }}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}\right)$, $\frac{\partial T}{\partial \mathit{\theta }}\left({\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}\right)$, ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}$ and ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}$ through a combination of both analytical and numerical methods. Finding ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{1}}}$ and ${\mathit{\theta }}_{{\mathrm{c}}_{\mathrm{2}}}$ ultimately requires a root search, as the integrals cannot be evaluated analytically. However, these values may be approximated to a high precision. For the root search, we used a Broyden's method, as computing the Jacobian is expensive for solutions with multiple critical latitudes. Note that the system of equations (Eqs. 28–33) can have more than one solution, indicating multiple equilibria.

2.3 With a continent

The method presented may be extended to include one or more zonally symmetric continents. The following analysis includes one such continent with a meridional extent of $l=\frac{\mathit{\pi }}{\mathrm{4}}$, stretching from latitude ${\mathit{\theta }}_{l_{\mathrm{1}}}$ to latitude ${\mathit{\theta }}_{l_{\mathrm{2}}}$. Figure 1 illustrates a planet with a continent of this kind. The analysis was repeated three times for three different continent configurations,

$$\begin{array}{}\text{(34)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathit{\theta }}_{l_{\mathrm{1}}}& {\displaystyle }={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}-{\displaystyle \frac{l}{\mathrm{2}}}-\mathit{\epsilon },\\ {\displaystyle }{\mathit{\theta }}_{l_{\mathrm{2}}}& {\displaystyle }={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}+{\displaystyle \frac{l}{\mathrm{2}}}-\mathit{\epsilon }.\end{array}\end{array}$$

where ε=0, ε=0.1 and ε=0.5. Parameter values on the continent may be altered to distinguish land from ocean and better capture the thermal response of the lithosphere. Here, the heat capacity C, the critical temperature for ice formation T_s and the albedo a(T) were altered on the part of the domain corresponding to the continent. This has the effect of changing the ice dynamics and subsequently the ice-albedo feedback on the continent. The extent of the continent l is kept constant, and no ice–ocean–land feedback is included in the model. Mathematical details on the application of the presented method to model configurations with a continent are omitted for brevity. Instead, we provide some results of our analysis. Interested readers are referred to Samuelsberg and Jakobsen (2023) for a full derivation of these solutions. The multiple branch structure of the model is displayed in Fig. 2 through bifurcation diagrams of the system with a continent configuration as in Eq. (34), where ε=0, ε=0.1 and ε=0.5. The control parameter is the scaled TSI Q. Stability properties of the stationary solutions are assessed using the numerical perturbation scheme outlined in Appendix B.

Figure 2 Bifurcation diagrams. Annual mean equilibrium surface temperatures plotted against the control parameter Q for model configurations with a continent as in Eq. (34), where ε=0, ε=0.1 and ε=0.5. Model parameters are those in Table 1. Solid lines indicate stable solutions and dotted lines are unstable solutions. The upper, stable branches in the bifurcation diagrams contain ice-free solutions, and the lower, stable branches contain solutions with a full ice cover (Snowball Earth states). Intermediate branches contain solutions with a partial ice cover.

Table 1 Model parameters used in the presented work. A, B, C, C_land and T_s are taken from North et al. (1981). D is taken from Kaper and Engler (2013), and s(θ) is taken from McGehee and Lehman (2012).

3 Discussion

We have applied the BIM to solve the stationary form of a North-type EBM with a zonally symmetric continent in three different configurations. The introduction of a continent resulted in the emergence of new equilibrium states. EBMs have a rich multiple solution structure (North, 1990). In zero-dimensional models that incorporate the ice-albedo feedback, three solutions exist; as one extends to one-dimensional models, the multiple branch structure becomes more complicated. North showed that there exist up to five solutions for a range of TSI values in the globally averaged model, one of which is the famously unstable small-ice-cap solution (North, 1984). From Fig. 2a, it is evident that up to seven equilibria may exist for a range of Q values in model configurations with a continent. The continent is initially placed in a north–south symmetrical configuration to investigate how the system is affected by symmetry. The meridional symmetry is subsequently violated by increasing ε. Although the number of equilibria for any given Q is at most seven for ε=0 and ε=0.1, the range over which seven equilibria can exist is reduced as ε is increased and has disappeared for ε=0.5. Prescribing the system with inherently symmetrical boundary conditions, i.e., ε=0, evidently introduces a very fine dynamic as seen in the looping branches in the bifurcation diagram of Fig. 2a, which disappears for non-zero ε. Furthermore, equilibria with a more complicated ice distribution, characterized by more than two critical latitudes, only appear for ε=0. Despite this, the range of TSI over which stable partial-ice-cover solutions can exist, i.e., the intermediate branches in Fig. 2, is much shorter for ε=0. The range of stable, intermediate branches is longest for ε=0.5. Additionally, for ε=0.5, there is a large range of TSI values over which bi-stability occurs between the ice-free solution and the partial-ice-cover solution, as seen in Fig. 2c. These observations are notable because, although it is well established that the climate has changed throughout geological history, the role of the land–sea distribution is not fully understood (Fluteau, 2003).

The robustness of the method was tested by allowing for parameter discontinuities at the land–sea boundaries. Mengel et al. (1988) and North and Kim (2017) have discussed the application of Fourier–Legendre series in EBMs with parameter discontinuities at the continent edges: a discontinuity in the albedo and heat capacity parameter C causes a solution expressed through Legendre modes to converge slowly. Presumably, similar discontinuities in other parameters will have the same effect. Moreover, North (1975a) discussed the potential for changing parameter values on finite zonal strips using spectral methods, but they did not address how this may result in several ice edges. The dynamics of EBMs are highly sensitive to model parameters (Soldatenko and Colman, 2019), and the introduction of a continent to the model can result in some unusual ice distributions. Figure 3 shows an equilibrium solution with six critical latitudes and two ice belts on the continent. Studying several critical latitudes is a natural extension of the BIM and follows the same general procedure. A similar ice belt has been observed by changing the obliquity of the model using the spectral method (Rose et al., 2017). For high obliquities, the traditional ice distribution of the two-critical-latitude solution is inverted, allowing for a similar analysis to the classical partial-ice-cover states, without additional ice edges. The ice belts overlaying the continent form a striped pattern, a phenomenon also observed in the related one-dimensional Daisyworld model with diffusion (Adams et al., 2003; Alberti et al., 2015). Adams et al. (2003) found that the daisies never coexisted; instead the equilibrium solution is characterized by zonal bands of single-species colonies reminiscent of Fig. 3. A similar striped pattern of daisy coverage was reported by Alberti et al. (2015) in a one-dimensional Daisyworld model with diffusion and a greenhouse effect.

Figure 3 An equilibrium solution (unstable) to Eq. (1) with a continent as in Eq. (34), where ε=0, Q=294 W m⁻² and other parameters are those in Table 1. Green vertical lines mark the continent borders, and blue vertical lines mark critical latitudes. The red horizontal lines mark the critical temperatures for ice formation.

Although the BIM is an effective approach for solving North-type EBMs, it has certain limitations. A step function albedo is used here, and the method is easily extended to any latitude dependence for the albedo. However, this method, like the spectral method, is limited to a step function-like temperature response at the ice edge. An arbitrary temperature-dependent albedo on either side of the ice edge renders the energy balance equation unsolvable through the presented method. An additional drawback of the presented method is that a root search is required to find the critical latitudes, ${\mathit{\theta }}_{{\mathrm{c}}_i}$. For solutions with a low number of critical latitudes, this poses no significant challenge. However, as the number of critical latitudes increases, so does the complexity of the root search and the necessity for a good starting point in the iteration. For solutions located in branches connected to the ice-free branch and the Snowball Earth branch in the bifurcation diagram, this is generally not a problem. Since these extreme solutions always exist within certain parameter regimes, we identify the bifurcation points and update the ansatz accordingly, repeating this process until the upper and lower branches connect. This approach ensures sufficient initial points for the root searches, allowing all coexisting states in connected branches to be identified. Additionally, for solutions with one, two or three critical latitudes, a graphical solution to the system of boundary integral equations (Eqs. 28–33) is possible, ensuring that any isolated branches are also detected. However, detecting isolated branches with multiple critical latitudes requires a more numerically expensive approach. Furthermore, the method can become tedious for complicated geographies and a high number of critical latitudes.

The BIM has a wide range of potential applications. While the model studied here only includes an ice-albedo feedback mechanism, the BIM can be generalized to solve one-dimensional EBMs with additional albedo feedbacks. A vegetation feedback has previously been introduced within an EBM framework for zero-dimensional models (Rombouts and Ghil, 2015; Alberti et al., 2018) and related one-dimensional models (Wood et al., 2008; Adams et al., 2003; Alberti et al., 2015). The method we have presented is applicable to North-type EBMs where the model is piecewise linear on subdomains, specifically if the explicit temperature dependence in the integrand of Eq. (10) can be eliminated through an appropriate partitioning of the domain. For instance, a vegetation feedback mechanism can be implemented by adding vegetation on the continent for temperatures within a certain growth regime, where the vegetation is modeled by altering the albedo. The BIM can be directly applied to North-type models with such simplified vegetation responses. Furthermore, including additional spatial dependence in the albedo represents another natural application of the BIM. For example, in studies involving North-type EBMs related to Snowball Earth events, alternative albedo parameterizations (Abbot et al., 2011) or spatial dependence for other model parameters can be easily implemented with the BIM. The ability to effectively handle spatially dependent model parameters also makes the BIM an appropriate tool for studying fragmented tipping in North-type EBMs (Bastiaansen et al., 2022).

4 Conclusions

In this paper, we have presented an analytical method for solving North-type EBMs. Solutions are expressed through explicit expressions, readily obtainable from quadrature methods. The presented method has some notable advantages compared to other analytical methods, e.g., North (1975a), for solving energy balance equations of this kind. It does not rely on truncating series expansions. Furthermore, the method remains straightforward and, computationally speaking, very fast, even for problems with partial land–sea geographies and parameter discontinuities at the boundaries separating land and sea. In addition, the BIM offers a formulaic framework for handling equilibrium solutions with several critical latitudes, ${\mathit{\theta }}_{{\mathrm{c}}_i}$, and unconventional ice distributions.

The development of new analytical methods of studying EBMs is motivated by the recognition of EBMs as useful tools for researchers in a variety of fields. The inherent simplicity of EBMs, characterized by few parameters, renders them particularly suitable for certain aspects of paleoclimatology (North and Kim, 2017; Abbot et al., 2011; Widiasih et al., 2024) and planetary science (Rose et al., 2017), where poorly constrained parameters and a diverse set of planetary conditions are frequently encountered. The BIM represents a systematic approach to solving North-type EBMs and excels under unconventional model configurations, particularly where emerging solutions describe climate states markedly different from the prevailing state of Earth. Analytical investigations of conceptual models continue to provide a valuable testing ground for ideas in climate science and insights into the complex dynamics involved as one ascends the climate model hierarchy.

Appendix A: Finding a Green's function

In this section, we find a Green's function, K(θ,ξ), for the operator in Eq. (4). A Green's function is any solution to Eq. (8). The two defining properties of the Dirac delta function are the following:

1. For any surface of interest, S, we must have

   $$\begin{array}{}\text{(A1)}& \underset{S}{\int }\mathrm{d}A{\mathit{\delta }}_{\mathit{\xi }}=\mathrm{1}.\end{array}$$

2. For any function, f(x), defined on S, we must have

   $$\begin{array}{}\text{(A2)}& \underset{S}{\int }\mathrm{d}A{\mathit{\delta }}_{\mathit{\xi }}f=f\left(\mathit{\xi }\right).\end{array}$$

For the line $\mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$ along the surface of a sphere with radius R, it can be shown that the term

$$\begin{array}{}\text{(A3)}& {\mathit{\delta }}_{\mathit{\xi }}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathit{\delta }(\mathit{\theta }-\mathit{\xi })}{\mathrm{2}\mathit{\pi }{R}^{\mathrm{2}}\sin \mathit{\theta }}},\end{array}$$

where δ(θ−ξ) is the usual delayed Dirac delta function on the line, will ensure that Eqs. (A1) and (A2) are satisfied. It is convenient to scale the Green's function we are seeking by a factor such that the right-hand side of Eq. (A3) becomes unity and that Eq. (8) becomes

$$\begin{array}{}\text{(A4)}& \mathcal{L}K={\displaystyle \frac{\mathit{\delta }(\mathit{\theta }-\mathit{\xi })}{\sin \mathit{\theta }}}.\end{array}$$

We will demand that Green's function is continuous across θ=ξ; therefore, we must have

$$\begin{array}{}\text{(A5)}& \underset{\mathit{\theta }\to {\mathit{\xi }}^{+}}{lim}K(\mathit{\theta },\mathit{\xi })=\underset{\mathit{\theta }\to {\mathit{\xi }}^{-}}{lim}K(\mathit{\theta },\mathit{\xi }).\end{array}$$

Integrating Eq. (A4) over a small interval centered on θ=ξ, we get

$$\begin{array}{}\text{(A6)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\int }_{\mathit{\xi }-\mathit{\epsilon }}^{\mathit{\xi }+\mathit{\epsilon }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }\mathit{\{}-{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}& {\displaystyle }\left(\sin \mathit{\theta }{\displaystyle \frac{\partial K}{\partial \mathit{\theta }}}\right)+\mathit{\beta }K\mathit{\}}=\\ {\displaystyle }& {\displaystyle }{\int }_{\mathit{\xi }-\mathit{\epsilon }}^{\mathit{\xi }+\mathit{\epsilon }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }{\displaystyle \frac{\mathit{\delta }(\mathit{\theta }-\mathit{\xi })}{\sin \mathit{\theta }}},\end{array}\text{(A7)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }-{\int }_{\mathit{\xi }-\mathit{\epsilon }}^{\mathit{\xi }+\mathit{\epsilon }}\mathrm{d}\mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}& {\displaystyle }\left(\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}K(\mathit{\theta },\mathit{\xi })\right)\\ {\displaystyle }& {\displaystyle }+\mathit{\beta }{\int }_{\mathit{\xi }-\mathit{\epsilon }}^{\mathit{\xi }+\mathit{\epsilon }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })=\mathrm{1},\end{array}\text{(A8)}& {\displaystyle }& {\displaystyle }-\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}K(\mathit{\theta },\mathit{\xi }){|}_{\mathit{\xi }-\mathit{\epsilon }}^{\mathit{\xi }+\mathit{\epsilon }}+\mathit{\beta }\underset{\mathit{\xi }-\mathit{\epsilon }}{\overset{\mathit{\xi }+\mathit{\epsilon }}{\int }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }K(\mathit{\theta },\mathit{\xi })=\mathrm{1}.\end{array}$$

Letting ε→0, we must have

$$\begin{array}{}\text{(A9)}& \underset{\mathit{\theta }\to {\mathit{\xi }}^{+}}{lim}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}K(\mathit{\theta },\mathit{\xi })-\underset{\mathit{\theta }\to {\mathit{\xi }}^{-}}{lim}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}K(\mathit{\theta },\mathit{\xi })=-{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\xi }}}.\end{array}$$

At θ≠ξ we evidently have

$$\begin{array}{}\text{(A10)}& \mathcal{L}K(\mathit{\theta },\mathit{\xi })=\mathrm{0}.\end{array}$$

We can therefore conclude that K(θ,ξ) must satisfy the necessary conditions of Eqs. (A5), (A9) and (A10). In order to find a Green's function that solves Eq. (A4), we need to find a basis of solutions for an equation of the form

$$\begin{array}{}\text{(A11)}& -{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}y\left(\mathit{\theta }\right)\right)+\mathit{\beta }y\left(\mathit{\theta }\right)=\mathrm{0}.\end{array}$$

Introducing a change of variables, x=cos θ, and a function, u, such that u(cos θ)=y(θ), it can be shown that Eq. (A11) can be written in the form

$$\begin{array}{}\text{(A12)}& (\mathrm{1}-x^{\mathrm{2}}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}u\left(x\right)+\mathrm{2}x{\displaystyle \frac{\partial }{\partial x}}u\left(x\right)-\mathit{\beta }u\left(x\right)=\mathrm{0}.\end{array}$$

Let λ be a number such that $-\mathit{\beta }=\mathit{\lambda }(\mathit{\lambda }+\mathrm{1})$. We may now write Eq. (A12) as a Legendre equation:

$$\begin{array}{}\text{(A13)}& (\mathrm{1}-x^{\mathrm{2}}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}u\left(x\right)+\mathrm{2}x{\displaystyle \frac{\partial }{\partial x}}u\left(x\right)+\mathit{\lambda }(\mathit{\lambda }+\mathrm{1})u\left(x\right)=\mathrm{0},\end{array}$$

which for some arbitrary real or complex value λ will have the known basis of solutions {P_λ, Q_λ}. We can therefore use the basis {P_λ(cos θ), Q_λ(cos θ)}, where $\mathit{\lambda }=\frac{\mathrm{1}}{\mathrm{2}}\left(\sqrt{\mathrm{1}-\mathrm{4}\mathit{\beta }}-\mathrm{1}\right)$, to construct the general solution to Eq. (A10):

$$\begin{array}{}\text{(A14)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }K(\mathit{\theta },\mathit{\xi })=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}a\left(\mathit{\xi }\right)P_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)+b\left(\mathit{\xi }\right)Q_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right),& \mathit{\theta }>\mathit{\xi },\\ c\left(\mathit{\xi }\right)P_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)+d\left(\mathit{\xi }\right)Q_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right),& \mathit{\theta }<\mathit{\xi }.\end{array}\right.\end{array}\end{array}$$

The coefficients $a\left(\mathit{\xi }\right),b\left(\mathit{\xi }\right),c\left(\mathit{\xi }\right)$ and d(ξ) can be determined through the conditions in Eqs. (A5) and (A9). Any choice of these coefficients satisfying Eqs. (A5) and (A9) will give a Green's function for the operator in Eq. (4). However, it makes sense for us to seek a Green's function that is non-singular in the domain $\mathit{\theta }\in [\mathrm{0},\mathit{\pi }]$: we want to develop a set of conditions for the coefficients a(ξ), b(ξ), c(ξ) and d(ξ) to ensure that Green's function (Eq. A14) is non-singular when θ→0 and θ→π. Let

$$\begin{array}{}\text{(A15)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{K}_{+}(\mathit{\theta },\mathit{\xi })& {\displaystyle }=a\left(\mathit{\xi }\right)P_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)+b\left(\mathit{\xi }\right)Q_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right),\\ {\displaystyle }{K}_{-}(\mathit{\theta },\mathit{\xi })& {\displaystyle }=c\left(\mathit{\xi }\right)P_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)+d\left(\mathit{\xi }\right)Q_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)\end{array}\end{array}$$

such that

$$\begin{array}{}\text{(A16)}& K(\mathit{\theta },\mathit{\xi })=\left\{\begin{array}{ll}{K}_{+}(\mathit{\theta },\mathit{\xi }),& \mathit{\theta }>\mathit{\xi },\\ K_{-}(\mathit{\theta },\mathit{\xi }),& \mathit{\theta }<\mathit{\xi }.\end{array}\right.\end{array}$$

Using a computer algebra system, we find a series expansion of K₊ around θ=π, and we recognize that there is a term in this expansion containing log (π−θ) with a coefficient $c_{+\mathrm{0}}\left(a\right(\mathit{\xi }),b(\mathit{\xi }\left)\right)$. We want to ensure that K₊ is non-singular at θ=π and therefore demand that

$$\begin{array}{}\text{(A17)}& c_{+\mathrm{0}}\left(a\right(\mathit{\xi }),b(\mathit{\xi }\left)\right)=\mathrm{0}\forall \mathit{\xi }\in [\mathrm{0},\mathit{\pi }].\end{array}$$

Similarly, we find a series expansion of K₋ around θ=0. In this expansion, there is a term containing log (θ) with a coefficient $c_{-\mathrm{0}}\left(c\right(\mathit{\xi }),d(\mathit{\xi }\left)\right)$, and we demand that

$$\begin{array}{}\text{(A18)}& c_{-\mathrm{0}}\left(c\right(\mathit{\xi }),d(\mathit{\xi }\left)\right)=\mathrm{0}\forall \mathit{\xi }\in [\mathrm{0},\mathit{\pi }].\end{array}$$

Solving the system of equations (Eqs. A5, A9, A17 and A18) for a(ξ), b(ξ), c(ξ) and d(ξ), we find the following Green's function:

$$\begin{array}{ll}\text{(A19)}& {\displaystyle }& {\displaystyle }K(\mathit{\theta },\mathit{\xi })={\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\frac{P_{\mathit{\lambda }}\left(\cos \mathit{\xi }\right)\left(\mathit{\pi }\cot \right(\mathit{\pi }\mathit{\lambda }\left)P_{\mathit{\lambda }}\right(\cos \mathit{\theta })-\mathrm{2}{Q}_{\mathit{\lambda }}(\cos \mathit{\theta }\left)\right)}{\mathrm{2}(\mathrm{1}+\mathit{\lambda })\left(P_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }+\mathrm{1}}\right(\cos \mathit{\xi })-P_{\mathit{\lambda }+\mathrm{1}}(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)\right)},& \mathit{\theta }>\mathit{\xi },\\ \frac{P_{\mathit{\lambda }}\left(\cos \mathit{\theta }\right)\left(\mathit{\pi }\cot \right(\mathit{\pi }\mathit{\lambda }\left)P_{\mathit{\lambda }}\right(\cos \mathit{\xi })-\mathrm{2}{Q}_{\mathit{\lambda }}(\cos \mathit{\xi }\left)\right)}{\mathrm{2}(\mathrm{1}+\mathit{\lambda })\left(P_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }+\mathrm{1}}\right(\cos \mathit{\xi })-P_{\mathit{\lambda }+\mathrm{1}}(\cos \mathit{\xi }\left)Q_{\mathit{\lambda }}\right(\cos \mathit{\xi }\left)\right)},& \mathit{\theta }<\mathit{\xi }.\end{array}\right.\end{array}$$

Green's function (Eq. A19) is non-singular and bounded at θ=0 and θ=π. The derivative of Green's function ( Eq. A19) is also bounded at the boundary and tends to zero $\forall \mathit{\xi }\in [\mathrm{0},\mathit{\pi }]$.

Appendix B: Stability analysis

In this section, we are going to test the stability of the stationary solutions found using the BIM. Applying the notation from Sect. 2, we may write the time-dependent form of Eq. (1) as

$$\begin{array}{}\text{(B1)}& \mathit{\gamma }{\partial }_{t}T+\mathcal{L}T=h(T,\mathit{\theta }),\end{array}$$

where $\mathit{\gamma }=\frac{C}{t_{\mathrm{0}}D}$. Stationary solutions are denoted $T_{\mathrm{0}}=T(\mathit{\theta },t=\mathrm{0})$, such that ℒT₀=h. We wish to investigate whether a small perturbation, δ, away from the equilibrium, T₀, will grow in time. The perturbed solution, $T(\mathit{\theta },t)=T_{\mathrm{0}}\left(\mathit{\theta }\right)+\mathit{\delta }(\mathit{\theta },t)$, inserted into Eq. (B1) yields

$$\begin{array}{}\text{(B2)}& \mathit{\gamma }{\partial }_t\mathit{\delta }+\mathcal{L}{T}_{\mathrm{0}}+\mathcal{L}\mathit{\delta }=h(T_{\mathrm{0}}+\mathit{\delta },\mathit{\theta }).\end{array}$$

A first-order expansion around (T₀,θ) of the right-hand side may be expressed as

$$\begin{array}{}\text{(B3)}& h(T_{\mathrm{0}}+\mathit{\delta },\mathit{\theta })\approx h(T_{\mathrm{0}},\mathit{\theta })+h_T(T_{\mathrm{0}},\mathit{\theta })\mathit{\delta }.\end{array}$$

Here we let the term

$$\begin{array}{}\text{(B4)}& a\left(T\right)=a_{\mathrm{1}}+{\displaystyle \frac{a_{\mathrm{2}}-a_{\mathrm{1}}}{\mathrm{2}}}(\mathrm{1}+\tanh (-\mathit{\sigma }(T+\mathrm{1})\left)\right),\end{array}$$

where the slope parameter σ=50, be a smooth function replicating the behavior of the step function albedo such that the derivative $h_T(T_{\mathrm{0}},\mathit{\theta })=\frac{\partial h}{\partial T}(T_{\mathrm{0}},\mathit{\theta })$ may be found analytically for a given T₀. By substituting Eq. (B3) into Eq. (B2), we get

$$\begin{array}{}\text{(B5)}& {\partial }_t\mathit{\delta }=\mathcal{H}\mathit{\delta },\end{array}$$

where

$$\begin{array}{}\text{(B6)}& \mathcal{H}(\cdot )={\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }}}\left[h_T(T_{\mathrm{0}},\mathit{\theta })(\cdot )-\mathcal{L}(\cdot )\right].\end{array}$$

Suppose that the perturbation δ has the form

$$\begin{array}{}\text{(B7)}& \mathit{\delta }(\mathit{\theta },t)=e^{\mathit{\lambda }t}{\mathit{\delta }}_{\mathrm{0}}\left(\mathit{\theta }\right).\end{array}$$

This turns Eq. (B5) into the following eigenvalue problem:

$$\begin{array}{}\text{(B8)}& \mathit{\lambda }{\mathit{\delta }}_{\mathrm{0}}=\mathcal{H}{\mathit{\zeta }}_{\mathrm{0}}.\end{array}$$

Real and positive λ values will evidently cause the perturbation to grow exponentially, resulting in an unstable but stationary solution T₀. The eigenvalues were subsequently approximated using a numerical scheme that solves the associated eigenvalue problem:

$$\begin{array}{}\text{(B9)}& \mathit{\lambda }{\mathit{\delta }}_{\mathrm{0}}=\mathbf{H}{\mathit{\delta }}_{\mathrm{0}},\end{array}$$

where the set of linear equations

$$\begin{array}{}\text{(B10)}& \mathit{\{}\mathit{\lambda }{\mathit{\delta }}_{\mathrm{0}}^i={\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }}}[h_T(T_{\mathrm{0}}^i,{\mathit{\theta }}_i){\mathit{\delta }}_{\mathrm{0}}^i-\widehat{\mathcal{L}}{\mathit{\delta }}_{\mathrm{0}}^i]{\mathit{\}}}_{i=\mathrm{0}}^N\end{array}$$

gives rise to the coefficient matrix H. Here $T_{\mathrm{0}}^i$ is the stationary solution evaluated on a uniform spatial grid, and $\widehat{\mathcal{L}}(\cdot )$ is a finite difference approximation of the differential operator ℒ. It can be shown that a second-order centered difference approximation for a smooth function δ₀, evaluated on a uniform grid θ_i, where ${\mathit{\delta }}_{\mathrm{0}}^i={\mathit{\delta }}_{\mathrm{0}}\left({\mathit{\theta }}_i\right)$, is

$$\begin{array}{ll}\text{(B11)}& {\displaystyle }\widehat{\mathcal{L}}{\mathit{\delta }}_{\mathrm{0}}^i& {\displaystyle }=\mathit{\beta }{\mathit{\delta }}_{\mathrm{0}}^i{\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{2}({\mathit{\delta }}_{\mathrm{0}}^{i-\mathrm{1}}-\mathrm{2}{\mathit{\delta }}_{\mathrm{0}}^i+{\mathit{\delta }}_{\mathrm{0}}^{i+\mathrm{1}})+({\mathit{\delta }}_{\mathrm{0}}^{i-\mathrm{1}}-\mathrm{4}{\mathit{\delta }}_{\mathrm{0}}^i+\mathrm{3}{\mathit{\delta }}_{\mathrm{0}}^{i+\mathrm{1}})d\mathit{\theta }\cot {\mathit{\theta }}_i}{\mathrm{2}d{\mathit{\theta }}^{\mathrm{2}}}}.\end{array}$$

At either end of the grid, forward and backward approximations are needed. These are

$$\begin{array}{ll}\text{(B12)}& {\displaystyle }{\widehat{\mathcal{L}}}_f{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{0}}& {\displaystyle }=\mathit{\beta }{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{0}}{\displaystyle }& {\displaystyle }+{\displaystyle \frac{-\mathrm{2}({\mathit{\delta }}_{\mathrm{0}}^{\mathrm{0}}-\mathrm{2}{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{1}}+{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{2}})+(\mathrm{5}{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{0}}-\mathrm{8}{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{1}}+\mathrm{3}{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{2}})d\mathit{\theta }\cot {\mathit{\theta }}_{\mathrm{0}}}{\mathrm{2}d{\mathit{\theta }}^{\mathrm{2}}}}\end{array}$$

and

$$\begin{array}{ll}\text{(B13)}& {\displaystyle }{\widehat{\mathcal{L}}}_b{\mathit{\delta }}_{\mathrm{0}}^N& {\displaystyle }=\mathit{\beta }{\mathit{\delta }}_{\mathrm{0}}^N{\displaystyle }& {\displaystyle }+{\displaystyle \frac{-\mathrm{2}({\mathit{\delta }}_{\mathrm{0}}^{N-\mathrm{2}}-\mathrm{2}{\mathit{\delta }}_{\mathrm{0}}^{N-\mathrm{1}}+{\mathit{\delta }}_{\mathrm{0}}^N)+({\mathit{\delta }}_{\mathrm{0}}^{N-\mathrm{2}}-{\mathit{\delta }}_{\mathrm{0}}^N)d\mathit{\theta }\cot {\mathit{\theta }}_N}{\mathrm{2}d{\mathit{\theta }}^{\mathrm{2}}}}\end{array}$$

for the forward and backward approximations, respectively. Furthermore, as t grows the perturbation in Eq. (B7) must adhere to the same constraints as the solution, i.e., the boundary conditions in Eqs. (15) and (16). A discrete formulation of these is ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{0}}={\mathit{\delta }}_{\mathrm{0}}^{\mathrm{2}}$ and ${\mathit{\delta }}_{\mathrm{0}}^N={\mathit{\delta }}_{\mathrm{0}}^{N-\mathrm{2}}$, which we ensure is enforced. We build the matrix H and examine the associated eigenvalues for a large number of points in the bifurcation diagram. Stability properties are subsequently inferred from the ensemble of stationary solutions within the same branch. The presented stability analysis agrees with the slope-stability theorem put forth by Cahalan and North (1979).