Polar Jet Stream Fluctuations in an Energy Balance Model

We investigate the effect of increased longwave radiative forcing (a proxy for increased greenhouse gas concentration) on the zonally averaged location of the eddy-driven jet stream in a latitude dependent, two-layer Energy Balance Model. The model includes separate terms for atmospheric and surface albedos, and takes into account reflections of shortwave radiation between the surface and atmospheric layers. We introduce the notion of a cloud factor function, which depends on temperature gradients, to simulate the eddy-driven jet. An increase in longwave radiative forcing initially results in a poleward movement of the jet stream's mean latitude, but as the forcing increases, the location of the jet stream becomes quasi-periodic and its mean location moves equatorward.


Introduction
The atmosphere and the ocean stabilize Earth's climate from uneven solar insolation by transporting heat from the equator to the poles.Energy balance models (EBMs), first introduced by Budyko [6], Sellers [40], include heat transport terms that reproduce zonally and annually averaged temperature profiles from this transport.These idealized climate models have been extensively studied (e.g.North [32,33,34]) and a wide range of modifications and additional forcings have been introduced in order to provide insights into causal relationships of components of Earth's climate, for example, [17,11,19,20,30,23,45,42,4,12] among many other studies.
In this paper, we use an energy balance model to investigate the dynamics of the polar jet stream of an aqua-planet in response to increasing greenhouse gas concentrations, with a focus on the role of cloud fraction and albedo.Both observations and climate model studies indicate that the general circulation pattern of the atmosphere is altered by anthropogenic warming, e.g., [3,14,18,22,26,27,51,52].Among these are two studies that employed EBMs to investigate the link between shifts of the midlatitude storm tracks to the shifts of the Hadley cell edge: Mbengue and Schneider [29] (hereafter MS18) and and Siler, Roe, Armour [41] (hereafter SRA18).
MS18 [29] defined the storm track in a one layer EBM as the latitude of maximum absolute value of the temperature gradient.In that model, the diffusion coefficient was increased within the Hadley cell, relative to the diffusion coefficient outside the cell, and the Hadley cell edge (or terminus) was interactive and also depended on the convective lapse rate γ in the tropics, which was treated as a parameter.The model predicts that storm tracks shift in tandem as the Hadley cell edge is moved poleward by decreasing γ.Their results also indicate that strengthening meridional temperature gradient at the Hadley cell terminus can reduce the distance between the Hadley cell edge and the storm tracks, resulting in storm tracks that do not parallel shifts of the Hadley cell terminus.SRA18 [41] studied a single layer Moist Energy Balance perturbation model.Assuming a reference climate determined by reanalysis or averages of climate models, their perturbation model determines a change in temperature and in evaporation minus precipitation, E − P , as a function of latitude, from forcings such as increased greenhouse gas concentrations.The extratropical latitude of the minimum value of E − P serves as the proxy for mid-latitude storm tracks.In the case of spatially uniform radiative forcing, SRA18 [41] found that downgradient energy transport implies a poleward expansion of the subtropics where E − P > 0, and a poleward shift in the extratropical minimum of E − P , consistent with a poleward shift of storm-track latitudes.
The idealized model considered in this paper is a latitude dependent, two-layer energy balance model that includes separate terms for atmospheric and surface albedos, and takes into account reflections of shortwave radiation between the surface and atmospheric layers, and includes heat diffusion terms for each layer.The novel feature of our model is what we refer to as a "cloud factor function", a function which depends on temperature gradients, and which dynamically simulates the eddy-driven or polar jet stream.More specifically, at any fixed time, the cloud factor function, C f (θ), is a dimensionless quantity that represents the fraction of the zonally averaged planetary albedo at latitude θ attributable to clouds.We use it to construct the atmospheric albedo as a function of latitude at each time step in our model (see Section 2.1 below).
The thermal wind equations link the horizontal temperature gradient to the polar-front jet and suggest that the location of the jet may be identified with the location of the maximum magnitude of the extratropical temperature gradient; this proxy was utilized in MS18 [29].Similar to MS18 [29], we interpret the latitude where this occurs as the averaged location of the eddy-driven jet, and define our cloud factor function to achieve a maximum value at that location at each time step in our numerical scheme.This allows us to track location of the jet as it moves dynamically until the system reaches equilibrium.
We must point out that the Hadley cell edge is not interactive in our model.We hold it fixed at 30 • latitude in our numerical experiments.However, this location can easily be modiflied, and the qualitative behavior of our model is robust with respect to this location.Despite this constraint, our model identifies a driver of jet stream fluctuations which has the potential to be incorporated into more complex climate models that include Hadley cell dynamics.
The cloud factor function, C f (θ), is constructed so that its minimum corresponds to the Hadley cell boundary and so that the lowest extratropical latitude, at which C f (θ) reaches a prescribed maximum value, identifies the mean location of the polar jet stream.We assume that the cloud factor function is given by a cubic Hermite spline.This spline is defined by specified values of the cloud factor function at four latitudes -the equator, the Hadley cell edge, the polar jet stream, and the pole.These values are fixed, but one of the latitudes -the polar jet stream latitude -is a function of the temperature gradient.There is thus one degree of freedom in the cloud factor function.The location of the jet is determined by the gradient of the average of the atmospheric and surface temperatures.
This paper is organized as follows.Section 2 is divided into subsections that describe the components of our model, including standard forcings, but which focus primarily on the couplings between the cloud factor function and the surface and atmospheric albedos.We also describe how the latitude, where the maximum magnitude of the temperature gradient occurs at each time step of our computations, alters the cloud factor function for the next time step.Section 3 describes the results of numerical experiments for changes in the location of the eddy-driven jet as radiative forcing increases, such as from increasing greenhouse gas concentrations.In Section 4, we compare the behavior of our model with other investigations of jet stream response to increasing greenhouse gas concentrations and offer concluding remarks.In addition, there are three appendices.Appendix A gives an explicit formula for the cloud factor function; Appendix B provides a concise description of the numerical scheme used in our computations; and Appendix C displays output data.

Model Description
Our EBM consists of an ocean covered surface layer and an overlying atmospheric layer.Throughout, we let x = sin ϕ, where ϕ is latitude4 , so that −1 ≤ x ≤ 1, but because our aqua-planet is symmetrical, we will generally display data only for the northern hemisphere, 0 ≤ x ≤ 1.
Let T s and T a represent the zonally averaged temperatures of the surface and atmosphere respectively, expressed in degrees Celsius.Here, T a is a measure of the free tropospheric temperature, say at 500 hPa, but as in [38] we express it as an equivalent surface air temperature, assuming a constant lapse rate (depending on the value of parameters used in the model) 5 .The time evolution of the temperatures are solutions to coupled differential equations of the form, where a is the radius of Earth, C a , C s are respectively specific heats of the atmosphere and surface, F out is the longwave radiative heat flux to space, and F up is the net flux of longwave radiation, latent heat, and sensible heat from the ocean to the atmosphere.The last terms in each equation represent meridional diffusive heat transport (given explicitly in Eqs (16a) and (16b) below).
Although multiple processes are involved in heat transport and although they vary across regions and time scales, Stone (1978) [46] demonstrated that the magnitude of the annual mean total meridional heat transport is insensitive to the details of dynamics of the atmosphere-ocean system.
As described below, these terms will be chosen to match the corresponding terms in the two layer energy balance model of Rose and Marshall [38,39].By contrast, the remaining two terms, F ↓ atm and F ↓ ground , in Eqs.(1a) and (1b) represent incoming solar radiation flux and both depend on the atmospheric albedo, α a , and ground albedo, α g .
To model the dependence of F ↓ atm and F ↓ ground on α a , and α g , we follow Qu and Hall [37] and Donohoe and Battisti [10].We assume an atmospheric layer within which the radiation undergoes three processes: reflection by a factor α a , transmission by a factor T sw (the transmissivity of shortwave radiation), and absorption by a factor A sw = 1 − α a − T sw .
Summing up the infinite number of transmissions and reflections between the atmosphere, the ground, the shortwave flux abosrbed by the ground F ↓ ground and the shortwave flux absorbed by the atmosphere F ↓ atm and the radiative flux to space from the top of the atmosphere are given by, (2) where s(x) is the annual weight function for incoming solar radiation (dimensionless, unit global mean) which, following [38,39], is given in terms of the second order Legendre polynomial P 2 (x) as, with s 2 = −0.48.
We note that F ↓ ground + F ↓ atm + F ↑ T OA = S 0 s(x)/4, that is, the sum of the various components of absorbed and reflected radiation equals the total quantity of incoming solar radiation.Eq. ( 4) does not appear in our EBM, but it shows that the planetary albedo can be identified as, where and T e α g can be considered as the contribution from the ground albedo α g to the planetary albedo modulated by the interactions with the atmosphere.

Cloud Factor Function
In order to assign latitudinal values to the ground and atmospheric albedos, α g and α a , we first introduce a cloud factor function, C f is a function of latitude θ and of the location θ(t) of the maximum of the absolute value of the temperature gradient.The function C f is related to the zonally averaged albedo at latitude θ attributable to clouds (see Eq.( 9) below).
The cloud factor function, whose general features are motivated by Figure 1, is explained in detail in Subsection 2.5 below.
It is difficult to measure cloud cover in the polar regions due to a number of factors, including thin and low lying clouds and polar conditions that create an unusual amount of near surface hazes and fogs [8].Cloud fraction in global climate models and atmospheric reanalyses vary widely [5], and clouds are among the main sources of uncertainty in modeling the Arctic climate [25].Because of these problems, there is an uncertainty in cloud cover over the polar regions.Vavrus et al. [50] conclude maximum cloudiness occurs over open water in the summer time with cloud fraction values of 81%.Palm et al. [35] agree that maximum cloudiness occurs over open water in the summer time but report model cloud fraction values of 90% or more.Both conclude that the average polar cloud fraction is increasing as the sea ice extent has been decreasing.
In [31] and references therein, Norris examined climate variability and found a positive cloud feedback on sea surface temperatures (SST), in the North Pacific during the boreal summer, where increased cloud amount acts to cool the ocean by decreasing surface insolation, and decreased SST favors greater marine stratiform cloudiness amount.This suggests a steep drop in temperature associated with high cloud cover.In addition, Figures 1, 2, and 3 in SRA18 [41] indicate maximum precipitation at the minimum of E − P with high values of precipitation poleward.This suggests high cloud cover poleward of the jet location.Taking these findings into consideration, we construct C f (θ, θ) through the use of cubic Hermite splines. 8The cloud factor function is incorporated into our climate model as described in Section 2.5.The graph of the cloud factor function is initially constrained to take extremal values at 0 • , 30 • , 50 • , and 90 • latitude, the locations of the equator, the Hadley cell edge, the polar jet stream, and the pole.Specifically, the coordinates are (0, 0.9), (30, 0.1), (50, 0.8) and (90, 0.8) so as to represent high cloudiness at the equator as well as poleward of 50 • degrees, and low cloudiness at 30 • degrees.However, as we explain in Sect.2.5, the graph will change with the time steps in the numerical runs of our model.A sample graph is shown in Figure 2. We note that our numerical experiments exhibit the same qualitative behavior, as we describe in this paper, even when the cloud factor function is modified so that the cloud cover varies in the region poleward of the jet or takes a different constant value in that region.
As pointed out in [28], the Southern Hemisphere polar jet is located at 50 • latitude.So this is a plausible choice for an initial location of the polar jet prior to radiative forcings that we will impose.We note that the EBM of SAR18 [41] locates the initial northern hemisphere minimum value of E − P (evaporation minus precipitation and the proxy in that EBM for storm track location) above 60 • latitude (see Figures 2f and 3b in [41]).

Albedo Functions
Our modeling of the atmospheric albedo α a and the ground albedo α g begins with an initial approximate estimate of the planetary albedo.As a reference frame and a guide, Figure 3 shows the zonal mean planetary albedo partitioned between atmospheric and surface components.In our model, we first approximate the total planetary albedo by choosing a reference planetary albedo α p0 of the form The coefficients α 0 p0 and α 1 p0 are chosen along with parameters for the ground albedo in Eq (11) so that the average planetary albedo approximates Earth's average planetary albedo, and in order to specify initial equilibrium locations of maximal absolute values of the temperature gradient (for further elaboration, see the third paragraph in Section 3). Figure 4 shows a plot of α p0 for this choice of parameters: α 0 p0 = 0.25 and α 1 p0 = 0.38.We emphasize that at no time step in our computational scheme does the function in Eq. ( 8) represent the planetary albedo in our model, which instead will vary in time in a way that depends on the global temperature distribution.We use α p0 , along with the cloud factor function C f , to define the atmospheric contribution to the planetary albedo as: where α clear is the clear sky (cloud free) albedo of the atmosphere which we take as constant, α clear = 0.149 [43].An initial sample plot of the atmospheric albedo is given in Figure 5.We note that α a depends on C f and, in turn, C f depends on temperature gradients, so α a depends on temperature gradients.We can now define the atmospheric transmittance of short wave radiation (SWR) in terms of α a as, where A sw = 0.05 is the atmospheric absorption of SWR [20].We note that T sw depends on C f , making it interactive.
Following other researchers (for example [21]), we model the ground albedo as a function of the surface temperature using the hyperbolic tangent function as follows, α g = 0.40 − 0.34 tanh(T s + 8). (11)

Albedo Constraint
The fraction of incoming solar energy sent back to space from Earth is about 29% [43] with roughly 88% of that coming from the atmospheric contribution and the remainder due to the modulated surface albedo [10,37].We therefore tune our model so that our initial atmospheric and modulated ground albedos are close to these values.They cannot be constrained in model runs because the atmospheric and ground albedo contributions in our model are dynamic and therefore fluctuate.
The total planetary albedo ᾱp is given by, ᾱp = 1 2 where, as before, x is the sine of latitude, α p (x) is the zonally averaged albedo at x given by Eq. ( 6), and s(x) is the annual weight function for incoming solar radiation given by equation 5.The planetary atmospheric albedo ᾱa is defined as, ᾱa = 1 2 where α a (x) is the zonally averaged atmospheric albedo at x. Therefore, we define the total planetary effective ground albedo by,

The Model
Our model is based on the energy balance equations given in this section.We begin by linearizing the terms F up , F out in Eqs (1a) and (1b) and write, Collecting the remaining terms from the preceding sections, the system of coupled PDEs for the zonally and column averaged two layer climate system becomes, Table 1 lists the parameter values for the constants in Eqs.(16a) and (16b).These are the same values in Rose and Marshall [38,39], except that our ocean heat capacity is greater by a factor of 10 in order to simulate a greater ocean depth.The value of C s in Table 1 together with the nominal values of heat capacity and density for water (as opposed to seawater) of 4184 J/kg/deg C and 1000 kg/m 3 assigns an ocean depth of approximately 24 meters.This value is shallow compared to observations of Earth's mixed layer depth [9], but the absence of land in our aqua-planet model is a compensating feature.At any rate, the qualitative behavior of our model is largely independent of the numerical value chosen for C s .
Parameter Units Numerical Value a m 6.373 Table 1: Parameter values for the EBM.
The initial (t = 0) temperature profile is specified below, and the dynamic feature of the cloud factor function are explained in the next section.
The system of equations Eq.( 16a) and (16b) is defined for −1 ≤ x ≤ 1, where x < 0 is the Southern Hemisphere and x > 0 is the Northern Hemisphere.But since the Southern and Northern hemispheres are symmetric (including our initial conditions), we need only consider the solution from 0 ≤ x ≤ 1.

Cloud Function Dynamics and Polar Jet Stream
The response of the eddy-driven jet to arctic amplification [16] and changing meridional temperature gradients has been analyzed extensively (e.g., [3,2,14,27,29,41,51] and references therein).With the thermal wind equations in mind, we identify the mean latitudinal position of the jet stream, at any time t, with the location of the maximum value of a meridional temperature gradient given by, We motivate this choice as follows.Let T (z) = T (θ, z) be the zonally averaged temperature at altitude z and fixed latitude θ.The vertically averaged temperature T = T (θ) at θ is given by, where h is the height of the troposphere, and the constant Γ is the zonally averaged lapse rate at θ. Thus, If we interpret T s = T (0) and T (0) − Γh, to be the atmospheric temperature at altitude h then from Eq.( 19), The gradient of T is then given by Eq.( 17).
We couple the temperature gradient (17) with the cloud factor function in the following way.We solve the model equations in Sect.2.4 numerically by time-stepping out to equilibrium (or quasi-periodicity).C f is updated at every timestep by setting it to 0.8 at the latitude of maximum temperature gradient.This choice, together with the cubic Hermite spline functional form and the specified equatorial, polar, and Hadley cell edge values, uniquely determines C f at all latitudes.Among other things, it entails that C f = 0.8 at all latitudes poleward of the maximum temperature gradient.
For example, the graph in Figure 2 corresponds to a maximum meridional temperature gradient occurring at 50 • latitude.Since the atmospheric albedo α a depends on C f (c.f.Eq.( 9)), it is updated at each time step.Similarly, the ground albedo α g (which is a function of latitude) is updated at each time step according to the values of the surface temperature T s in the previous time step (see Eq.( 11)).Numerical approximation details are described in Appendix B.

Numerical Results
In this section, we present numerical results from our model in response to increases in radiative forcing, such as from increased greenhouse gas concentrations.Following [38], to simulate this, we decrease the parameter A out which controls the flux of outgoing longwave radiation (OLR) from the top of the atmosphere.Our focus is on how the latitudinal locations of the maximum modulus of temperature gradient are affected by these increases in radiative forcing.We interpret those latitudes as the averaged locations of the polar jet stream.
Since the coupled partial differential equations of the model are non autonomous, equilibrium temperature and temperature gradient values for each experiment must be found by numerically running them out to equilibrium 9 .The results of this section take as initial temperature distributions the final equilibrium temperatures obtained by Rose and Marshall [38] (in their Figure 2), but the model behaviors are insensitive to the choice of initial temperature distributions.
To set a reference climate, we take A out = 214 Wm −2 .In equilibrium, this results in a climate with a planetary albedo, ᾱp = 0.298, and average temperatures given by T s = 14.4 • C and T a = 15.5 • C. The maximum absolute value of the atmospheric temperature gradient occurs at 55.4 • latitude.This is our proxy for the average latitude of the jet stream.The temperature and gradient distributions for A out = 214 Wm −2 are displayed in Figure 6(a) and Table 2.
By decreasing the parameter A out , we introduce a longwave radiative forcing in the model.Meridional profiles of T s and T a and temperature gradient plots for A out = 214, 213, 212, and 211 Wm −2 are shown in Figure 6.
Equilibrium is reached for the first three forcings, A out = 214, 213, 212 Wm −2 .However, for A out = 211 Wm −2 , the maximum absolute value of the atmospheric temperature gradient begins to exhibit oscillatory behavior.This is indicated by the red dots in Figure 7.As A out continues to decrease to values below 211 Wm −2 (so that radiative forcing increases), the modulus of the temperature gradient given by Eq.( 17) does not peak at a singular latitude, but instead produces a collection of nearly equal large values within an interval of latitudinal coordinates.As a physical interpretation, this suggests oscillatory behavior of the jet stream, and this is shown graphically in Figure 8, based on the data in Appendix C. Additional detail is shown in Figure 9, which displays plots of the temperature gradient within a narrower range of latitudes, and illustrates the formation of approximate plateaus of maximum values of |∂(T a + T s )/∂θ|.We note that numerical experiments show that the same oscillatory behavior appears when the time step is reduced to half days and quarter days (instead of days).
Tables 2 and 3 show the mean latitudinal locations of the jet, standard deviations from the means of the jet locations, along with temperature and albedo data, as A out decreases from 214 to 202 Wm −2 .The standard deviations reveal the extent of oscillations of the jet.As shown in Table 2, oscillations increase as A out decreases to 208 Wm −2 .The movement of the jet location, as the forcing increases, is initially poleward, but as the forcing increases further, the mean jet location begins to move equatorward.Table 3 shows that the mean jet location continues to move equatorward as the forcing increases (i.e., as A out decreases), and the standard deviation data indicates that the jet oscillations decrease and nearly cease at A out = 206 Wm −2 , and lower values, so that the mean jet location is nearly constant for each of those values.
The picture that emerges is that significant oscillations of the jet occur only for the band of A out values between 211 and 207 Wm −2 , and the mean location of the jet increases poleward from 55.4 • latitude for A out = 214 Wm −2 to 62.3 • latitude for A out = 212 Wm −2 , and thereafter moves equatorward.
The climate sensitivity of our model can be determined from the temperature data in Tables 2 and 3.For the purpose of comparison, we first note that the IPCC's AR6 estimate for Earth's modern Effective Radiative Forcing (ERF) for a doubling of atmospheric CO 2 is 3.93 ± 0.47 Wm −2 , and 3.73 ± 0.44 Wm −2 for the stratospherically adjusted radiative forcing.The equilibrium climate sensitivity is estimated to be 3 • C [13].
Tables 2 and 3 show that the climate sensitivity of our model varies with temperature and forcing.This is not unprecedented.In their study of climate sensitivity in the context of high temperature and large radiative forcings, Caballero and Huber [7] which amounts to a warming of 3.1 • C from the IPCC's reported effective radiative forcing of 3.9 Wm −2 for a doubling of atmospheric CO 2 .

Discussion
Our results may be compared with observations and predictions from more elaborate models.Using the Coupled Model Intercomparison Project (CMIP5) and assuming the representative concentration pathway 8.5 (RCP8.5)scenario, Barnes and Polvani [3] found that all jets migrate poleward in the twenty-first century.Using reanalysis, Manney and Hegglin [27] found that the southern polar jet has shown a robust poleward shift, while the northern polar jet has shifted equatorward in most regions and seasons.Liu et al. showed in [24] that, in a simulation of the Last Glacial Maximum, NCAR's CCSM4 model indicates that, in the Southern Hemisphere, the ice line advances equatorward while the jet shifts poleward.In [14] Francis and Vavrus found evidence to support a linkage between rapid Arctic warming and more frequent high-amplitude, wavy jet-stream configurations (though they considered zonally asymmetric aspects of the flow which our model does not simulate), and in [22] Karamperidou, Cioffi, and Lall considered meridional surface temperature gradients and found them to be determinants of large-scale atmospheric circulation patterns.
The behavior of our model shares qualitative features with these investigations.An increase in radiative forcing, as from increased greenhouse gas concentrations, results in an initial poleward movement of the polar jet, followed by a equatorward shift of averaged locations and quasi-periodic oscillations, under greater forcings.Our results may also be compared to those of MS18 [29] and SAR18 [41], both of which used EBMs to demonstrate the influence of changing Hadley cell boundaries on the location of mid-latitude storm tracks.Our results do not contradict those findings but suggest that the latitudinal distribution of clouds may play a significant role as well.
More broadly, the cloud factor function in our model may be regarded as a prototype for further investigations.The cubic Hermite spline used to define the cloud factor function in this article depends only on a small number of fixed values, those at the equator, the Hadley cell boundary, the pole, and at the location of the maximum absolute value of the temperature gradients (see Subsection 2.1).But additional data points, including interactive data points in more elaborate models that incorporate physical processes influencing cloud cover at other latitudinal locations, might improve the climate sensitivity of the model considered here and add further insight into the dynamics of the polar jets.

Appendices A Cloud Factor Function Formula
The formula for the cloud factor function is shown here when the extratropical maximum occurs at θ latitude.

B Solution Methodology For The Initial Boundary Value Problem
The initial boundary value problem (IBVP) ( 16) falls in the class of linear evolution problems for which various numerical methods have been developed.We have employed in this paper an implicit finite difference method (FDM) based on the Crank-Nicholson scheme [1,47].This scheme has the desirable property of being inherently stable.More specifically, we subdivide the spatial variable interval [0,1] uniformly in I subintervals (x i , x i+1 ), i = 0, ..., I where x i = i∆x; ∆x being the spatial step size that is set to be 10 −3 (See Figure 10).Similarly, we consider for the time variable t, the equidistant sequence t n = n∆t; n = 0, 1, ..., N , where the time step ∆t is set to be 1 day and N is chosen large enough for the temperature to reach the asymptotic regime, i.e, the equilibrium of the solution of the IBVP (16).For the simplicity of the publication, we introduce the auxiliary variable T to denote either the temperature of the atmospheric layer, T a or the temperature of the surface layer, T s .We then approximate T (x i , t n ) by T n i where T n i is the solution of the algebraic system resulting from the adopted finite difference scheme.
The derivatives that occur in the IBVP ( 16) are approximated as follows.First, we have distributed the spatial derivative and then we have used the following second order approximation, and The first order time derivative is replaced by a second order approximation using the Crank-Nicholson relations [1,47] ∂T ∂t (x i , t n+ A schematic interpretation or cone of dependance of the adopted FDM discretization is depicted in Figure 10.It shows the implicit nature of this scheme.It also reveals that the evaluation of the temperature at the boundaries T n 0 (resp.T n I ) requires the values of T n −1 (resp.T n I+1 ).These "fictitious" values are set to be T n −1 = T n 0 and T n I+1 = T n I ; n = 0, ..., N .This choice results from the first order approximation of the boundary condition, IBVP (16).Note that the algebraic system ( 27) can be expressed in a compact representation as follows, Where A and B are block diagonal matrices whose entries are explicitly given in equations C.1 -C.14, pages 88 -92 in [36].The vector T consists of the temperature values for the atmosphere followed by the surface temperature values.The components of the vector b consist of all terms not linear in temperature.
The temperature gradients reported in Figures 6, 7, and 9 have been evaluated with the software package (numpy.gradient)[48].This routine computes the gradient using second order accurate central differences in the interior points and either first or second order accurate one-side differences at the boundaries.

Figure 1 :
Figure 1: Zonal mean cloud fraction from CMIP3 models and compared to observations (International Satellite Cloud Climatology Project, ISCCP).7

Figure 2 :
Figure 2: Cubic Hermite spline cloud factor C f plotted as a function of latitude from equator to pole with the first extratropical maximum θ = 50 • latitude.In general, the location of the the first extratropical maximum is interactive and varies in time.See equation 22 in Appendix A.

Figure 5 :
Figure 5: Initial atmospheric albedo α a plotted as a function of latitude from equator to pole for cubic Hermite spline cloud factor C f (θ) using α p0 = 0.25 + 0.38x 4 and the graph in Figure 2.

Figure 6 :
Figure 6: Meridional profiles of T s and T a with temperature gradients.Increased greenhouse gas concentrations are modeled by decreasing values of A out .The maximum value of |∂/∂θ(T a + T s )| (scaled by a factor of 0.15 for display purposes) represents the average latitude of the polar jet under the indicated forcings of A out .In plots (a) through (c), the jet moves poleward monotonically as A out decreases, but the gradient in plot (d) begins to form an approximate plateau with oscillatory equilibrium location.The graph shown in subfigure (d) is for a time step at which the maximum value of |∂/∂θ(T a + T s )| occurs at latitude 64.027.

Figure 7 :
Figure 7: Meridional profiles of T s and T a with temperature gradients for A out = 211 Wm −2 at two different time steps beyond 8725 days.A out = 211 Wm −2 is the largest integer value of A out at which oscillatory behavior of the maximum value of |∂/∂θ(T a + T s )| occurs.

Figure 9 :
Figure 9: Sample maxima of |∂(T a + T s )/∂θ| indicated with a red dot, along with displays of approximate plateaus of maximum values, at particular time steps for values of A out < 212 Wm −2 for which the jet is oscillatory.With each time step, the maximum on each plot shifts to a different latitude.

Figure 10 :
Figure 10: A schematic interpretation of the FDM approximation.The value of T n+1 i

Table 3 :
Model Data for High Forcings and Decreasing Jet Oscillations gave evidence that hothouse climate states may have different climate sensitivities per doubling of CO 2 than Earth's present state.In their study of early Paleogene and possible future high temperature modern climates, the temperature gain with each doubling of CO 2 was not constant according to their model, but instead increased with increasing CO 2 concentrations.By contrast, climate sensitivity of our model varies strongly and nonmonotonically with temperature and A out , encompassing values that are both above and below reasonable estimates of the modern Earth's climate sensitivity, but also Earthlike sensitivity at high temperatures and forcings.Unit increases in forcing, from A out = 208 to 207 and A out = 207 to 206 Wm −2 each result in an increase of the atmospheric temperature T a by only 0.1 • C.But unit increases from lower values of A out , corresponding to higher temperatures, shown in Table3result in increases of 0.9 and 0.8 • C (corresponding to climate sensitivities of 3.5 • C and 3.1 • C respectively, assuming the IPCC's reported effective radiative forcing of 3.9 Wm −2 for a doubling of CO 2 ).Data in Table2reveals an unrealistically high climate sensitivity for the larger consecutive values of A out compared to Earth's modern climate, and unrealistically low climate sensitivity for smaller consecutive values of A out .An increase in forcing from A out = 214 to 213 Wm −2 results in an increase of T a by 2.4 • C, but as A out decreases, the temperature increases decline.For consecutive large values of A out , it is likely that positive shortwave feedbacks, created by the interactive clouds and surface albedo, nearly cancel out the negative feedback associated with the increase of outgoing longwave radiation with increasing temperatures.We note, however, that over the full range of A out values, the average climate sensitivity is evidently closer to Earth-like climate sensitivity.Comparing the data for A out = 214 and A out = 202 Wm −2 , the ratio of temperature (T a ) IBVP(16)is then replaced by the following algebraic system,