Climate bifurcations in a Schwarzschild equation model of the Arctic atmosphere
 ^{1}Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark
 ^{2}Dept. of Mathematics & Statistics, University of Guelph, Guelph, Canada
 ^{3}Faculty of Science, Ontario Tech University, Oshawa, Canada
 ^{1}Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark
 ^{2}Dept. of Mathematics & Statistics, University of Guelph, Guelph, Canada
 ^{3}Faculty of Science, Ontario Tech University, Oshawa, Canada
Correspondence: Allan R. Willms (awillms@uoguelph.ca)
Hide author detailsCorrespondence: Allan R. Willms (awillms@uoguelph.ca)
A column model of the Arctic atmosphere is developed including the nonlinear positive feedback responses of surface albedo and water vapour to temperature. The atmosphere is treated as a grey gas and the flux of longwave radiation is governed by the twostream Schwarzschild equations. Water vapour concentration is determined by the Clausius–Clapeyron equation. Representative concentration pathways (RCPs) are used to model carbon dioxide concentrations into the future. The resulting 9D twopoint boundary value problem is solved under various RCPs and the solutions analysed. The model predicts that under the highest carbon pathway, the Arctic climate will undergo an irreversible bifurcation to a warm steady state, which would correspond to annually icefree conditions. Under the lowest carbon pathway, corresponding to very aggressive carbon emission reductions, the model exhibits only a mild increase in Arctic temperatures. Under the two intermediate carbon pathways, temperatures increase more substantially, and the system enters a region of bistability where external perturbations could possibly cause an irreversible switch to a warm, icefree state.
Climate change is causing rapid temperature increases in the polar regions. A fundamental question is whether these temperature increases are reversible. If humanity fails to prevent a substantial warming of the planet in the next few decades, which is appearing to be more and more likely, will it be possible in the future to reverse our effects on climate enough to restore lower temperatures, or will we have passed a tipping point beyond which return to the present state is impossible? We address this question in particular for the Arctic, where the observed climate change is the most dramatic.
The Earth's climate is an extremely complex system. Modelling efforts range from simple models attempting to isolate the most pertinent features to very complicated numerical models trying to capture as many details as possible. The model presented here is close to the simple end of this spectrum, although not as simple as some, in that it is a 9D nonlinear twopoint boundary value problem. The advantage of relatively simple models is that they allow more direct analysis of cause and effect, which is often obscured in highly complicated models.
The term “tipping point” is used by different researchers in various ways; see Russill (2015) and Lenton et al. (2008) for some definitions and discussion of the term. In all cases, however, tipping points are associated with large qualitative changes in a system due to relatively small changes in the parameters or “forcings” that drive the system. In the present paper, tipping points arise as a result of saddlenode and cusp bifurcations in the mathematical model. The mathematical theory of bifurcations is welldeveloped (Kuznetsov, 2004) and employed here. Figure 1 illustrates the typical behaviour associated with these bifurcations. In Fig. 1a there are two saddlenode bifurcations resulting in a parameter interval of bistability; that is, two stable states coexist for an interval of μ values. If the system is on the lower branch of stable states, then, as μ increases through a critical value μ_{crit}, there is an abrupt jump to the upper branch of stable states. In contrast, as μ decreases, the jump back to the lower state does not occur until a much smaller critical value of μ. This phenomenon is called hysteresis. If present in the Earth's climate system, it implies that once the upward jump occurs, it may be very difficult to achieve the reverse jump back to the original climate state. Figure 1b illustrates the situation of a cusp bifurcation, where the two saddlenode bifurcations in Fig. 1a have coalesced, for example, as another parameter of the system is varied. In this situation a small change in μ will also cause a large change in the system F, although it will be smooth and reversible. In Fig. 1c, even though the saddlenode bifurcations may be outside the parameter interval of interest, abrupt large transitions in the system can result from a small noise or perturbationinduced change to the system even when the parameter value remains constant. It is the presence of saddlenode bifurcations in a mathematical model, even if not occurring precisely at the system's current parameter value, that is the root cause of all of the behaviours shown in Fig. 1.
For a tipping point to be present, the underlying mathematical model will be characterized by nonlinearity, generally in the form of a positive feedback that accelerates change once change has begun. For the Arctic, one of the primary positive feedbacks is the surface albedo. When the Arctic Ocean is frozen, the surface reflects a significant portion of the insolation back into space, but open water absorbs much more heat from the Sun. Timing of the melt in the spring has significant impact (Zheng et al., 2021). An earlier melt means considerably more heat is absorbed by open water, raising the water temperature and delaying freezeup in the autumn. The freezeup date for the Beaufort, Chukchi, Laptev, and Kara seas, for example, has been getting later by 6–11 d per decade since 1979 (Stroeve et al., 2014). September sea ice extent has been decreasing at an accelerating rate. The linear trend from 1979 to 2001 is −7 % per decade, but including data up to 2013, the linear trend is −14 % per decade (Stroeve et al., 2014). Thus the observational evidence indicates that the processes behind this phenomenon are not linear at all but nonlinear.
Past studies on general circulation models (GCMs) have given mixed results regarding the presence of multiple stable states for ice conditions in the Arctic. Some indicate that there appears to be a continuous transition from perennial ice cover to annually icefree that is reversible (Schröder and Connolley, 2007; Tietsche et al., 2011; Armour et al., 2011). Other studies have shown evidence for nonlinear behaviour in sea ice loss, especially in the transition from seasonally icefree to annually icefree (Winton, 2006, 2008; Ridley et al., 2008). On the other hand, smaller conceptual models generally show bistability and abrupt transitions in sea ice cover (Thorndike, 1992; MüllerStoffels and Wackerbauer, 2011; Eisenman and Wettlaufer, 2009; Björk and Söderkvist, 2002; Abbot et al., 2011; Merryfield et al., 2008; Flato and Brown, 1996). The most common result from all these models seems to be that sea ice will likely transition from perennial to seasonally icefree in a continuous, reversible manner, but significant warming beyond that point will likely cause an abrupt change to annually icefree (Bathiany et al., 2016). See the introduction in Eisenman (2012). The model we present here is an annually averaged model with no seasonal component. It is not a model of sea ice in particular but rather a column model of the atmosphere that incorporates a nonlinear albedo response to surface temperature. Bistability in our model with both warm and cold solutions corresponds to annually averaged icecovered or icefree situations.
The Arctic climate model presented here is motivated by three observations. First is the observation that the climate changes taking place on the Earth today are most dramatic in the high Arctic. Therefore, it is prudent to put a special focus on understanding Arctic climate change. Second, irreversible change is inevitably the result of nonlinear geophysical processes. So, while this model is kept very simple, it does include significant nonlinear phenomena that can lead to tipping points. Third, the 3D spherical shell of the atmosphere of the Earth is rotationally symmetric about the polar axis if annually and zonally averaged. Due to the rotation of the Earth, Hadley, Ferrel and polar cells form in the global circulation. If perfect rotational symmetry is assumed, the polar axis becomes flowinvariant, and this remains approximately true for the real Earth. Thus, a 1D model restricted to the polar axis can be expected to give useful information about climate in a neighbourhood of the pole. The study of a rotationally symmetric spherical shell model by Lewis and Langford (2008) gives support to this hypothesis. A vertical column of atmosphere at other points on Earth would have a horizontal component of velocity, invalidating the type of analysis used here. Globally averaged climate models do reduce to one (vertical) dimension, but they give little information specific to the Arctic.
The present model builds on the simple energy balance slab model of Dortmans et al. (2019), which was applied to paleoclimate transitions, and the model of Kypke et al. (2020), which was applied to anthropogenic climate change. The primary improvement of the present model is a more physically accurate description of the atmosphere. Instead of using a slab to represent a uniform atmosphere with absorption properties similar to the real atmosphere, here we use the Schwarzschild twostream equations to model absorption in the atmosphere explicitly as a function of altitude (Pierrehumbert, 2010, p. 191).
A bifurcation analysis is performed on the model, tracking the steadystate solutions as carbon dioxide levels increase. The question of reversibility is a question of whether the current cold state simply warms but persists. The disappearance of this cold state through a saddlenode bifurcation would result in an abrupt change in climate that may be practically irreversible. The simpler model of Kypke et al. (2020) showed this behaviour under certain CO_{2} representative concentration pathway scenarios. We seek here to determine whether the present, more accurate model also displays this behaviour.
Section 2 and Appendix A provide a detailed derivation of the model. The model parameter values and calibration of some of them to empirical data are presented in Appendix B. Although much of the detail is relegated to the Appendices, the authors feel this detail constitutes an essential part of the contribution of the paper, providing clarity, justification of choices, and the information necessary for replication. Hence we consider them essential reading. Section 3 presents the results, and the conclusions are in Sect. 4.
The model is developed from first principles and has the following features.

The atmosphere is a 1D column at the North Pole with physical properties that vary with altitude, from the surface to the tropopause.

The incoming solar radiation is annually averaged and undergoes reflection and absorption in the atmosphere as well as at the Earth's surface.

The surface albedo is a nonlinear function of the surface temperature.

A wellmixed surface boundary layer is included.

The Earth emits longwave radiation as a black body.

The atmosphere is considered to be a grey gas.

The Schwarzschild twostream equations govern the absorption and emission of both upward and downwarddirected longwave radiation in the atmosphere.

The atmospheric absorption of longwave radiation is due to three factors: water vapour, CO_{2} concentration, and clouds.

Water vapour concentration is governed by the nonlinear Clausius–Clapeyron equation.

Transfer of latent and sensible heat from the surface to the atmosphere is modelled.

Both ocean and atmospheric meridional heat transports to the Arctic are dictated by empirical values.

In the Arctic, there is a slow downward movement of air in the column corresponding to the polar circulation cell near the pole (Lewis and Langford, 2008; Langford and Lewis, 2009; Lutgens and Tarbuck, 2019). This is achieved via mass transport of air into the column in its upper portion and out of the column near the bottom.

The radiation absorption coefficients are calibrated by fitting the model to global average data.

The functional forms of the mass transport and atmospheric heat transport are used to calibrate the model to an empirical Arctic temperature profile.
The annually and zonally averaged Earth atmosphere is rotationally symmetric around the polar axis, which is invariant under the flow. Therefore, if one considers a column of the atmosphere near the North Pole, it is reasonably approximated by a 1D model with altitudevarying quantities. This approximation becomes exact in the limit as the diameter of the column shrinks to zero. Alternatively, one can view the model as a meridional and zonal average over a cylinder centred at the North Pole. Further, although the Arctic Ocean is not zonally symmetric, in the above view, the contribution of ocean heat transport can be reasonably captured as a scalar quantity. Thus our model is more precisely a model of the North Pole rather than the Arctic. Nonetheless, we do use some empirical data for the region north of 70^{∘} to calibrate the model for two reasons: (1) data further north are not readily available, and (2) the data we use are not likely to alter too much if they were measured closer to the pole. The values for atmospheric heat transport and ocean heat transport are two that may change significantly as one moves north from 70^{∘}, and we therefore analyse the behaviour of the model over a wide range for these parameter values.
The model domain is a vertical cylinder of crosssectional area A (m^{2}) and circumference C_{b} (m). The atmosphere is assumed to be uniform in the crosssectional direction, so that the model's dependent variables can be interpreted as crosssectional averages that vary with the 1D vertical coordinate $z\in [\mathrm{0},{z}_{\mathrm{T}}]$, where z_{T} (m) is the height of the tropopause. This domain is divided into a surface boundary layer with height z_{B} (m), z_{B}≪z_{T}, and the troposphere proper, $z\in [{z}_{\mathrm{B}},{z}_{\mathrm{T}}]$. The model consists of a set of initial value problems (IVPs), with an independent spatial variable z on [0,z_{B}], which can be solved analytically, and a twopoint boundary value problem (BVP) on [z_{B},z_{T}] that depends on the solutions to the IVPs. The model has equations governing the vertical wind speed, w (m s^{−1}), the air density, ρ (kg m^{−3}), the upward and downward longwave radiation, I_{+} and I_{−} (W m^{−2}), the downward shortwave radiation, I_{S} (W m^{−2}), the latent and sensible heat transport, F_{C} (W m^{−2}), and the temperature T (K). Any reflection of shortwave radiation from either the surface or the atmosphere is ignored and is simply considered to have left the system. The model is depicted in Fig. 2 and is derived and explained in detail in the following subsections. Many of the details of the model, including its nondimensionalization, the vanishing conduction limit, and modelling choices used for various functional forms, are in Appendix A. Calibration of the model parameters to empirical data is described in detail in Appendix B, and the reader is referred to Tables B1 and B2 of that appendix for the values of the parameters.
2.1 Mass, momentum, and energy balance
The model equations in the troposphere are developed from the fundamental transport theorem in one spatial dimension:
where f is the density of some “property”, χ is the flux of that property, and S is a source/sink term. The timederivative term will be taken as zero since only the steadystate solution is considered. The properties subject to this equation are mass, momentum, and energy. To model the Arctic, the cylinder is centred at the North Pole, and, since the atmospheric polar cell has slow downward movement near the pole, it is assumed that w<0.
2.1.1 Mass
If the property f in Eq. (1) is the mass density, ρ (kg m^{−3}), then the flux is χ=ρw. There is mass flux across the vertical boundary of the cylinder, M_{b}(z) ($\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$), which is assumed to be immediately spread out evenly across the layer, and hence the mass flux across the vertical boundary in the model is really a mass source term in the interior giving
Thus, at steady state the mass balance equation is
The mass flux through the vertical boundary into the column is written as
where M_{tot} ($\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$) is a nonnegative constant and $\mathit{\varphi}\left(x\right):[\mathrm{0},\mathrm{1}]\mapsto \mathbb{R}$ is a dimensionless function that represents the portion of inward mass flux across the vertical boundary of the column at the given altitude. Positive ϕ indicates inward flow. The ratio of the crosssectional area of the column to the area of its side, $A/\left({C}_{\mathrm{b}}\right({z}_{\mathrm{T}}{z}_{\mathrm{B}}\left)\right)$, in the definition of M_{b}(z) is included as a useful simplifying convenience. The (positive inward) mass fluxes across the bottom and top boundaries of the column are given by M_{tot}Φ_{B} and M_{tot}Φ_{T}, respectively, where Φ_{B} and Φ_{T} are dimensionless constants in the interval $[\mathrm{1},\mathrm{1}]$. These quantities must satisfy
The first of these conditions dictates that there is no net mass entering the system, and the second is a normalization condition so that M_{tot} represents both the total mass entering the column per unit crosssectional area and the magnitude of the total mass leaving the column per unit crosssectional area. In order to obtain downward vertical flow throughout the cylinder, it shall be assumed that Φ_{B}<0, Φ_{T}>0, ϕ≤0 in the lower part of the cylinder and ϕ≥0 in the upper part. With these definitions, the mass balance equation becomes
2.1.2 Momentum
Now take the property f in Eq. (1) to be the momentum density, ρw. The vertical flux is χ=ρw^{2}, and the source term S has two components, one due to contact forces (stress) (van Groesen and Molenaar, 2017, p. 56) and one due to internal body forces (gravity):
where P (N m^{−2}) is the pressure and g (m s^{−2}) is the gravitational acceleration. It is assumed that mass entering the cylinder from the vertical boundary has no vertical momentum. Thus the momentum balance equation at steady state is
(In the case of no flow (w=0), the above would read $\mathrm{d}P/\mathrm{d}z=\mathit{\rho}g$, which is the hydrostatic equation.)
2.1.3 Energy
Finally, consider the case where the property in Eq. (1) is the total energy density given by
which corresponds to the sum of kinetic energy, gravitational potential energy, and internal heat energy densities. Here, c_{v} ($\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$) is the specific heat capacity of the air. The flux has two components, one due to advection and one due to conduction:
where k ($\mathrm{W}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$) is the thermal conductivity. The source/sink S has eight terms, one due to work done by contact forces, two due to mass entering or leaving across the vertical boundary (one of these accounts for gravitational potential energy and the other internal heat energy; there is no addition to kinetic energy since the mass appearing has no velocity), three due to radiation (shortwave downward and longwave upward and downward), one due to latent and sensible heat transport, and one due to atmospheric heat transport. It is important to distinguish the difference between the mass transport across the boundary and the atmospheric heat transport. It is assumed that the mass moving across the vertical boundary is at the same temperature as the mass inside at each altitude. Since mass transfer across the vertical boundary is inward in the upper portion, where the temperature is cooler, and outward in the bottom portion, mass transport results in a net transport of heat out through the vertical boundary, but this will be small since the mass flux, M_{b}(z), is small. The reason M_{b}(z) is small is that it generates the average slow movement of air downward near the North Pole (about 1 mm s^{−1}), due to the circulation of the polar cell. This slow averaged circulation of air does not account for the atmospheric heat transport. The main transport of heat in the atmosphere is via turbulent mixing captured in our model by a source term, F_{A}(z) (W m^{−3}), whose functional form is discussed in Sect. A3. Thus, S is given by
The governing equations for the longwave radiation intensities are the twostream Schwarzschild equations and for the shortwave radiation a standard absorption equation:
where k_{S} (m^{2} kg^{−1}) is the shortwave absorption coefficient, σ ($\mathrm{W}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{4}}$) is the Stefan–Boltzmann constant, and κ (m^{−1}) is the longwave absorption coefficient with terms corresponding to absorption by clouds, carbon dioxide, and water vapour:
Here, k_{Cl} (m^{−1}), k_{C} (m^{2} kg^{−1}), and k_{W} (s^{2} kg^{−1}) are absorption coefficients that will be calibrated, μ (ppm) is the CO_{2} concentration expressed as the ratio of moles of CO_{2} to moles of dry air, ${M}_{{\mathrm{CO}}_{\mathrm{2}}}$ and M_{A} (kg mol^{−1}) are the molar masses of CO_{2} and dry air, respectively, $\mathit{\delta}\left(\right(z{z}_{\mathrm{B}})/({z}_{\mathrm{T}}{z}_{\mathrm{B}}\left)\right)$ is the relative humidity at altitude z, and P^{sat}(T) (N m^{−1}) is the saturated water vapour partial pressure at temperature T. The dependence of this last quantity on T is given by the Clausius–Clapeyron equation
where P^{sat}(T_{R}) is the pressure at a reference temperature T_{R} (which we take to be 273.15 K), L_{v} (m^{2} s^{−2}) is the latent heat of vaporization for water, and ${R}_{\mathrm{W}}=R/{M}_{\mathrm{W}}$ ($\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}$) is the gas constant for water, R ($\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{mol}}^{\mathrm{1}}$) is the universal gas constant, and M_{W} (kg mol^{−1}) is the molar mass of water. The corresponding density is ${\mathit{\rho}}_{\mathrm{W}}^{\mathrm{sat}}\left(T\right)={P}^{\mathrm{sat}}\left(T\right)/\left({R}_{\mathrm{W}}T\right)$ by the ideal gas law. The vertical heat transport (latent and sensible heat) is assumed to be governed by a simple exponential decay:
where b (m^{−1}) is a suitable decay constant. Substituting these expressions for the total energy density and its flux and sources into Eq. (1) and combining the expressions from Eqs. (5) to (10), the energy balance equation at steady state for $z\in [{z}_{\mathrm{B}},{z}_{\mathrm{T}}]$ is given by
The nonlinear effects of both water vapour and carbon dioxide concentration on longwave radiation absorption in the atmosphere are contained within the factor κ(ρ,T), defined by Eq. (8). In earlier work (Dortmans et al., 2019; Kypke et al., 2020), these two effects were studied separately before combining them. It was shown there that, if the atmosphere becomes warmer, then the concentration of water vapour increases due to the Clausius–Clapeyron relation, and this accelerates the greenhouse warming of the atmosphere well beyond that due to carbon dioxide alone. This is an important positive feedback in the model.
In order to complete the system, a constitutive relation between the density ρ and the pressure P is needed, for which we use the ideal gas law,
where ${R}_{\mathrm{A}}=R/{M}_{\mathrm{A}}$ ($\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$) is the gas constant for air.
The mass, momentum, and energy balance equations, Eqs. (3), (4), and (11), along with the Schwarzschild equations (Eqs. 5 and 6) and the equations governing shortwave absorption (Eq. 7) and sensible and latent heat transport (Eq. 10) are the differential equations for the BVP for $z\in [{z}_{\mathrm{B}},{z}_{\mathrm{T}}]$, with dependent variables w, ρ, I_{+}, I_{−}, I_{S}, F_{C}, T, and $\mathrm{d}T/\mathrm{d}z$. Equations (8), (9), and (12) define certain quantities in these differential equations in terms of these dependent variables. The forms of the functions F_{A}(z) and ϕ(z) are prescribed; the process of choosing these functions is described in detail in Appendix A3.2 and A3.3, respectively. Before discussing the boundary conditions for the BVP, it is necessary to consider the surface boundary layer.
2.2 Surface boundary layer
The model includes a boundary layer extending from z=0 to z=z_{B}. It is assumed that this layer is wellmixed so that temperature T_{B}=T(z_{B}), density ρ_{B}=ρ(z_{B}), and relative humidity δ_{B}=δ(z_{B}) in this layer are constant. The temperature of the surface, T_{S}, can in general be larger or smaller than T_{B}.
The primary reason for including a boundary layer is a numerical one. As shown in Appendix A2, the model is numerically stiff due to the thermal conductivity of air being very small. To remove the stiffness, a limit to vanishing conduction is taken, and this results in an algebraic expression for the temperature gradient that includes the vertical wind speed as a factor in the denominator. As the vertical wind speed must be zero at the Earth's surface, there is a singularity in the temperature gradient there. The introduction of the surface boundary layer avoids this singularity.
The total mass crossing from the atmosphere into the boundary layer per unit time is M_{tot}Φ_{B}A. This quantity is negative, since Φ_{B}<0, indicating flow out of the atmosphere and into the boundary layer. This mass exits through the vertical boundary of the layer with an assumed constant mass flux K, at each z, and conservation of mass dictates
(K<0 indicates the flux is outward.) This exiting mass carries gravitational potential energy. The change in potential energy in a slab of height Δz at height z in the boundary layer is C_{b}ΔzKgz, so that the total change in potential energy over the boundary layer is
Consistent with the modelling assumption that mass flux across the vertical boundary conveys no momentum or kinetic energy to the system, the loss of mass out of the vertical boundary of the boundary layer also has no effect on the momentum or kinetic energy. Further, since the temperature in the boundary layer, T_{B}, is assumed to be equal to the temperature of the atmosphere at z=z_{B}, it follows that there is also no net energy change in the boundary layer due to advection of internal energy – the internal energy entering via advection at the top of the layer is equal to the internal energy leaving the layer through the vertical boundary.
Consider now the energy balance at the Earth's surface. There is energy transport from the surface to the boundary layer in the form of sensible and latent heat, which is modelled, as per Pierrehumbert (2010, p. 396–398), as
where C_{D} is a dimensionless drag coefficient, U (m s^{−1}) is the horizontal wind speed, and P^{sat}(T) is given by Eq. (9). Along with this there is energy input to the surface from the Sun, I_{S}(0), some of which is reflected by the surface albedo, longwave radiation both inward, I_{−}(0), and outward, I_{+}(0), and ocean heat transport, F_{O} (W m^{−2}). Therefore, the energy balance at the surface is
where
is the surface albedo, here modelled as a sigmoid function increasing from α_{c} at cold temperatures to α_{w} at warm temperatures, with the midway point being at the reference temperature T_{R} (freezing point) and with a steepness of transition determined by the dimensionless constant ω.
Now consider the energy balance for the combined surface and boundary layer (one could alternatively consider just the boundary layer without the surface, but the chosen formulation results in a slightly smaller equation). Input energy to this combined surface and boundary layer includes ocean heat transport and shortwave and longwave radiation entering at z_{B}. Output energy includes upward longwave radiation at z_{B}, the shortwave radiation reflected from the surface, and the latent and sensible heat F_{C} at z_{B}. Further, there are kinetic and gravitational potential energy fluxes and heat conduction in/out of the layer through its top at z_{B}, and there is gravitational potential energy loss through the vertical boundary given by Eq. (13). Therefore, the energy density balance for the combined surface and boundary layer is
Since temperature, pressure, and relative humidity are constant in the boundary layer, the radiation equations may be solved analytically inside the layer in order to relate the radiation terms at z=0 to those at z=z_{B}. The simple ordinary differential equation (ODE) for F_{C} is also easily solved in the boundary layer. The initial (independent spatial variable z=0) condition for the upward longwave radiation, I_{+}(0), is that it is equal to the black body radiation of the surface, $\mathit{\sigma}{T}_{\mathrm{S}}^{\mathrm{4}}$. The initial condition for that latent heat, F_{C}(0), is given by Eq. (14). Initial conditions for I_{−} and I_{S} are not necessary since only a relation between the values of these functions at 0 in terms of their value at z_{B} is required. The IVPs for I_{+} and F_{C} and the ODEs for I_{−} and I_{S} in the boundary layer are
and their solutions, via standard means, give
Equations (18) and (19) provide two boundary conditions for the BVP on the troposphere. The energy balance equations (Eqs. 15 and 17) along with Eqs. (16), (20), and (21) provide two further boundary conditions.
2.3 Boundary conditions for the BVP
There are eight unknown dependent variables: w, ρ, I_{+}, I_{−}, I_{S}, F_{C}, T, and $\mathrm{d}T/\mathrm{d}z$; in addition, the surface temperature, T_{S}, is an unknown constant (independent of z) that is determined through Eq. (14), while the pressure, P, can be written in terms of the others via Eq. (12). The boundary conditions for the system on the interval [z_{B},z_{T}] are wind speed at z_{B} given by the requirement that the advected mass Aw(z_{B})ρ_{B} equals the mass flux M_{tot}Φ_{B}A, pressure at the surface equal to the standard pressure, P_{0}, upward longwave radiation at z_{B} given by Eq. (18), vertical heat transport at z_{B} given by Eq. (19), the energy balance equations (Eqs. 15 and 17) with expressions from Eqs. (20) and (21) substituted in, no downward longwave radiation at z_{T}, shortwave radiation at z_{T} equal to the insolation, Q, less what is reflected by the clouds, Q_{R}, and a local critical point for T at z_{T}, which, respectively, correspond to the following equations:
where F_{C0} is given by Eq. (14). The last three terms of Eq. (27) may be simplified using Eq. (22) so that they read
There are nine boundary conditions, but there are only eight dependent variables for which we have differential equations. The discrepancy is explained by the presence of T_{S}, which is an additional scalar unknown. The nine boundary conditions determine eight conditions for the differential equations as well as the value for T_{S}. One way of treating this is simply to extend the system of differential equations to include the equation
In addition, as described in Appendix A, to avoid numerical stiffness, we take the limit as the heat conduction of air, k, tends to zero. This effectively reduces the size of the model by one dimension.
The model is nondimensionalized and put in standard form as detailed in Appendix A. The parameter values and their calibration to empirical data are provided in Appendix B.
For the Arctic parameter values given in Appendix B
and for a given CO_{2} concentration, μ, the model can
be solved numerically.
We used MATLAB's builtin BVP solver “bvp5c
” to solve individual instances of the (nondimensionalized) model and AUTO for continuation calculations. The results of the model for μ=390 ppm are shown in Fig. 3.
The altitude dependence of the mass flux ϕ and the atmospheric meridional heat
transport, F_{A}, were calibrated to an empirical Arctic temperature profile
from Cronin and Jansen (2016) as detailed in Appendix B.
From this figure we see that the model fits the temperature profile
very well, with some discrepancy near the surface.
The overall atmospheric heat transport, Fig. 3c,
indicates that, for our model, the upper half of the troposphere receives the most input of heat, while the bottom fifth actually has a net outward heat transport.
This removal of heat near the bottom is likely the cause of the discrepancy between
our model values and the Cronin and Jansen data.
The negative values of F_{A} near the surface are due to our modelling choice for F_{A}.
It is possible that alternative modelling formulations for F_{A} could yield
a better fit to the Cronin and Jansen data; however, the forms that we tried (including the ones reported in Appendix B and a few others) did not improve upon the one used here.
Nonetheless, given the relative simplicity of our model, we feel the fact that it can fit
the data as well as it does is remarkable.
For μ=390, the surface temperature from the solved model is −19.7 ^{∘}C. Starting from this solution, we numerically continued the solution to other values of μ, resulting in the Sshaped curve in Fig. 4a. Between approximately μ=464 and μ=859, there are three solutions. The lower and upper solutions are stable cold and warm solutions, while the middle branch is unstable. Currently the Arctic is on the cold branch of this curve. The model predicts that, as CO_{2} levels rise, the equilibrium surface temperature in the Arctic will increase gradually at first. However, when μ exceeds 859 ppm, the model displays a saddlenode bifurcation where the stable cold solution is annihilated together with the unstable solution. At this point the climate would rapidly approach the warm stable solution where surface temperatures are significantly higher. The temperatures on the warm branch may seem unreasonably high. Although our model includes specific details of some of the natural phenomena governing the Arctic climate behaviour and is calibrated with real data, it is primarily a qualitative rather than quantitative model. It is the fact that the model predicts the qualitative feature of a saddlenode bifurcation that is the important result, not the precise temperature of the model's warm solution. Part of the reason the model has relatively high temperatures on the warm branch is because it has constant values for F_{O} and ${F}_{\mathrm{A}}^{\mathrm{tot}}$. It is likely that, should the Arctic's annually averaged temperature rise dramatically, both F_{O} and ${F}_{\mathrm{A}}^{\mathrm{tot}}$ would be affected in a downward direction, which would reduce the warm equilibrium temperatures to some degree. Another important emergent feature of the model is that there is bistability for μ in the range [464,859]. Thus, even though the Arctic may in the future be on the cold branch in this range, it is possible that a strong enough climate disturbance could push the climate out of the basin of attraction of the cold solution and into the basin of attraction of the warm solution. The necessary strength of such a disturbance decreases as one moves closer to the upper end of this range. Our model does not incorporate the seasonal variation of solar input to the Arctic but rather uses an annually averaged value. Thus the equilibria in our model are annual averages. The seasonal variation of insolation would effectively result in oscillations around the annual average. These oscillations themselves may make a significant contribution to the “disturbance” needed to push the system out of the basin of attraction of the cold branch.
The International Panel on Climate Change (IPCC) has published various CO_{2} emission scenarios for the future based on possible levels of global action to suppress such emissions (Intergovernmental Panel on Climate Change, 2013, Box TS.6; van Vuuren et al., 2011). These scenarios are called Representative Concentration Pathways (RCPs) and are numbered based on the radiative forcing in the year 2100 due to anthropogenic emissions compared to the year 1750. The original four published RCPs are RCP2.6, RCP4.5, RCP6.0, and RCP8.5, representing strong mitigation (2.6) through “ignore the problem” (8.5) responses by world governments. Each RCP indicates likely CO_{2} concentration levels in the atmosphere out to the year 2100. From 2100 to 2200, the scenarios assume a “constant composition commitment”, which essentially freezes emission levels and eventually leads to a constant CO_{2} level in the atmosphere for all but RCP8.5. These RCPs are plotted in Fig. 4b. The CO_{2} levels for the four different pathways at the year 2200 are continued as dashed lines into Fig. 4a. It is clear from the figure that RCP8.5 leads to CO_{2} concentrations that far surpass the saddlenode bifurcation, whereas the other three RCPs do not. This result is in agreement with our simpler model (Kypke et al., 2020). Both RCP6.0 and RCP4.5 end at levels within the bistable range, and indeed all the RCPs except RCP2.6 are within that range from about the year 2050 onward. The black dasheddotted line extending from the bifurcation in the upper panel to RCP8.5 in the lower panel illustrates that RCP8.5 reaches the bifurcation near the year 2092.
The curve of equilibria in Fig. 4a displays hysteresis: CO_{2} levels rising past 859 ppm will cause a jump from the cold equilibrium state to the warm state, but a return to the cold state will not happen until CO_{2} levels are brought below 464 ppm, where the saddlenode bifurcation of the warm equilibrium is located (left bend of the S curve). If CO_{2} levels follow a trajectory similar to RCP8.5, Arctic climate may change drastically in less than 100 years, but a return to the current cold state may be essentially impossible for thousands of years afterward, assuming humankind can develop and implement the required technology to reduce atmospheric carbon levels sufficiently.
Total atmospheric heat transport, ${F}_{\mathrm{A}}^{\mathrm{tot}}$, and ocean heat transport, F_{O}, are inputs to the model that are empirically based and that are not altered by the state of the model itself. We investigated the persistence of the saddlenode bifurcation and hysteresis in the presence of the uncertainty of these two values. Figure 5 shows two 2D bifurcation diagrams plotting the locations of the saddlenode bifurcations in the model (the left and right bends of the S curve in Fig. 4). Figure 5a shows the bifurcations as both ${F}_{\mathrm{A}}^{\mathrm{tot}}$ and μ vary. Figure 5b shows a similar diagram where F_{O} replaces ${F}_{\mathrm{A}}^{\mathrm{tot}}$. The right curve in each of the panels is a curve of saddlenode bifurcations corresponding to the right bend in the S curve in Fig. 4; they are transitions from a cold state to a warm state. Similarly, the left curve in each of the panels is a curve of saddlenode bifurcations that correspond to the left bend in the S curve and are the transitions from a warm state to a cold state. The shaded area between the two curves is the region of bistability. The curves meet in a cusp bifurcation, but it is important to note that this only occurs when CO_{2} levels are mathematically negative. Thus the two saddlenode bifurcations and the hysteresis phenomenon will be present for all physically possible CO_{2} levels and for all reasonable values of ${F}_{\mathrm{A}}^{\mathrm{tot}}$ and F_{O}. These bifurcation curves also make it evident that a transition from a cold state to a warm state occurs if either ${F}_{\mathrm{A}}^{\mathrm{tot}}$ or F_{O} are increased, even if CO_{2} levels are constant. If ${F}_{\mathrm{A}}^{\mathrm{tot}}$ and F_{O} increase along with carbon dioxide levels in the near future, which seems reasonably likely, then the saddlenode bifurcation will occur at lower levels of CO_{2}. For example, a 10 % increase in ${F}_{\mathrm{A}}^{\mathrm{tot}}$ to 110 W m^{−2} from 100 would cause the saddlenode bifurcation to occur at about μ=754 ppm, which would mean that both RCP8.5 and RCP6.0 would pass through the saddlenode bifurcation.
Although the model presented here is clearly a simplification of the climate, made possible by the near invariance of the vertical flow on the polar axis, we believe it captures some of the most important aspects relevant for Arctic climate change. The model predicts that if humanity keeps carbon emission levels close to RCP6.0 or lower, then the Arctic will not likely undergo a sudden dramatic rise in annual average temperature. However, if carbon emissions are much worse than RCP6.0, such a change is likely, and the cause is a saddlenode bifurcation of the stable cold equilibrium. Such a change would clearly have catastrophic effects on the Arctic environment, leading to massive global effects. These results are in agreement with those of Årthun et al. (2021), who, from a study of various CMIP6 models of Arctic climate, predict that under a low emissions scenario sea ice loss will be seasonal but that for a high emissions scenario it will be yearround for all areas of the Arctic. Further, the hysteresis displayed by the model indicates that a change of this nature may be practically irreversible. Although some comfort may be taken that only the worst of the four carbon pathway scenarios ends in such a catastrophic change, the model shows that both RCP6.0 and RCP4.5 lie in the region of bistability from the year 2070 onward. In a bistable situation, external disturbances could cause the system to jump into the basin of attraction of the warm equilibrium, effectively bringing about the catastrophic change prior to the bifurcation. The model is too simple to allow for any reasonable measure of the likelihood of such an occurrence, but the important thing is that the model exhibits bistability in the parameter region where the system is likely to reside within 50 years' time. This bistability has been shown to persist regardless of the values of the two biggest uncertainties in the model, the atmospheric and ocean heat transport to the Arctic from the midlatitudes.
Our model addresses the equilibrium state only and represents the Arctic temperature as an annual average. The real Arctic climate undergoes massive seasonal changes, which effectively means that the system is actually oscillating around the equilibrium temperatures of our model. Such temperature oscillations may be sufficient to effectively push a system located on a cold solution in the bistable regime to “above” the unstable solution and so into the basin of attraction of the warm solution.
Seasonal variations in Arctic temperatures and sea ice are studied in many works, including Eisenman and Wettlaufer (2009), who argue that a tipping point to a yearround icefree Arctic is not likely to occur while the Arctic is icecovered for a significant portion of the year, but once the Arctic is seasonally icefree for a sufficient number of months in a year, such a tipping point becomes more likely. A possible enhancement to our model would be to include seasonal solar variation and an ice model as in Eisenman and Wettlaufer (2009).
This Appendix provides details regarding the model. The model is written in a nondimensional form in Sect. A1, defining nondimensional variables and parameters. This system is then transformed in Sect. A2 to a standard form of nine firstorder ordinary differential equations and corresponding boundary conditions. This standard form makes it evident that the system is numerically stiff due to the fact that the thermal conductivity, k, of air is very small. To remove this stiffness, the limit as k→0 is applied, which reduces the system by one dimension. Section A3 discusses our choice of the functional forms for the dependence of relative humidity δ(z), atmospheric meridional heat transport F_{A}(z), and mass flux M_{b}(z) on altitude.
A1 Nondimensionalization
Define the nondimensional variables
where P_{0} (N m^{−2}) is the standard atmospheric pressure, T_{R} is the freezing point of water in Kelvin, ${\mathit{\rho}}_{\mathrm{0}}={P}_{\mathrm{0}}/\left({R}_{\mathrm{A}}{T}_{\mathrm{R}}\right)$ is the density at standard pressure and freezing temperature, and $\mathit{\sigma}{T}_{\mathrm{R}}^{\mathrm{4}}/\left({c}_{v}{\mathit{\rho}}_{\mathrm{0}}{T}_{\mathrm{R}}\right)$ represents the vertical velocity required to move a parcel of air with standard density at freezing temperature so that the power transferred is equal to the radiative power for a black body at the same temperature. This comes to about 1 mm s^{−1}.
Applying the change in variables to the troposphere BVP, Eqs. (3)–(7), (10), (11), and (31), we get
for $\widehat{z}\in [\mathrm{0},\mathrm{1}]$. The boundary conditions, Eqs. (22)–(30), become
where
are nondimensionalized functions describing the sensible/latent heat flux from the surface, the absorption of longwave radiation due to greenhouse gases, the atmospheric heat transport, and the surface albedo, and α_{w}, α_{c}, ω, Φ_{B}, and
are nondimensional constants. Values of the physical parameters are given in the tables in Appendix B. It turns out that all of the above nondimensional constants are close to order 1 (range 0.04 to 26) except $\mathit{\zeta}=\mathrm{5.6}\times {\mathrm{10}}^{\mathrm{3}}$, which is the relative boundary layer thickness, $\mathit{\u03f5}=\mathrm{2.35}\times {\mathrm{10}}^{\mathrm{6}}$, which is the nondimensional conductance, and $H=\mathrm{7.96}\times {\mathrm{10}}^{\mathrm{12}}$, which is a factor in the kinetic energy term. As shown in the next section, the fact that ϵ is small causes system stiffness.
A2 Standard form and vanishing conduction limit
The BVP given by Eqs. (A2)–(A19) contains derivatives of products of some of the variables. It can be put in standard ODE form via algebraic manipulations. First we expand the derivative on the lefthand side of Eq. (A3) and use Eq. (A2) to simplify it:
Thus, Eq. (A3) may be replaced with
where we have used Eq. (A8) to replace the derivative of y_{7}. We can also expand the derivative in Eq. (A2) to get the equivalent equation
Equations (A25) and (A26) are a linear system in $\mathrm{d}{y}_{\mathrm{1}}/\mathrm{d}\widehat{z}$ and $\mathrm{d}{y}_{\mathrm{2}}/\mathrm{d}\widehat{z}$, namely
Solving this yields
Now expanding the derivatives on the righthand side of Eq. (A7) and simplifying, we obtain
In summary, the BVP for the troposphere in standard form is given by
with boundary conditions
As mentioned above, of the nondimensional constants, all are close to order 1 except $\mathit{\zeta}=\mathrm{5.6}\times {\mathrm{10}}^{\mathrm{3}}$, $H=\mathrm{7.96}\times {\mathrm{10}}^{\mathrm{12}}$, and $\mathit{\u03f5}=\mathrm{2.3}\times {\mathrm{10}}^{\mathrm{6}}$ (with z_{T}=9000). The constant ζ only appears in the boundary conditions. The constant H occurs in the system only in summations with other nonderivative terms; hence the fact that it is small only means that those terms contribute little. It does not cause stiffness. The relatively small value of the constant ϵ, however, does cause stiffness in Eq. (A37) of the system. Because ϵ is small, the variable y_{8}, which is the dimensionless rate of vertical temperature change $\mathrm{d}T/\mathrm{d}z$, will approach the y_{8}nullcline very rapidly (which means at a very short z distance from either boundary). To simplify numerical computation, we take the limit as ϵ goes to zero, which is equivalent to saying that conduction is negligible. In this limit, Eq. (A37) is an algebraic expression from which we can isolate y_{8}:
With this expression for y_{8}, the system is reduced by eliminating Eq. (A37) and boundary condition Eq. (A47). Since y_{1}, the dimensionless vertical velocity, is a factor in the denominator, it is necessary that y_{1} be nonzero throughout the atmosphere in order not to introduce a singularity. For this reason we select Φ_{T}>0, Φ_{B}<0, and ϕ such that it is nonnegative in the upper atmosphere and nonpositive in the lower atmosphere.
A3 Modelling choices for functional forms
In this section we describe the various functional forms that we used for relative humidity, atmospheric heat transport, and mass flux. For the latter two, several different forms were tried, and these are detailed below. Calibration to empirical data, described in Appendix B, was used to select specific functional forms for the heat transport and mass flux.
A3.1 Relative humidity
The relative humidity is modelled as a linear function decreasing with altitude from a higher surface value, δ_{B}, to a lower value at the tropopause, δ_{T}. Specifically,
A3.2 Atmospheric heat transport
Atmospheric heat transport is primarily due to largescale turbulent mixing of the column with its environment. This mixing is not modelled explicitly, but, instead, it is incorporated into the model via the function F_{A}(z), which represents the thermal energy supplied to the column by turbulent mixing. The integral of F_{A}(z) over the atmosphere thickness represents the total atmospheric heat transport in/out of the system. So, given a set amount of such energy, ${F}_{\mathrm{A}}^{\mathrm{tot}}$, we have
This provides one restriction on the function F_{A}(z), but its precise form is a modelling choice. However, the exact form must be chosen with care, because with certain choices of F_{A}(z), the boundary conditions can only be satisfied with unrealistic solutions. In particular, if the values of F_{A}(z) near the tropopause (z near z_{T}, that is, $\widehat{z}$ near 1) are too small, then the temperature drops precipitously toward absolute zero; if they are too big, the temperature turns around and climbs rapidly. Thus, in order to automate an appropriate choice of F_{A}(z), we have proceeded as follows.

Choose ${\widehat{F}}_{\mathrm{A}}\left(\mathrm{1}\right)$ such that the temperature gradient at the tropopause is zero, that is, reimpose boundary condition Eq. (A47). Thus use Eq. (A48) with y_{8}(1)=0 to solve for ${\widehat{F}}_{\mathrm{A}}\left(\mathrm{1}\right)$. Let ${\widehat{F}}_{\mathrm{A}\mathrm{1}}$ denote this value of ${\widehat{F}}_{\mathrm{A}}\left(\mathrm{1}\right)$, and let F_{A1} denote the corresponding dimensional value of F_{A}(z_{T}), that is, ${F}_{\mathrm{A}\mathrm{1}}={\widehat{F}}_{\mathrm{A}\mathrm{1}}\mathit{\sigma}{T}_{\mathrm{R}}^{\mathrm{4}}/({z}_{\mathrm{T}}{z}_{\mathrm{B}})$.

Assume that F_{A}(z) is of the form
$$\begin{array}{}\text{(A51)}& {F}_{\mathrm{A}}\left(z\right)={F}_{\mathrm{Ab}}\left(z\right)+{\displaystyle \frac{{F}_{\mathrm{A}}^{\mathrm{tot}}}{{z}_{\mathrm{T}}{z}_{\mathrm{B}}}}\mathit{\psi}\left({\displaystyle \frac{z{z}_{\mathrm{B}}}{{z}_{\mathrm{T}}{z}_{\mathrm{B}}}}\right),\end{array}$$where the base function F_{Ab}(z) is a linear function passing through zero at the midpoint of the atmosphere and equal to F_{A1} at z_{T}, that is,
$$\begin{array}{rl}& {F}_{\mathrm{Ab}}\left(z\right)={F}_{\mathrm{A}\mathrm{1}}{\displaystyle \frac{\mathrm{2}z({z}_{\mathrm{T}}+{z}_{\mathrm{B}})}{{z}_{\mathrm{T}}{z}_{\mathrm{B}}}},\\ & \u27fa{\widehat{F}}_{\mathrm{Ab}}\left(\widehat{z}\right)={\widehat{F}}_{\mathrm{A}\mathrm{1}}(\mathrm{2}\widehat{z}\mathrm{1}).\end{array}$$This base portion of F_{A} contributes no net heat to the column; it is simply a factor that essentially moves heat around in the column in order to ensure the temperature gradient at the top is zero.
The remaining portion of F_{A} is the actual atmospheric heat transport entering the column from outside. We assume that $\mathit{\psi}\left(\widehat{z}\right)$ satisfies
$$}\mathit{\psi}\left(\mathrm{1}\right)=\mathrm{0}\phantom{\rule{2em}{0ex}}\text{and}\phantom{\rule{2em}{0ex}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int}}\mathit{\psi}\left(x\right)\mathrm{d}x=\mathrm{1},$$so that the value of F_{A} is not altered at the tropopause and so that ${F}_{\mathrm{A}}^{\mathrm{tot}}$ (W m^{−2}) represents the total energy flux of atmospheric heat transport entering the column. Thus, the nondimensional function is
$${\widehat{F}}_{\mathrm{A}}\left(\widehat{z}\right)={\widehat{F}}_{\mathrm{A}\mathrm{1}}(\mathrm{2}\widehat{z}\mathrm{1})+{\widehat{F}}_{\mathrm{A}}^{\mathrm{tot}}\mathit{\psi}\left(\widehat{z}\right),$$where
$${\widehat{F}}_{\mathrm{A}}^{\mathrm{tot}}={\displaystyle \frac{{F}_{\mathrm{A}}^{\mathrm{tot}}}{\mathit{\sigma}{T}_{\mathrm{R}}^{\mathrm{4}}}}.$$ 
As a numerical issue, since the boundary condition value ${\widehat{F}}_{\mathrm{A}\mathrm{1}}$ is needed in the computation of the vector field and since the MATLAB solver we are using does not have a way of making boundary condition information available to the vector field computation function, we circumvented this issue by adding another variable to the problem y_{10}, with differential equation $\mathrm{d}{y}_{\mathrm{10}}/\mathrm{d}\widehat{z}=\mathrm{0}$ and boundary condition ${y}_{\mathrm{10}}\left(\mathrm{1}\right)={\widehat{F}}_{\mathrm{A}\mathrm{1}}$.

The choice of the functional form of $\mathit{\psi}\left(\widehat{z}\right)$ is somewhat open, and we tested the following two forms:
$$\begin{array}{}\text{(A52)}& {\displaystyle}\mathit{\psi}\left(\widehat{z}\right)& {\displaystyle}={g}_{\mathrm{1}}(\mathrm{1}\widehat{z},{L}_{\mathit{\psi}}),\text{(A53)}& {\displaystyle}\mathit{\psi}\left(\widehat{z}\right)& {\displaystyle}={g}_{\mathrm{2}}(\mathrm{1}\widehat{z},{L}_{\mathit{\psi}}),\end{array}$$where the functions g_{1} and g_{2} are defined as
$$\begin{array}{ll}{\displaystyle}{g}_{\mathrm{1}}(x,L)=& {\displaystyle}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{L\mathit{\pi}}{\mathrm{1}\mathrm{cos}\left(L\mathit{\pi}\right)}}\mathrm{sin}\left(L\mathit{\pi}x\right),\\ \text{(A54)}& {\displaystyle}{g}_{\mathrm{2}}(x,L)=& {\displaystyle}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{\mathrm{2}L\mathit{\pi}}{\mathrm{2}L\mathit{\pi}\mathrm{sin}\left(\mathrm{2}L\mathit{\pi}\right)}}\left(\mathrm{1}\mathrm{cos}\left(\mathrm{2}L\mathit{\pi}x\right)\right),\end{array}$$and where L is a parameter in (0,1] free to be chosen. (The functions g_{1} and g_{2} will also be utilized in the modelling of $\mathit{\varphi}\left(\widehat{z}\right)$.) The primary difference between these two forms is that the first has a nonzero slope at x=0, while the latter has a zero slope there.
A3.3 Mass flux
The function ϕ dictates the mass flux across the vertical boundary (negative outward), and, along with the fluxes Φ_{B} and Φ_{T} across the bottom and top of the column, drives the vertical air movement in the column. The only general restrictions on these fluxes are given by Eq. (2). To model the situation in the Arctic, we want a downward flow of air with a vertical wind speed, w, of the order of 1 mm s^{−1} in the column. A reasonable assumption at the tropopause would be to set w=0; however, since our model has a singularity when w=0, we impose a small wind speed at the tropopause by ensuring Φ_{T} is positive. At the surface boundary layer we impose Φ_{B}≤0. Further, to simplify matters and to ensure a downward flow throughout the column, we assume that
where z_{c} is some point in [0,1]. The following forms were tested for $\mathit{\varphi}\left(\widehat{z}\right)$:
piecewise constant,
piecewise linear,
piecewise sine,
piecewise cosine,
where g_{1} and g_{2} are defined by Eq. (A54). (In the case that z_{c}=0 or z_{c}=1, it is understood that only the nonempty interval for ϕ in the above definitions is used and that it is closed at both ends.)
This Appendix lists the parameter values used in the model and discusses how some of them were calibrated to empirical data. Section B2 gives the calculation of the average annual insolation for the Earth north of 70^{∘} N.
B1 Parameter values and calibration
Values of the model parameters are given in Tables B1 and B2. The parameters in Table B1 are physical constants. The parameters in Table B2 are those that have been assigned from empirical data and knowledge or whose values have resulted from fitting the model to empirical data. The model was applied to the globally averaged situation for the purposes of calibration of some parameters (see below) and then also applied to the Arctic.
^{a} Values used for calibration only. ^{b} For the purposes of calibrating the model to global data, both α_{c} and α_{w} were set to the empirical value of $\mathrm{24}/\mathrm{185}$. When calibrating the model's mass flux and atmospheric heat transport parameters to Arctic data, both α_{c} and α_{w} were set to the empirical value of $\mathrm{2}/\mathrm{3}$.
Here we provide justification and explanation of our choice of parameter values in Table B2. The height of the boundary layer was set to z_{B}=50 m. The model is not very sensitive to this parameter. The height of the tropopause is about 9 km at the poles and 17 km at the Equator, so we used the lower value for the Arctic and a middle value of 14 km for the global average. The globally averaged insolation, the atmospheric solar reflection, and the average surface albedo are all obtained from Wild et al. (2013). For the Arctic, the insolation is the annual average for the region north of 70^{∘} N, the calculation of which is shown in Sect. B2. For the average global situation, there is no ocean or atmospheric heat transport; the corresponding values for the Arctic come from Mayer et al. (2019) and Serreze et al. (2007). The Arctic atmospheric reflection and surface albedo come from Kalnay et al. (1996), National Centers for Environmental Protection/National Weather Service/NOAA/U.S. Department of Commerce (2021), and Mayer et al. (2019). Relative humidity is low at the top of the troposphere, so in both the global and Arctic cases it was set to 10 %. The surface relative humidity was set to 75 % for the global average and 70 % for the Arctic. The surface horizontal wind speed, U, which is a factor in the sensible and latent heat transport from the surface, was set to 10 m s^{−1}. The exact value is not too important since U always appears multiplied by the drag coefficient factor C_{D}, which we calibrate to data below. For the global model it is appropriate to assume that the average vertical wind speed is zero, but since the model requires a nonzero wind speed, we set ${M}_{\mathrm{tot}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, and we set ${\mathrm{\Phi}}_{\mathrm{B}}=\mathrm{1}$ and Φ_{T}=0.2 and assumed $\mathit{\varphi}\left(\widehat{z}\right)$ is given by Eq. (A57), with z_{c}=0 and L_{ϕT}=1. (L_{ϕB} is irrelevant since z_{c}=0.) These settings make the wind speed relatively constant and of the order of 10^{−3} mm s^{−1}, which is far enough away from the singularity to avoid convergence issues but small enough so all of the convectionrelated terms in the model become negligible. Mass flux for the Arctic situation was set with trial and error to $\mathrm{8.0}\times {\mathrm{10}}^{\mathrm{4}}$ $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, which gave vertical wind speeds in the column of the order of 0.5 mm s^{−1}. Since ${F}_{\mathrm{A}}^{\mathrm{tot}}=\mathrm{0}$ for the global case, the form of ψ and therefore also the parameter L_{ψ} are not relevant.
Calibration of the other model parameters was done in two steps. First, the absorption coefficients, k_{S}, k_{C}, k_{W}, and k_{Cl}, the decay for sensible and latent heat transport, b, and the drag coefficient, C_{D} (which is a multiplicative factor of both B_{2} and B_{3}), were calibrated using global average energy fluxes obtained from Wild et al. (2013). In addition, this calibration attempted to match estimates from Schmidt et al. (2010), indicating that 25 % of absorption is due to carbon dioxide, 25 % due to clouds, and 50 % due to water vapour. These fractions are determined from the model via
The relevant data are given in Table B3.
Using these parameter settings, we minimized the sum of squares of the differences between the data from Table B3 (after nondimensionalization) with the model outputs allowing the parameters k_{S}, k_{C}, k_{W}, k_{Cl}, b, and C_{D} to vary. In the minimization calculation, the terms associated with the contribution values (last three columns of Table B3) were given a heuristic weight of 0.01, since these values are less reliable than the other data. The resulting calibrated values for these parameters are given in Table B2; the results of the fitting are given in Table B3 and Fig. B1. As can be seen in Table B3, the minimization achieved very good agreement with the globally averaged data.
Using the calibrated values for k_{S}, k_{C}, k_{W}, k_{Cl}, b, and C_{D} obtained from the first step, the second calibration step was to select the parameters for the functions $\mathit{\varphi}\left(\widehat{z}\right)$ and $\mathit{\psi}\left(\widehat{z}\right)$ to attempt to match the annual temperature profile for the Arctic from Cronin and Jansen (2016, Fig. 1). Data from that figure of their paper are reproduced in Table B4. The Arctic values of the geographicdependent parameters from Table B2 were used. Using each of the four forms for ϕ, given by Eqs. (A55)–(A58), and each of the two forms for F_{A}, given by Eqs. (A52) and (A53), the sum of the square of the differences between the model and the data in Table B4 was minimized by allowing the parameters z_{c}, Φ_{B}, L_{ϕT}, L_{ϕB}, and L_{ψ} to vary. The fit quality is shown in Fig. B2. From this figure it is evident that the first two forms for ϕ, namely Eqs. (A55) and (A56), do not give adequate fits. The other two forms for ϕ are similar, and both forms of ψ only give small changes. The best fit is the third form for ϕ and the second form for ψ, that is, Eqs. (A57) and (A53), and so these forms were chosen for the model. The calibrated parameters for these forms of ϕ and ψ are given in Table B2, and the corresponding functions ϕ and ψ are shown in Fig. B3.
For the above calibrations, the albedo values α_{c} and α_{w} in Eq. (A23) were set to the same constant value ($\mathrm{24}/\mathrm{185}$ for calibration to global data, $\mathrm{2}/\mathrm{3}$ for calibration to Arctic data) in order to ensure the empirical albedo value was matched regardless of the surface temperature. Now that the calibrations are complete, we wish to apply the model to the Arctic using different values of the carbon dioxide level. As such, we must choose values of α_{c}, α_{w}, and ω that allow for the albedo to change appropriately with surface temperature. The empirical albedo of $\mathrm{2}/\mathrm{3}$, which was used in the Arctic calibration, resulted in a surface temperature of T_{S}=253.4 K. Since this temperature is well below freezing, the albedo at this temperature should be near the maximum albedo, so α_{c} was set to 0.667. Second, α_{w} was set to 0.1, corresponding to the fact that, north of 70^{∘} latitude, the Earth is mostly oceancovered (albedo 0.06 for open water) and partly land (albedo 0.1–0.4). The value of ω was then calculated from Eq. (16) so that $\mathit{\alpha}\left(\mathrm{253.4}\right)=\mathrm{2}/\mathrm{3}$. This resulted in ω=0.01942.
B2 Insolation
This section presents the calculation of the insolation used in the Arctic model, where, in particular, the insolation is taken as an annual average over the region north of 70^{∘} N.
Select a Cartesian coordinate system $(x,y,z)$ for the solar system with the Sun at the origin, with the z axis perpendicular to the Earth's orbital plane and with the positive x axis defined by the direction in the orbital plane from the Sun to the centre of the Earth when the Earth's North Pole is furthest from the Sun (Northern Hemisphere winter solstice). Let (r,θ) be the usual polar coordinates for the centre of the Earth on the orbital plane. Approximate the incoming solar radiation to the Earth as parallel rays travelling from the direction $q=[\mathrm{cos}\mathit{\theta},\mathrm{sin}\mathit{\theta},\mathrm{0}{]}^{T}$ with energy flux S_{0}=1366 W m^{−2} (the solar constant). Let ϕ and ψ be the latitude and “longitude” of a location on the Earth's surface, where we assume that ψ=0 is aligned with the positive x axis of the solar coordinate system, not with some fixed location on the Earth's surface. If the Earth's axis of rotation was parallel to the z axis (no tilt), then the unit outward normal to the Earth's surface in the solar coordinate system would be
However, the Earth's axis is actually tilted by an angle $\mathit{\beta}=(\mathrm{23.5}/\mathrm{180})\mathit{\pi}$ in the negative sense in the (x,z) plane. Applying this tilt gives the unit outward normal as
The insolation striking the Earth's surface at (ϕ,ψ) is then ${S}_{\mathrm{0}}max(q\cdot n,\mathrm{0})$, where the maximum is due to the fact that the dot product is negative for points on the dark side of the Earth, away from the Sun, and hence the insolation there is zero.
Let S be the region of the Earth between latitudes ϕ_{1} and ϕ_{2}. The area element is dS=R^{2}cos ϕdψdϕ, where R is the radius of the Earth. Therefore the average annual (θ runs 0 to 2π) insolation in a region, S, of the Earth is
Numerical integration of the above formula yields Q=341 W m^{−2} for the entire globe, which agrees well with Wild's value of 340 (the difference is likely due to some ambiguity in the precise value of S_{0}), and yields Q=185 W m^{−2} for the Arctic region north of 70^{∘} latitude. If one considers the limit as the region of interest shrinks to size zero around the North Pole, the value of the insolation approaches the limit 173.85 W m^{−2}.
Analysis code is available from the authors on request.
All of the data used in this paper are publicly available from the references listed.
The research project was conceptualized by WFL. Funding acquisition was by WFL, GML, and ARW. The methodology was developed by KLK, GML, and ARW. KLK did the formal analysis, investigation and validation. Software was primarily written by KLK with support from GML and ARW. Figure visualization was done by KLK and ARW. ARW wrote the original draft; all the authors were involved in reviewing and editing.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We were supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). Kolja L. Kypke thanks the Ontario Ministry of Colleges and Universities and the University of Guelph for a Queen Elizabeth II Graduate Scholarship in Science and Technology.
This research has been supported by the Natural Sciences and Engineering Research Council of Canada (grant nos. 400450 and 006257).
This paper was edited by Vicente PerezMunuzuri and reviewed by Marek Stastna and one anonymous referee.
Abbot, D. S., Silber, M., and Pierrehumbert, R. T.: Bifurcations leading to summer Arctic sea ice loss, J. Geophys. Res.Atmos., 116, https://doi.org/10.1029/2011JD015653, 2011. a
Armour, K. C., Eisenman, I., BlanchardWrigglesworth, E., McCusker, K. E., and Bitz, C. M.: The reversibility of sea ice loss in a stateoftheart climate model, Geophys. Res. Lett., 38, L16705, https://doi.org/10.1029/2011GL048739, 2011. a
Årthun, M., Onarheim, I. H., Dörr, J., and Eldevik, T.: The Seasonal and Regional Transition to an IceFree Arctic, Geophys. Res. Lett., 48, e2020GL090825, https://doi.org/10.1029/2020GL090825, 2021. a
Bathiany, S., Notz, D., Mauritsen, T., Raedel, G., and Brovkin, V.: On the Potential for Abrupt Arctic Winter Sea Ice Loss, J. Climate, 29, 2703–2719, https://doi.org/10.1175/JCLID150466.1, 2016. a
Björk, G. and Söderkvist, J.: Dependence of the Arctic Ocean ice thickness distribution on the poleward energy flux in the atmosphere, J. Geophys. Res.Oceans, 107, 371–3717, https://doi.org/10.1029/2000JC000723, 2002. a
Cronin, T. W. and Jansen, M. F.: Analytic RadiativeAdvective Equilibrium as a Model for HighLatitude Climate, Geophys. Res. Lett., 43, 449–457, https://doi.org/10.1002/2015GL067172, 2016. a, b, c, d, e
Dortmans, B., Langford, W. F., and Willms, A. R.: An energy balance model for paleoclimate transitions, Clim. Past, 15, 493–520, https://doi.org/10.5194/cp154932019, 2019. a
Eisenman, I.: Factors Controlling the Bifurcation Structure of Sea Ice Retreat, J. Geophys. Res., 117, D01111, https://doi.org/10.1029/2011JD016164, 2012. a
Eisenman, I. and Wettlaufer, J. S.: Nonlinear Threshold Behavior During the Loss of Arctic Sea Ice, P. Natl. Acad. Sci. USA, 106, 28–32, https://doi.org/10.1073/pnas.0806887106, 2009. a, b, c
Flato, G. M. and Brown, R. D.: Variability and climate sensitivity of landfast Arctic sea ice, J. Geophys. Res.Oceans, 101, 25767–25777, https://doi.org/10.1029/96JC02431, 1996. a
Intergovernmental Panel on Climate Change: Climate Change 2013: The Physical Science Basis, in: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G. K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge Univ. Press, Cambridge, UK and New York, NY, USA, https://doi.org/10.1017/CBO9781107415324, 2013. a
Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deavan, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K. C., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R., and Joseph, D.: The NCEP/NCAR 40Year Reanalysis Project, B. Am. Meteorol. Soc., 77, 437–471, https://doi.org/10.1175/15200477(1996)077<0437:TNYRP>2.0.CO;2, 1996. a
Kuznetsov, Y.: Elements of Applied Bifurcation Theory, 3rd edn., SpringerVerlag, New York, ISBN: 9780387219066, 2004. a
Kypke, K. L., Langford, W. F., and Willms, A. R.: Anthropocene climate bifurcation, Nonlin. Processes Geophys., 27, 391–409, https://doi.org/10.5194/npg273912020, 2020. a, b, c
Langford, W. F. and Lewis, G. M.: Poleward expansion of Hadley cells, Can. Appl. Math. Quart., 17, 105–119, 2009. a
Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping Elements in the Earth's Climate System, P. Natl. Acad. Sci. USA, 105, 1786–1793, https://doi.org/10.1073/pnas.0705414105, 2008. a
Lewis, G. M. and Langford, W. F.: Hysteresis in a rotating differentially heated spherical shell of Boussinesq fluid, SIAM J. Appl. Dyn. Syst., 7, 1421–1444, https://doi.org/10.1137/070697306, 2008. a, b
Lutgens, F. K. and Tarbuck, E. J.: The Atmosphere: An Introduction to Meteorology, 14th edn., Pearson Education Inc., Boston, USA, ISBN: 9780134758589, 2019. a
Mayer, M., Tietsche, S., Haimberger, L., Tsubouchi, T., Mayer, J., and Zuo, H.: An Improved Estimate of the Coupled Arctic Energy Budget, J. Climate, 32, 7915–7934, https://doi.org/10.1175/JCLID190233.1, 2019. a, b
Merryfield, W. J., Holland, M. M., and Monahan, A. H.: Multiple Equilibria and Abrupt Transitions in Arctic Summer Sea Ice Extent, American Geophysical Union (AGU), 151–174, https://doi.org/10.1029/180GM11, 2008. a
MüllerStoffels, M. and Wackerbauer, R.: Regular network model for the sea icealbedo feedback in the Arctic, Chaos, 21, 013111, https://doi.org/10.1063/1.3555835, 2011. a
National Centers for Environmental Protection/National Weather Service/NOAA/U.S. Department of Commerce: NCEP/NCAR Global Reanalysis Products, 1948–continuing, updated monthly, Research Data Archive [data set], https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html (last access: November 2021), 1994. a
Pierrehumbert, R. T.: Principles of Planetary Climate, Cambridge University Press, Cambridge, UK, ISBN: 9780521865562, 2010. a, b
Ridley, J., Lowe, J., and Simonin, D.: The demise of Arctic sea ice during stabilisation at high greenhouse gas concentrations, Clim. Dynam., 30, 333–341, https://doi.org/10.1007/s0038200702914, 2008. a
Russill, C.: Climate Change Tipping Points: Origins, Precursors, and Debates, WIREs Clim. Change, 6, 427–434, https://doi.org/10.1002/wcc.344, 2015. a
Schmidt, G. A., Ruedy, R. A., Miller, R. L., and Lacis, A. A.: Attribution of the presentday total greenhouse effect, J. Geophys. Res., 115, D20106, https://doi.org/10.1029/2010JD014287, 2010. a, b
Schröder, D. and Connolley, W. M.: Impact of instantaneous sea ice removal in a coupled general circulation model, Geophys. Res. Lett., 34, L14502, https://doi.org/10.1029/2007GL030253, 2007. a
Serreze, M. C., Barrett, A. P., Slater, A. G., Steele, M., Zhang, J., and Trenberth, K. E.: The largescale energy budget of the Arctic, J. Geophys. Res., 112, D11122, https://doi.org/10.1029/2006JD008230, 2007. a
Stroeve, J. C., Markus, T., Boisvert, L., Miller, J., and Barrett, A.: Changes in Arctic melt season and implications for sea ice loss, Geophys. Res. Lett., 41, 1216–1225, https://doi.org/10.1002/2013GL058951, 2014. a, b
Thorndike, A. S.: A toy model linking atmospheric thermal radiation and sea ice growth, J. Geophys. Res.Oceans, 97, 9401–9410, https://doi.org/10.1029/92JC00695, 1992. a
Tietsche, S., Notz, D., Jungclaus, J. H., and Marotzke, J.: Recovery mechanisms of Arctic summer sea ice, Geophys. Res. Lett., 38, L02707, https://doi.org/10.1029/2010GL045698, 2011. a
van Groesen, E. and Molenaar, J.: Continuum Modeling in the Physical Sciences, SIAM, Philadelphia, USA, ISBN: 9780898716252 2017. a
van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J.F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J., and Rose, S. K.: The Representative Concentration Pathways: An Overview, Clim. Change, 109, 5–31, https://doi.org/10.1007/s105840110148z , 2011. a
Wild, M., Folini, D., Schär, C., Loeb, N., Dutton, E. G., and KönigLanglo, G.: The global energy balance from a surface perspective, Clim. Dynam., 40, 3107–3134, https://doi.org/10.1007/s0038201215698, 2013. a, b, c
Winton, M.: Does the Arctic sea ice have a tipping point?, Geophys. Res. Lett., 33, L23504, https://doi.org/10.1029/2006GL028017, 2006. a
Winton, M.: Sea Ice–Albedo Feedback and Nonlinear Arctic Climate Change, American Geophysical Union (AGU), 111–131, https://doi.org/10.1029/180GM09, 2008. a
Zheng, L., Cheng, X., Chen, Z., and Liang, Q.: Delay in Arctic Sea Ice FreezeUp Linked to Early Summer Sea Ice Loss: Evidence from Satellite Observations, Remote Sens., 13, 2162, https://doi.org/10.3390/rs13112162, 2021. a