Climate bifurcations in a Schwarzschild equation model of the Arctic atmosphere

. 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 ﬂux of longwave radiation is governed by the two-stream 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 two-point 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 ice-free 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, ice-free state.

bifurcations is well-developed (Kuznetsov, 2004) and employed here. Figure 1 illustrates the typical behavior associated with 25 these bifurcations. In Fig. 1(a) there are two saddle-node 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 stable state, then, as µ increases through a critical value µ crit , there is an abrupt jump to the upper stable state. 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 30 the original climate state. Fig. 1(b) illustrates the situation of a cusp bifurcation, where the two saddle-node bifurcations in (a) 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. 1(c), even though the saddle-node bifurcations may be outside the parameter interval of interest, abrupt large transitions in the system can result from a small noise-or perturbation-induced change to the system even when the parameter value remains constant. It is the presence of 35 saddle-node bifurcations in a mathematical model, even if not occurring precisely at the system's current parameter value, that are the root causes of all of the behaviors 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 40 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 freeze up in the autumn. The freeze up date for the Beaufort, Chukchi, Laptev, and Kara Seas, for example, has been getting later by 6-11 days 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% 45 per decade (Stroeve et al., 2014). Thus the observational evidence indicates that the processes behind this phenomenon are not linear at all, but nonlinear.
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 flow-invariant, 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, 70 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 models of Dortmans et al. (2019), which was applied to paleoclimate transitions, and 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 75 atmosphere with absorption properties similar to the real atmosphere, here we use the Schwarzschild two-stream equations to model absorption in the atmosphere explicitly as a function of altitude (Pierrehumbert, 2010, pg. 191).
A bifurcation analysis is performed on the model, tracking the steady-state 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 saddle-node bifurcation would result in an abrupt change in climate that may be practically irreversible. 80 The simpler model of Kypke et al. (2020) showed this behaviour under certain CO 2 representative concentration pathway scenarios. We seek here to determine if 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 is presented in Appendix B. Section 3 presents the results, and the conclusions are in Section 4.

85
The model is developed from first principles and has the following features.
-The atmosphere is a one-dimensional 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.

90
-The surface albedo is a nonlinear function of the surface temperature.
-A well-mixed surface boundary layer is included.
-The Earth emits longwave radiation as a black body.
-The atmosphere is considered to be a gray gas.
-The Schwarzschild two-stream equations govern the absorption and emission of both upward and downward directed 95 longwave radiation in the atmosphere.
-The atmospheric absorption of longwave radiation is due to three factors: water vapor, CO 2 concentration, and clouds.
-Water vapor 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 transport to the Arctic are dictated by empirical values.

100
-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 analytically, and a two-point 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 115 and is simply considered as leaving the system.

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 time derivative 120 term will be taken as zero since only the steady-state 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.

Mass
If the property f in (1) is the mass density, ρ (kg m −3 ), then the flux is χ = ρw. There is mass flux across the vertical boundary 125 of the cylinder, M b (z) ( kg m −2 s −1 ), which is assumed to be immediately spread out evenly across the layer, hence the mass flux across the vertical boundary in the model is really a mass source term in the interior giving S = mass entering into cylinder layer of width ∆z volume of layer Thus, at steady-state the mass balance equation is The mass flux through the vertical boundary into the column is written as where M max (kg m −2 s −1 ) is a nonnegative constant and φ(x) : [0, 1] → 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 cross sectional area of the column to the area of its side, 5 https://doi.org/10.5194/npg-2022-2 Preprint. Discussion started: 19 January 2022 c Author(s) 2022. CC BY 4.0 License.
The first of these conditions dictates that there is no net mass entering the system, and the second is a normalization condition 140 so that M max represents both the total mass entering the column per unit cross sectional area, and the magnitude of the total mass leaving the column per unit cross sectional 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 145

Momentum
Now take the property f in (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, pg. 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 dP dz = −ρg, which is the hydrostatic equation.)

155
Finally consider the case where the property in (1) is the total energy density given by e = 1 2 ρw 2 + ρgz + c v ρT, which corresponds to the sum of kinetic energy, gravitational potential energy, and internal heat energy densities. Here c v (J kg −1 K −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 (W m −1 K −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 terms due to radiation (shortwave downward, and longwave upward and downward), one due to latent and sensible heat transport, 165 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 170 that it generates the average slow movement of air downward near the north pole (about 1 mm/s), 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 Section A3. Thus S is given by

175
The governing equations for the longwave radiation intensities are the two-stream Schwarzschild equations and for the shortwave radiation a standard absorption equation: 180 where k S (m 2 kg −1 ) is the shortwave absorption coefficient, σ (W m −2 K −4 ) is the Stefan-Boltzmann constant, and κ (m −1 ) is the long wave 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 CO2 and M A (kg mol −1 ) are the molar masses of 185 CO 2 and dry air, respectively, δ z−z B z T −z B 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 190 of vaporization for water, and R W = R/M W (J K −1 kg −1 ) is the gas constant for water, R (J K −1 mol −1 ) is the universal gas constant, and M W (kg mol −1 ) is the molar mass of water. The corresponding density is ρ sat W (T ) = P sat (T )/(R W T ) 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. Therefore the energy balance equation at steady-state is d dz 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 A = R/M A (J kg −1 K −1 ) is the gas constant for air.
The mass, momentum, and energy balance equations, (3), (4), and (11), along with the Schwarzschild equations (5)-(6), and the equations governing shortwave absorption (7) and sensible and latent heat transport (10), are the differential equations for the BVP on [z B , z T ]. Equations (8), (9), and (12) define certain quantities in these differential equations. Before discussing the 205 boundary conditions for the BVP it is necessary to consider the surface boundary layer.

Surface Boundary Layer
The model includes a boundary layer extending from z = 0 to z = z B . It is assumed that this layer is well-mixed so that temperature T B = T (z B ) and density ρ 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 .

210
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 introdcution of the surface boundary layer avoids this singularity.

215
The total mass crossing from the atmosphere into the boundary layer per unit time is M max Φ 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 of potential energy 220 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 225 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 230 the form of sensible and latent heat, which is modelled, as per Pierrehumbert (2010, pgs. 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 (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 235 surface is where is the surface albedo, here modelled as a sigmoid function increasing from α c at cold temperatures to α w at warm temperatures, 240 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 short and longwave radiation entering at z B . Output 245 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 (13). Therefore the energy density balance for the combined surface and boundary layer is Since temperature and pressure 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 with those at z = z B . The simple ODE for F C is also easily solved in it is equal to the black body radiation of the surface, σT 4 S . The initial condition for that latent heat, F C (0), is given by (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)-(19) provide two boundary conditions for the BVP on the troposphere. The energy balance equations (15), (17) along with equations (16), (20), and (21) provide two further boundary conditions.

270
There are eight unknown dependent variables: w, ρ, I + , I − , I S , F C , T , and dT dz ; the pressure P can be written in terms of the others via (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 max Φ B A, pressure at the surface equal to the standard pressure, upward longwave radiation at z B given by (18), vertical heat transport at z B given by (19), the energy balance equations (15) and (17) with expressions from (20) and (21) where F C0 is given by (14). The last three terms of Eqn. (27) may be simplified using Eqn. (22) so that they read There are nine boundary conditions, but there are only eight dependent variables for which we have differential equations. The 290 discrepancy is explained by the presence of T S , which is an additional scalar variable. 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 295 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 is provided in Appendix B.

Results
For the Arctic parameter values given in Appendix B and for a given CO 2 concentration, µ, the model can be solved numeri- 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.

310
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.   RCPs are RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.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 RCP 8.5. These RCPs are plotted in the lower panel of follow a trajectory similar to RCP 8.5, Arctic climate may change drastically in less than 100 years, but a return to the current 350 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 tot A , and ocean heat transport, F O , are inputs to the model that are empirically based and are not altered by the state of the model itself. We investigated the persistence of the saddle-node bifurcation and hysteresis in the presence of the uncertainty of these two values. Figure 4 shows two two-dimensional bifurcation diagrams plotting the The right curve in each of the panels is a curve of saddle-node bifurcations corresponding to the right bend in the S-curve in Fig. 3; they are transitions from a cold state to a warm state. Similarly, the left curve in each of the panels is a curve of saddle-node bifurcations that correspond to the left bend in the S-curve, and are the transitions from a warm state to a cold 360 state. The narrow 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 saddle-node bifurcations and the hysteresis phenomenon will be present for all physically possible CO 2 levels, and for all reasonable values of F tot A and F O . These bifurcation curves also make it evident that a transition from a cold to a warm state occurs if either F tot A or F O are increased, even if CO 2 levels are constant. If F tot A and F O increase along with carbon dioxide levels in the near future, which 365 seems reasonably likely, then the saddle-node bifurcation will occur at lower levels of CO 2 . For example, a ten percent increase of F tot A to 110 W m −2 from 100 would cause the saddle-node bifurcation to occur at about µ = 754 ppm, which would mean that both RCP8.5 and RCP6.0 would pass through the saddle-node bifurcation.

Conclusions
Although the model presented here is clearly a simplification of the climate, made possible by the near invariance of the vertical 370 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 RCP 6.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 RCP 6.0, such a change is likely, and the cause is a saddle-node 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 375 Å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 for a high emissions scenario it will be year-round 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 RCP 6.0 and RCP 4.5 lie in the region of bistability from the year 2070 onward. In a bistable situation, external 380 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 year's time. This bistability has been shown to persist regardless of the values of the two biggest "unknowns" in the model, the atmospheric and ocean heat transport to the Arctic from the mid-latitudes.

385
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. Fig. 3 suggests that such temperature oscillations would not likely need to be very large to 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 year-round ice-free Arctic is not likely to occur while the Arctic is ice-covered for a significant portion of the year, but once the Arctic is seasonally ice-free 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 like in Eisenman and Wettlaufer (2009).

395
This appendix provides details regarding the model. The model is written in a nondimensional form in Section A1, defining nondimensional variables and parameters. This system is then transformed in Section A2 to a standard form of nine first order 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 400 dependence of relative humidity, atmospheric meridional heat transport, and mass flux on altitude.

A3.2 Atmospheric heat transport
Atmospheric heat transport is primarily due to large scale turbulent mixing of the column with its environment. This mixing is 500 not modelled explicitly, but the thermal energy supplied to the column by it is represented by the function F A (z). 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 tot A we have This provides one restriction on the function F A (z), but its precise form is not prescribed.

505
Getting reasonable results from the model seems to require a delicate balance for F A , especially near the tropopause. If it is too small, then the temperature drops precipitously toward absolute zero; if it is too big, the temperature turns around and climbs rapidly. This behaviour is sensitive to F A (z) for z near z T , that is,ẑ near 1. In order to automate the choice of F A (z) we have proceeded as follows: 1. ChooseF A (1) such that the temperature gradient at the tropopause is zero, that is, re-impose boundary condition (A47).

Assume that F
where the base function F Ab (z) is a linear function passing through zero at the midpoint of the atmosphere and equal to 515 F A0 at z T , that is, 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 assure the temperature gradient at the top is zero.
so that the value of F A is not altered at the tropopause, and so that F Amax (W m −2 ) represents the total energy flux of atmospheric heat transport entering the column. Thus the non-dimensional function iŝ 3. As a numerical issue, since the boundary condition valueF A0 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 dy10 dẑ = 0 and boundary condition y 10 (1) =F A0 . is somewhat open, we tested the following two forms: where the functions g 1 and g 2 are defined as

535
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 φ(ẑ).) The primary difference between these two forms is that the first has a non-zero 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 540 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 (2). To model the situation in the Arctic, we want a downward flow of air with a vertical wind speed, w, on 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 we assume that where z c is some point in [0,1]. The following forms were tested for φ(ẑ): where g 1 and g 2 are defined by (A52). (In the case that z c = 0 or z c = 1, it is understood that only the non-empty interval for φ in the above definitions is used and that it is closed at both ends.)

Appendix B: Model Parameters and Calibration
This appendix lists the parameter values used in the model and discusses how some of them were calibrated to empirical data.
Section B1 gives the calculation of the average annual insolation for the Earth north of 70 N.
Values of most of the model parameters are given in Tables B1 and B2. The parameters in Table B1 are physical constants except for the last two. The surface horizontal wind speed, U , which is a factor in the sensible and latent heat transport from 560 the surface, was set to 10 m/s. The exact value is not too important since U always appears multiplied by the factor C D , which we calibrate to data below. The height of the boundary layer was set to z B = 50 m. The model is not very sensitive to this parameter.
The parameters in Table B2 are those that depend on geographic location. 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. The height of the 565 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 Section B1. For the average global situation there is no ocean or atmospheric heat transport; also come from Mayer et al. (2019). Relative humidity is low at the top of the troposphere so in both the global and Arctic cases was set to 10%. The surface relative humidity was set to 75% for the global average and 70% for the Arctic. 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 max = 2.0 × 10 −6 kg m −2 s −1 , and we set Φ B = −1, Φ T = 0.2, and assumed φ(ẑ) is given by (A55) 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 on the 575 order of 10 −3 mm/s, which is far enough away from the singularity to avoid convergence issues, but is small enough so all of the convection-related terms in the model become negligible. Mass flux for the Arctic situation was set with trial and error to 8.0×10 −4 kg m −2 s −1 , which gave vertical wind speeds in the column on the order of 0.5 mm/s. Since F tot A = 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 2 , and k Cl , 580 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 T ot = 1 0κ (y 2 , y 7 )y 4 dẑ, C contrib = 1 T ot 1 0 G C µy 2 y 4 dẑ,

585
Cl contrib = 1 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 2 , 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 590 weight of 0.01, since these values are less reliable than the other data. The resulting calibrated values for the parameters are given in Table B4; the results of the fitting are given in Table B3 and Figure 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 2 , k Cl , b, and C D obtained from the first step, the second calibration step was to select the parameters for the functions φ(ẑ) and ψ(ẑ) to attempt to match the annual temperature profile for the Arctic from 595 Cronin and Jansen (2016, Fig. 1). Data from that figure of their paper is reproduced in Table B5. The Arctic values of the geographic-dependent parameters from Table B2 were used. Using each of the four forms for φ, given by Eqns. (A53)-(A56), and each of the two forms for F A , given by Eqns. (A50)-(A51), the sum of the square of the differences between the model and the data in Table B5 was minimized by allowing the parameters z c , Φ B , L φT , L φB , and L ψ to vary. The fit quality is shown in Figure B2. From this figure it is evident that the first two forms for φ, namely equations (A53) and (A54) do not give adequate 600 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, equations (A55) and (A51) and so these forms were chosen for the model. The calibrated parameters for these forms of φ and ψ are given in Table B6 and the corresponding functions φ and ψ are shown in Figure B3.
For the above calibrations, the albedo values α c and α w were treated as the same constant 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    7.60e-04 Figure B2. Temperature profiles for best fits for each of the φ (rows) and ψ (columns) forms. The subfigure labels φ = m for m = 1, 2, 3, 4, correspond to Eqns. (A53)-(A56), respectively, and the labels ψ = n for n = 1, 2, correspond to Eqns. (A50)-(A51), respectively. The data from Cronin and Jansen (2016) are the red circles. The numbers in the bottom left are the residuals for the fits.
model to the Arctic allowing for the albedo to change with surface temperature. The values of α c , α w , and ω were chosen as follows. The empirical albedo of 2/3 was used in the Arctic calibration, which 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 ocean    Table B7.

B1 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 615 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 traveling from the direction q = −[cos θ, sin θ, 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 620 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 The insolation striking the Earth's surface at (φ, ψ) is then S 0 max (q · n, 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.
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.
Competing interests. The authors have no competing interests.
Acknowledgements. We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). KK