Polar Jet Stream Fluctuations in an Energy Balance Model

TL;DR

This study uses a simplified climate model to analyze how increased greenhouse gases affect the position of the jet stream, revealing initial poleward shifts followed by equatorward movement at higher forcing levels.

Contribution

It introduces a novel cloud factor function in an energy balance model to simulate eddy-driven jet stream fluctuations under increased radiative forcing.

Findings

Initial poleward shift of jet stream with increased forcing

Jet stream location becomes quasi-periodic at higher forcing levels

Mean jet stream position moves equatorward with stronger forcing

Abstract

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.

Tables

Table 1. Table 1: Parameter values for the EBM as in Rose and Marshall [ 38 , 39 ] .

Parameter	Units	Numerical Value
a	m	$6.373\times {10}^6$
$S_0$	$W{m}^{-2}$	1367
$s_2$		-0.48
$C_a$	$J{m}^{-2}{}_{}{}^{\circ }C_{}^{-1}$	${10}^7$
$C_s$	$J{m}^{-2}{}_{}{}^{\circ }C_{}^{-1}$	${10}^7$
$K_a$	$m^2{s}^{-1}$	$2.7\times {10}^6$
$K_s$	$m^2{s}^{-1}$	$5.2\times {10}^5$
$B_{up}$	$W{m}^{-2}{}_{}{}^{\circ }C_{}^{-1}$	15
$A_{up}$	$W{m}^{-2}$	238
$B_{out}$	$W{m}^{-2}{}_{}{}^{\circ }C_{}^{-1}$	1.7
$A_{out}$	$W{m}^{-2}$	211

Table 2. Table 2: Average Locations of Max Temperature Gradients, Model Two

$A_{\text{o}ut}$	Mean Location of Unstable Jet	Standard Deviation
207 Wm-2	61.0∘ latitude	3.22
206 Wm-2	57.7∘ latitude	5.55
205 Wm-2	55.0∘ latitude	6.60
204 Wm-2	51.8∘ latitude	6.83

Table 3. Table 3: Model Two Latitudes of Polar Jet defined by max ⁡ | ∂ / ∂ θ ​ ( T a + T s ) | 𝜃 subscript 𝑇 𝑎 subscript 𝑇 𝑠 \max|\partial/\partial\theta(T_{a}+T_{s})| for t n ∈ [ 8725 , 8760 ] superscript 𝑡 𝑛 8725 8760 t^{n}\in[8725,8760] , the asymptotic regime of its solutions.

$A_{\text{o}ut}=207$	$A_{\text{o}ut}=206$	$A_{\text{o}ut}=205$	$A_{\text{o}ut}=204$
62.252	50.896	62.375	48.417
62.252	62.252	49.464	47.138
62.252	62.252	62.375	46.469
62.252	51.169	49.641	62.375
62.252	62.252	48.331	48.331
62.252	50.896	62.375	46.970
62.252	62.252	49.200	62.375
52.940	62.252	62.375	47.988
62.252	51.169	49.464	46.719
62.252	62.252	62.375	62.375
62.252	50.896	49.641	48.417
62.252	62.252	48.331	47.138
62.252	62.252	62.375	62.375
52.655	51.169	49.200	48.590
62.252	62.252	62.375	47.307
62.252	50.896	49.464	46.970
62.252	62.252	62.375	62.375
62.252	62.252	49.641	48.159
52.467	51.169	48.331	46.803
62.252	62.252	62.375	62.375
62.252	50.896	49.200	48.417
62.252	62.252	62.375	47.138
52.279	62.252	49.464	62.375
62.252	51.169	62.375	48.504
62.252	62.252	49.641	47.307
62.252	50.896	48.331	46.886
52.186	62.252	62.375	62.375
62.252	62.252	49.200	48.245
62.252	51.169	62.375	46.886
62.252	62.252	49.464	62.375
52.186	50.896	62.375	48.417
62.252	62.252	49.641	47.054
62.252	62.252	48.331	62.375
62.252	51.169	62.375	48.504
52.092	62.252	49.200	47.307
62.252	50.896	62.375	46.886

Equations

$$C_a\frac{\partial T_a}{\partial t}$$

$$C_s\frac{\partial T_s}{\partial t}$$

$$F_{ground}^{\downarrow }=\frac{(1-{\alpha }_g){\mathcal{T}}_{sw}}{(1-{\alpha }_a{\alpha }_g)}\frac{S_{0}s(x)}{4}$$

$$F_{atm}^{\downarrow }=(1-{\alpha }_g-{\mathcal{T}}_{sw})(1+\frac{{\alpha }_g{\mathcal{T}}_{sw}}{1-{\alpha }_a{\alpha }_g})\frac{S_{0}s(x)}{4}$$

$$F_{TOA}^{\uparrow }=({\alpha }_a+\frac{{\mathcal{T}}_{sw}^2{\alpha }_g}{(1-{\alpha }_a{\alpha }_g)})\frac{S_{0}s(x)}{4}$$

$$s(x)=1+s_2{P}_2(x),$$

$${\alpha }_p={\alpha }_{p,atm}+{\alpha }_{p,ground}={\alpha }_a+{\mathcal{T}}_e{\alpha }_g$$

$${\mathcal{T}}_e=\frac{{\mathcal{T}}_{sw}^2}{1-{\alpha }_a{\alpha }_g},$$

$${\alpha }_{p0}={\alpha }_{p0}^0+{\alpha }_{p0}^1{x}^4,$$

$${\alpha }_a={\mathcal{C}}_f({\alpha }_{p0}-{\alpha }_{clear})+{\alpha }_{clear},$$

$${\mathcal{T}}_{sw}=1-{\alpha }_a-A_{sw},$$

$${\alpha }_g=0.40-0.34\tanh (T_s+6.7).$$

$${\alpha }_p=\frac{1}{2}{\int }_{-1}^1{\alpha }_p(x)\text{s}(x)dx$$

$${\alpha }_a=\frac{1}{2}{\int }_{-1}^1{\alpha }_a(x)\text{s}(x)dx,$$

$${\overline{{\mathcal{T}}_e\alpha }}_g={\alpha }_p-{\alpha }_a$$

$$F_{out}=$$

$$F_{up}=$$

$$C_a\frac{\partial T_a}{\partial t}$$

$$-A_{out}-B_{out}{T}_a+\frac{C_a{K}_a}{a^2}\frac{\partial }{\partial x}\left[(1-x^2)\frac{\partial T_a}{\partial x}\right]$$

$$C_s\frac{\partial T_s}{\partial t}$$

$$+\frac{C_s{K}_s}{a^2}\frac{\partial }{\partial x}\left[(1-x^2)\frac{\partial T_s}{\partial x}\right]$$

$$\sqrt{1-x^2}\frac{\partial T_a}{\partial x}{\text{vline}}_{x=-1,0,1}=\sqrt{1-x^2}\frac{\partial T_s}{\partial x}{\text{vline}}_{x=-1,0,1}=0;t>0$$

$$\frac{\partial T_a(t,x)}{\partial x}.$$

$$\frac{1}{2}\frac{\partial }{\partial \theta }(T_a(t,\theta )+T_s(t,\theta )).$$

$$C_f(\theta )=\{\begin{array}{ll}0.000111111(\frac{30-\theta }{15}+1){\theta }^2+0.001(\frac{\theta }{15}+1)(30-\theta {)}^2& 0\le \theta \le 30\\ & \\ 0.000111111(1+\frac{\theta -30}{15})(60-\theta {)}^2+0.001(1+\frac{60-\theta }{15})(\theta -30{)}^2& 30\le \theta \le 60\\ & \\ 0.001(\frac{\theta -60}{15}+1)(90-\theta {)}^2+0.001(\frac{90-\theta }{15}+1)(\theta -60{)}^2& 60\le \theta \le 90\end{array}$$

$$\frac{\partial T}{\partial x}(x_i,t^n)\approx \frac{T_{i+1}^n-T_{i-1}^n}{2\Delta x}$$

$$\frac{{\partial }^{2}T}{\partial x^2}(x_i,t^n)\approx \frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}$$

$$\frac{\partial T}{\partial t}(x_i,t^{n+\frac{1}{2}})=\frac{1}{2}[\frac{\partial T}{\partial t}(x_i,t^{n+1})+\frac{\partial T}{\partial t}(x_i,t^n)]$$

$$\frac{\partial T}{\partial t}(x_i,t^{n+\frac{1}{2}})\approx \frac{T_i^{n+1}-T_i^n}{\Delta t}$$

$$\beta (T_{a_i}^{n+1}-T_{a_i}^n)$$

$$-(B_{up}+B_{out})T_{a_i}^{n+1}+B_{up}{T}_{s_i}^{n+1}+A_{up}-A_{out}$$

$$+(1-{\alpha }_{a_i}^n-{\mathcal{T}}_{s{w}_i}^n)(1+\frac{{\alpha }_{g_i}^n{\mathcal{T}}_{s{w}_i}^n}{(1-{\alpha }_{a_i}^n{\alpha }_{g_i}^n)})\frac{S_{0}s(x_i)}{4}$$

$$+{\beta }_i^{{}^{\prime }}(T_{a_{i+1}}^n-2T_{a_i}^n+T_{a_{i-1}}^n)-{\beta }_i^{{}^{\prime \prime }}(T_{a_{i+1}}^n-T_{a_{i-1}}^n)$$

$$-(B_{up}+B_{out})T_{a_i}^n+B_{up}{T}_{s_i}^n+A_{up}-A_{out}$$

$$+(1-{\alpha }_{a_i}^n-{\mathcal{T}}_{s{w}_i}^n)(1+\frac{{\alpha }_{g_i}^n{\mathcal{T}}_{s{w}_i}^n}{(1-{\alpha }_{a_i}^n{\alpha }_{g_i}^n)})\frac{S_{0}s(x_i)}{4}]$$

$$\gamma (T_{s_i}^{n+1}-T_{s_i}^n)$$

Full text

**Polar Jet Stream Fluctuations in an Energy

Balance Model

**

Cord Perillo111Department of Mathematics and Interdisciplinary Research Institute for the Sciences, California State University, Northridge, Northridge, CA 91330-8313. Email: cord.perillo@csun.edu., David Klein222Department of Mathematics and Interdisciplinary Research Institute for the Sciences, California State University, Northridge, Northridge, CA 91330-8313. Email: david.klein@csun.edu. Orcid ID: https://orcid.org/0000-0003-1964-4378, Rabia Djellouli 333Department of Mathematics and Interdisciplinary Research Institute for the Sciences, California State University, Northridge, Northridge, CA 91330-8313. Email: rabia.djellouli@csun.edu.

**Abstract. We investigate the effect of increased greenhouse gas concentrations on the zonally averaged location of the jet stream boundary of the polar cell 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 cell structure of Earth’s global atmospheric circulation. A simulation of increased greenhouse gas concentrations initially results in a poleward movement of the polar jet stream’s mean latitude, but as the concentration increases the jet moves equatorward. When the model takes both the atmospheric and the surface temperature gradients into account, under high concentrations the location of the jet stream becomes unstable and quasi-periodic as its mean location moves equatorward.

**

INDEX TERMS: energy balance model, polar jet fluctuations, temperature gradient, quasi-periodicity

1 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 [7], Sellers [41], 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, [11, 18, 20, 31, 24, 45, 42, 5] 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., [4, 13, 17, 22, 27, 28, 50, 51]. 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 [30] (hereafter MS18) and and Siler, Roe, Armour [40] (hereafter SRA18).

MS18 [30] 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 $\gamma$ 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 $\gamma$. 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 [40] 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 [40] found that down-gradient 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 cell structure of Earth’s global atmospheric circulation. More specifically, the cloud factor function, $\mathcal{C}_{f}(\theta)$, is a dimensionless quantity that represents the fraction of the zonally averaged planetary albedo at latitude $\theta$ attributable to clouds. $\mathcal{C}_{f}(\theta)$ is the mechanism by which our model takes into account the influence of Earth’s rotation on the atmospheric general circulation. 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 jet stream and suggest that the most stable regions of the jet stream are located where the magnitude of the temperature gradient is maximized, the same proxy utilized as in MS18 [30]. Similar to MS18 [30], we interpret the latitude where this occurs as the averaged southern boundary of the Northern Hemispheric polar cell, 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 the polar cell as it moves dynamically with each time step until the system reaches equilibrium.

We must point out that the Hadley cell edge is not interactive. 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 polar jet stream fluctuations, which has the potential to be incorporated into more complex climate models that include Hadley cell dynamics.

The cloud factor function, $\mathcal{C}_{f}(\theta)$ is constructed so that its minimum and maximum values can be interpreted as the primary cell boundaries of the global circulation. $\mathcal{C}_{f}(\theta)$, is a cubic Hermite spline continuously differentiable at all latitudes in a hemisphere with specified values at the equator, the Hadley cell boundary, the location of the polar jet stream, and the pole. We model it in two different ways. In the first case, called “Model One”, the location of the polar jet is determined by a gradient of the atmospheric temperature. In the second case, called “Model Two”, the location of the jet is determined by a gradient of the average of the atmospheric and surface temperatures. The two versions exhibit the same general qualitative behavior, but Model Two also exhibits instabilities and quasi-periodic behavior with interesting and suggestive physical interpretations (see Subsection 3.2).

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 polar cell boundary 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 for Model Two.

2 Model Description

Our EBM consists of an ocean covered surface layer and an overlying atmospheric layer. Throughout, we let $x=\sin\phi$, where $\phi$ is latitude444This formula assumes that $\phi$ is measured in radians. Later, in the context of temperature gradients, it will be calculated as $x=\sin(\pi\theta/180)$ where $\theta$ is given in degrees., so that $-1\leq x\leq 1$, but because our aqua-planet is symmetrical, we will generally display data only for the northern hemisphere, $0\leq x\leq 1$.

Let $T_{s}$ and $T_{a}$ represent the zonally averaged temperatures of the surface and atmosphere respectively. 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)555In particular we will vary the parameter $A_{\text{out}}$ in Eq. (16a) to simulate changes in greenhouse gas concentrations. The time evolution of the temperatures are solutions to coupled differential equations of the form,

[TABLE]

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 heat flux from the ocean to the atmosphere from all sources. 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}^{\downarrow}$ and $F_{ground}^{\downarrow}$, in Eqs. (1a) and (1b) represent incoming solar radiation flux and both depend on the atmospheric albedo, $\alpha_{a}$, and ground albedo, $\alpha_{g}$.

To model the dependence of $F_{atm}^{\downarrow}$ and $F_{ground}^{\downarrow}$ on $\alpha_{a}$, and $\alpha_{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 $\alpha_{a}$, transmission by a factor $\mathcal{T}_{sw}$ (the transmissivity of shortwave radiation), and absorption by a factor $A_{sw}=1-\alpha_{a}-\mathcal{T}_{sw}$.

Summing up the infinite number of transmissions and reflections between the atmosphere and the ground, the total downwelling flux to the ground $F_{ground}^{\downarrow}$, the net flux into the atmosphere $F_{atm}^{\downarrow}$, and the total upwelling flux at the top of the atmosphere (TOA), $F_{TOA}^{\uparrow}$ are given by,

[TABLE]

[TABLE]

[TABLE]

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,

[TABLE]

with $s_{2}=-0.48$.

From the last expression, the planetary albedo is identified as,

[TABLE]

where

[TABLE]

and $\mathcal{T}_{e}\alpha_{g}$ can be considered as the contribution from the ground albedo $\alpha_{g}$ to the planetary albedo modulated by the interactions with the atmosphere.

2.1 Cloud Factor Function

In order to assign latitudinal values to the ground and atmospheric albedos, $\alpha_{g}$ and $\alpha_{a}$, we first introduce a cloud factor function, $\mathcal{C}_{f}=\mathcal{C}_{f}(\theta)$, a dynamic function of latitude $\theta=\theta(x)$. The function $\mathcal{C}_{f}$ represents the fraction of the zonally averaged planetary albedo at latitude $\theta$ attributable to clouds (see Eq.(9) below). It is the mechanism by which our model takes into account the influence of planetary rotation on the accumulation of clouds by the global atmospheric circuclation.

Consistent with earlier studies, Liu et al. [26] found that the planetary rotation has a strong effect on the cell structure of the atmosphere and, in particular, the Hadley cell becomes narrower and weaker with a higher planetary rotation rate. Liu et al. [26] found that the magnitudes of cloud adjustments with rotation are comparable to, or exceed, the magnitude of the climatological ocean heat uptake, and concluded that in energetic terms cloud variations are arguably more important than ocean transport variations in meridional heat transport.

It is difficult to measure cloud cover in the polar regions due to a number of factors: thin and low lying clouds, low visibility between clouds and the underlying surface, and polar conditions create an unusual amount of near surface hazes and fogs [8]. Because of these problems, there is an uncertainty in cloud cover over the polar regions. Vavrus et al. [49] 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\%$. Both conclude that the average polar cloud fraction is increasing as the sea ice extent has been decreasing.

Taking these findings into consideration, we construct $\mathcal{C}_{f}(\theta)$ through the use of cubic Hermite splines.777Cubic Hermite splines are continuously differentiable at all points, including juncture points. We note that the use of linear splines instead of cubic splines results in qualititative similar final results. The 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 boundaries of an initial idealized cell structure of the general circulation. Specifically, the coordinates are $(0,0.9),(30,0.1),(50,0.9)$ and $(90,0.9)$ so as to represent high cloudiness at the equator as well as from 50∘ degrees and poleward, and low cloudiness at 30∘ degrees, but, as we explain in Sect. 2.5, the graph will change with the time steps in the numerical runs of our models. A sample graph is shown in Figure 2.

As pointed out in [29], the Southern Hemisphere polar jet is stable at $50^{\circ}$ 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 [40] 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^{\circ}$ latitude (see Figures 2f and 3b in [40]).

2.2 Albedo Functions

Our modeling of the atmospheric albedo $\alpha_{a}$ and the ground albedo $\alpha_{g}$ begins with an initial approximate estimate of the planetary albedo. As a reference frame and a guide, Figure 8 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 $\alpha_{p0}$ of the form

[TABLE]

The coefficients $\alpha_{p0}^{0}$ and $\alpha_{p0}^{1}$ are chosen along with parameters for the ground albedo in Eq (11) so that the equilibrium average planetary albedo given in Section 3 approximates Earth’s average planetary albedo, and in order to specify initial equilibrium locations of maximal absolute values of temperature gradients. Figure 4 shows a plot of $\alpha_{p0}$ for this choice of parameters: $\alpha_{p0}^{0}=0.25$ and $\alpha_{p0}^{1}=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 models, which instead will vary in time in a way that depends on the global temperature distribution. We use $\alpha_{p0}$, along with the cloud factor function $\mathcal{C}_{f}$, to define the atmospheric contribution to the planetary albedo as:

[TABLE]

where $\alpha_{clear}$ is the clear sky (cloud free) albedo of the atmosphere which we take as constant, $\alpha_{clear}=0.149$ [43]. An initial sample plot of the atmospheric albedo is given in Figure 5. We note that plots of $\alpha_{a}$ depend on $\mathcal{C}_{f}$, and in the sequel $\mathcal{C}_{f}$ will in turn depend on temperature gradients, as described below, so the atmospheric albedo will also depend on temperature gradients.

We can now define the atmospheric transmittance of short wave radiation (SWR) in terms of $\alpha_{a}$ as

[TABLE]

where $A_{sw}=0.05$ is the atmospheric absorption of SWR [20]. We note that $\mathcal{T}_{sw}$ depends on $\mathcal{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,

[TABLE]

2.3 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 restrict our atmospheric and modulated ground albedos to these approximate values. They are of necessity approximate because the atmospheric and ground albedo contributions in our model are dynamic and therefore fluctuate.

The total planetary albedo $\alpha_{p}$ is given by,

[TABLE]

where, as before, $x$ is the sine of latitude, $\alpha_{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 $\alpha_{a}$ is defined as,

[TABLE]

where $\alpha_{a}(x)$ is the zonally averaged atmospheric albedo at $x$. Therefore, we define the total planetary effective ground albedo by,

[TABLE]

2.4 The Model

The energy balance equations in this section will be incorporated into what we describe as Model One and Model Two in Section 3, but with different dynamics of the cloud factor function. The algorithms for the dynamics are given in Subsection 2.5, and they depend on the equations for the climate model developed in this section.

We begin by linearizing the terms $F_{up},F_{out}$ in Eqs (1a) and (1b) and write,

[TABLE]

Then, 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]

Table 1 lists the parameter values appearing in Eqs. (16a) and (16b). These are the same choices made by Rose and Marshall [38, 39].

The initial ($t=0$) temperature profile is specified below, and the dynamic feature of the cloud factor function is explained in the following subsection.

The system of equations given in Sect.2.4 is defined for $-1\leq x\leq 1$, where $x<0$ is the Southern Hemisphere and $x>0$ is the Northern Hemisphere. The Southern and Northern hemispheres are symmetric with respect to the boundary value problem. So, by symmetry, we need only consider the solution from $0\leq x\leq 1$.

2.5 Cloud Function Dynamics and Polar Jet Stream

The response of eddy-driven jets and storm track latitudes to arctic amplification and changing meridional temperature gradients has been analyzed extensively (e.g., [4, 3, 13, 28, 30, 40, 50] and references therein). With the thermal wind equations in mind, we identify the mean latitudinal position of the polar jet stream, at any time $t$, with the location of the maximum value of a meridional temperature gradient in two different ways, which we refer to as Model One and Model Two.

Model One - We identify the location of the polar jet at any time t as the maximum absolute value of

[TABLE]

This derivative is proportional to the rate of change of mid atmospheric temperature, $T_{a}$, with respect to the Euclidean distance, from the latitude with coordinate $x$, to the plane of the equator. 2. 2.

Model Two - In this case we identify the jet location with the maximum of the absolute value of the gradient of the average of the surface and atmospheric temperatures, i.e., as the absolute value of

[TABLE]

Other choices for temperature gradient, involving the atmospheric temperature, $T_{a}$, are possible. But unlike those of Models One and Two, the alternatives have been investigated and found to lead to less physical behavior of the cloud factor function in our model.

For each of the two models, the temperature gradients (17) and (18) are coupled in turn 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). At each time step $t^{n}=n\Delta t$ of the numerical calculation, the cloud factor function $\mathcal{C}_{f}$ is updated so that it takes the value $0.9$ for all $x$ greater than or equal to the value of $x$ which maximizes the absolute value of the meridional temperature gradient, in the previous time step, for the model in question. For example, the graph in Figure 2 corresponds to a maximum meridional temperature gradient occurring at $60^{\circ}$ latitude. Since the atmospheric albedo $\alpha_{a}$ depends on $\mathcal{C}_{f}$ (c.f. Eq.(9)), it is updated in this way. Similarly, the ground albedo $\alpha_{g}$ (which is a function of latitude) is updated at each time step according to the value of the surface temperature $T_{s}$ in the previous time step (see Eq.(11)). Numerical approximation details are described in Appendix B.

3 Numerical Results

In this section we present numerical results for our Model One and Model Two in response to increases in radiative forcing, such as from increased greenhouse gas concentrations. Following [38], to simulate increased greenhouse gas concentrations, 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 gradients, in equilibrium, given by Eqs. (17) and (18), 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 models are non autonomous, equilibrium temperature and temperature gradient values for each experiment must be found by numerically running them out to equilibrium 999For low values of $A_{\text{o}ut}$, Model Two does not reach equilibrium with a constant location of the temperature gradient. Instead the maximum temperature gradient becomes quasi-periodic, oscillating between different latitudes, as elaborated below. 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.

3.1 Model One

To set a reference climate for Model One, we take $A_{out}$ = 213.9 Wm*-2*. In equilibrium, this yields a climate with a high planetary albedo, $\alpha_{p}=0.337$, and cold average temperatures given by $T_{s}=5.9^{\circ}$C and $T_{a}=8.0^{\circ}$C. The maximum absolute value of the atmospheric temperature gradient occurs at $51.7^{\circ}$ latitude. This is our proxy for the average latitude of the jet stream. The temperature and gradient distributions are displayed in Figure 6(a).

By decreasing the parameter $A_{out}$, we introduce a radiative forcing to simulate an increase in greenhouse gas concentrations in the model. With a decrease to 213.5 Wm*-2* of $A_{out}$, the average jet location begins to move poleward and at equilibrium reaches $58^{\circ}$ latitude with global average temperatures of $T_{s}=9.4^{\circ}$C and $T_{a}=11.1^{\circ}$C and planetary albedo of .327. The temperature and gradient distributions are presented in Figure 6(b).

With $A_{out}=213$ Wm*-2* the global average temperatures and planetary albedo more closely resemble those of present day Earth, with $T_{s}=14.4^{\circ}$C, $T_{a}=15.6^{\circ}$C and a planetary albedo of .300. The boundary of the polar cell increases to $62.1^{\circ}$ latitude as displayed in Figure 6(c). As $A_{out}$ decreases to 206 Wm*-2*, the global average surface and atmospheric temperatures increase as does the boundary of the polar cell which reaches a maximum latitude of 80.0 ∘, depicted in Figure 6(d).

However, as $A_{out}$ decreases further, to values less than 206 Wm*-2*, the movement of the polar jet location is reversed. It begins to decrease. When $A_{out}=205$ Wm*-2*, the equilibrium location of the jet decreases to $78.5^{\circ}$ latitude, as shown in Figure Figure 6(e). For $A_{out}=203$ Wm*-2*, the jet location moves further equatorward to $77.7^{\circ}$ latitude. See Figure 6(f).

3.2 Model Two

Model Two exhibits a behavior not shared with Model One. With sufficient forcing (i.e., low values of $A_{out}$), the modulus of the temperature gradient given by Eq.(18) 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 the emergence of instability of the jet stream.

To explain this phenomenon, we start with $A_{out}$ = 211 Wm*-2*, and then decrease it by integer values. The choice of $A_{out}$ = 211 Wm*-2* results in an equilibrium climate with planetary albedo, $\alpha_{p}=0.324$, and global average temperatures given by $T_{s}=14.9^{\circ}$C and $T_{a}=12.2^{\circ}$C. The maximum absolute value of the atmospheric temperature gradient occurs at $49.1^{\circ}$ latitude, the proxy for the average latitude of the jet stream in this model. The temperature and gradient distributions are displayed in Figure 7(a). Thereafter, the location of the maximum modulus of the temperature gradient increases monotonically as $A_{out}$ decreases until $A_{out}$ = 208 Wm*-2* for which the average jet stream location is $60.6^{\circ}$ latitude, shown in Figure 7(c). At $A_{out}=207$, Figure 7(d) shows the beginning of the formation of an approximate plateau of maximum values of the modulus of the temperature gradient.

Numerical calculations, displayed in Appendix C, show that the latitudinal locations of the maximum modulus of the temperature gradient are eventually quasi-periodic, oscillating among a finite but increasing number of locations as $A_{out}$ is decreased. The means and standard deviations of the latitudes of maximal moduli of the temperature gradients are displayed in Table 2. Figures 8(a) through 8(d) show portions of the weighted temperature gradient plots close to the maximum values.

4 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 [4] found that all jets migrate poleward in the twenty-first century. Using reanalysis, Manney and Hegglin [28] 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 [25] 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 [13] Francis and Vavrus found evidence to support a linkage between rapid Arctic warming and more frequent high-amplitude, wavy jet-stream configurations, 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 models share qualitative features with these investigations. A simulation of increased greenhouse gas concentrations results in an initial poleward movement of the polar jet, followed by a equatorward shift under greater forcings. In Model Two, the location of jet, under high concentrations, exhibits quasi-periodic behavior and instability. Our results may also be compared to those of MS18 [30] and SAR18 [40], 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.

Appendix A Cloud Factor Function Formulas

The formula for a cloud factor function is shown here. The Hadley cell boundary is taken as $30^{\circ}$ latitude and the first extratropical maximum occurs at $60^{\circ}$ latitude. The formula for the graph displayed in Figure 2 is given by

[TABLE]

Appendix B Solution Methodology For The Initial Boundary Value Problem

The initial boundary value problem (IBVP) (16c) 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 [2, 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\Delta x$; $\Delta x$ being the spatial step size that is set to be $10^{-3}$ (See Figure 9). Similarly, we consider for the time variable t, the equidistant sequence $t^{n}=n\Delta t$; $n=0,1,...,N$, where the time step $\Delta t$ is set to be 1 and N is chosen large enough for the temperature to reach the asymptotic regime, i.e, the equilibrium of the solution of the IBVP(16c). 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_{i}^{n}$ where $T_{i}^{n}$ is the solution of the algebraic system resulting from the adopted finite difference scheme.

The derivatives that occur in the IBVP (16c) are approximated as follows. First, we have distributed the spatial derivative and then we have used the following second order approximation,

[TABLE]

and

[TABLE]

The first order time derivative is replaced by a second order approximation using the Crank-Nicholson relations [2, 47]

[TABLE]

and

[TABLE]

sequentially, IBVP(16c) is then replaced by the following algebraic system,

[TABLE]

where

[TABLE]

A schematic interpretation or cone of dependance of the adopted FDM discretization is depicted in Figure 9. It shows the implicit nature of this scheme. It also reveals that the evaluation of the temperature at the boundaries $T_{0}^{n}$(resp. $T_{I}^{n}$) requires the values of $T_{-1}^{n}$(resp. $T_{I+1}^{n}$). These “fictitious” values are set to be $T_{-1}^{n}=T_{0}^{n}$ and $T_{I+1}^{n}=T_{I}^{n};\;n=0,...,N$. This choice results from the first order approximation of the boundary condition, IBVP (16c).

Note that the algebraic system (24b) can be expressed in a compact representation as follows,

[TABLE]

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 linear system (26) is solved using LPACK package (routine -gesv)[1] that is based on LU type decomposition [14].

The temperature gradients reported in Figures 6 - 8 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.

Appendix C Quasi-Periodic Oscillations of Polar Jet in Model Two for Large Forcings