Planetary geostrophic Boussinesq dynamics: Barotropic flow, baroclinic instability and forced stationary waves

Motions on planetary spatial scales in the atmosphere are governed by the planetary geostrophic equations. However, little attention has been paid to the interaction between the baroclinic and barotropic flows within the planetary geostrophic scaling. This is the focus of the present study, which utilizes planetary geostrophic equations for a Boussinesq fluid supplemented by a novel evolution equation for the barotropic flow. The latter is affected by meridional momentum flux due to baroclinic flow and drag by the surface wind. The barotropic wind, on the other hand, affects the baroclinic flow through buoyancy advection. Via a relaxation towards a prescribed buoyancy profile the model produces realistic major features of the zonally symmetric wind and temperature fields. We show that there is considerable cancellation between the barotropic and the baroclinic surface zonal mean zonal winds. Linear and nonlinear model responses to steady diabatic zonally asymmetric forcing are investigated, and the arising stationary waves are interpreted in terms of analytical solutions. We also study the problem of baroclinic instability on the sphere within the present model.


INTRODUCTION
Using scale considerations, Burger (1958) suggested that for atmospheric motions on planetary scales, that is, scales comparable with the radius of the Earth, the vorticity is quasi-stationary and the vorticity equation takes the form of a balance between the divergence of the horizontal wind and the advection of planetary vorticity.Later Phillips (1963) proposed for the description of planetary scale dynamics the planetary geostrophic equations (PGEs), or geostrophic motions of type 2. In the PGEs the pressure is hydrostatically balanced and the horizontal wind is in geostrophic balance where the full variations of the Coriolis parameter f are considered.The vertical velocity in the anelastic approximation of the PGEs results solely from variations of f and there is only one prognostic equation, namely for temperature.Because of their reduced complexity the PGEs are part of the atmospheric module in some Earth system models of intermediate complexity (e.g., Petoukhov et al., 2000;Totz et al., 2018), allowing numerically efficient long-term climate simulations (for examples see Ganopolski and Rahmstorf, 2001;Claussen et al., 2002;Petoukhov et al., 2005).
Only recently has the range of validity of the PGEs been revised using currently available reanalysis data (Egger and Hoinka, 2017) and simplified general circulation model simulations (Dolaptchiev and Klein, 2013).The latter authors found, from investigating spectrally decomposed fields, that the horizontal fluxes of relative and planetary vorticity are almost divergence-free on the planetary scale.Egger and Hoinka (2017) showed that the vertical velocity from the PGEs captures the stationary features in the tropospheric zonal perturbations.However, the standard deviation of the vertical velocity was greatly underestimated due to the lack of synoptic scale dynamics in the PGEs.The important effect of the synoptic eddies has been incorporated into the PGEs using multiple scale asymptotics (Pedlosky, 1984;Dolaptchiev and Klein, 2013;Boljka and Shepherd, 2018) and statistical-dynamical approaches (Petoukhov et al., 2003;Coumou et al., 2011;Totz et al., 2018).
Despite the popularity of the PGEs, little attention has been paid to the evolution of the barotropic flow under the planetary geostrophic scaling.As stated by Bresch et al. (2006), the PGEs do not represent a closed set of equations, and an additional evolution equation for the barotropic pressure has thus been proposed to close the system.Using asymptotic expansion, Dolaptchiev and Klein (2009) have generalized the closure to the case of fully compressible flow with variable Coriolis parameter.The derived closure has the form of a prognostic equation for the barotropic vorticity and is a dynamical alternative to other diagnostic closures (e.g., Petoukhov et al., 2000).This study is a first attempt to address the effect of the closure on the planetary geostrophic dynamics by utilizing numerical simulations of a Boussinesq fluid on the sphere.
In addition, in the present study the linear and nonlinear responses of the PGE model to steady diabatic forcing is considered.The arising stationary waves are interpreted in terms of analytical solutions.We also study the problem of baroclinic instability within the PGEs.This was first done by Wiin-Nielsen (1961), but to our knowledge the problem on the sphere has not yet been considered.The latter is in contrast to quasi-geostrophic or primitive equation dynamics, where a large body of theoretical work already exists on the topic (e.g., Hollingsworth, 1975;Simmons and Hoskins, 1976;Baines and Frederiksen, 1978).
This paper is organized as follows.In section 2 an asymptotic derivation of the PGEs and the equation for the barotropic dynamics is presented.The representation of diabatic and frictional effects, as well as a summary of the nonlinear and linear model equations, can be found in section 3. The nonlinear model simulations are discussed in section 4 for different model configurations.In section 5 analytical wave solutions are presented and compared with the linear/nonlinear numerical simulations, and the problem of baroclinic instability is also studied.Concluding discussions can be found in section 6.

ASYMPTOTIC DERIVATION
Using asymptotic analysis the PGEs were derived by Dolaptchiev and Klein (2009) from the full compressible fluid flow equations.Here for the first time the evolution equation for the barotropic pressure is studied within the PGEs.In order to simplify the analysis, we consider the hydrostatic Boussinesq equations as a starting point.These equations are isomorphic to the primitive equations in pressure coordinates (Vallis, 2006).Although the Boussinesq approximation is limited to vertical scales smaller than the scale height of the atmosphere (see, for example, the recent work by Egger and Hoinka, 2018 on the validity of the incompressibility assumption), the Boussinesq equations are widely used for studying the large-scale circulation (e.g., Held and Hou, 1980;Vallis, 2006).The non-dimensional governing equations for a Boussinesq fluid on the sphere are where − →  (, ) denotes the horizontal velocity vector,  the vertical velocity, Φ the pressure fluctuations divided by the reference density and b the buoyancy, and where with the Coriolis parameter f and the radial unit vector − →   .
The source terms due to friction and diabatic effects are denoted by − →   and S b , respectively.The vertical coordinate is indicated by z and in the following  is longitude and  is latitude.The horizontal Nabla operator is defined as and for the horizontal divergence of − → , where a is the radius of the Earth.For the non-dimensionalization of the variables a reference horizontal velocity of U = 10 m ⋅ s −1 and a planetary horizontal length scale of L = 10 7 m is used.With the above scales the Rossby number  = U/f 0 L, where f 0 denotes the Coriolis parameter at 45 • N, takes a value of about 10 −2 .The reference value for the vertical velocity W is set to the standard value W = HU/L, where H = 10 4 m is the scale height of the atmosphere.The normalized pressure Φ and the buoyancy are non-dimensionalized using P = LUf 0 and B = P/H, respectively, in order to assure geostrophic and hydrostatic balance to leading order.The potential temperature  can be computed from buoyancy using b = g(− 0 )/ 0 , where g is gravity acceleration and  0 is a constant reference potential temperature.Note that the small parameter  and the characteristic scales used in this article differ from those described in Dolaptchiev andKlein (2009, 2013).In the latter studies a unified asymptotic approach is utilized, where the characteristic quantities for non-dimensionalization are valid for a variety of flow regimes.To keep the present asymptotic analysis concise, we begin here with the characteristic scales described above, which are appropriate for planetary scale motions.
We assume that each dependent variable from Equations 1-4 can be represented as an asymptotic series in terms of : where  = (, , , , Φ,   , − →   ).Substituting the ansatz above in Equations 1-4, we obtain as leading-order non-trivial asymptotic equations the PGEs where the zero superscript in all dependent variables has been dropped.Next, the flow is separated into barotropic and baroclinic parts.For instance, we write for the horizontal wind where the baroclinic part is indicated by a prime and the barotropic one is defined as with z a denoting the height of the atmosphere.Averaging Equation 8 vertically and applying rigid-lid boundary conditions, one obtains By taking the curl of Equation 6, one obtains the vanishing divergence of the planetary vorticity flux which for the barotropic component reads From Equations 12 and 14 it follows that ⟨v⟩ z vanishes and ⟨u⟩ z does not depend on longitude: From Equations 6, 8 and 10 the baroclinic wind satisfies where in the last equation we have used w = 0 at z = 0.For a given buoyancy field the baroclinic part of Φ can be found by integrating the hydrostatic balance (Equation 7): Thus, Equation 9gives a prediction of b from which − →  ′ and  can be determined.However, in order to determine the evolution of the vertically averaged zonal wind field ⟨u⟩ z , we need to consider the next-order asymptotic equations.From Equations 1 and 3 we obtain for the (1) zonal momentum equation and the () continuity equation Averaging Equation 21over  and z and applying vanishing vertical velocity at the top and at the bottom yields Here the vertical and zonal mean of  (1) is defined as 20), averaging the result zonally and vertically and making use of Equation 22, we derive a vorticity equation for the barotropic component of the flow where the vorticity is defined as  = − →   ⋅  × − →  and the source term on the right-hand side is given by .Expressing all nonlinear terms on the left-hand side of Equation 24 in terms of the meridional momentum flux yields (25) From the last equation one can determine the evolution of ⟨u⟩ z ; due to Equation 12 one can introduce a stream function In Equation 27, Δ =  2 denotes the horizontal Laplace operator in spherical coordinates and Ψ does not depend on  due to Equation 15.
Using asymptotic analysis, Boljka and Shepherd (2018) showed that there is a connection between the planetary scale barotropic flow equation and the preservation of angular momentum.An alternative to Equation 25 can be derived by multiplying Equation 20 with acos and averaging again zonally and vertically in order to obtain an equation for the vertical mean of the axial angular momentum M = acosu: where   =  cos  (0)  .In the last equation the angular momentum M lacks the planetary component (Ωa 2 cos 2 ), since the zonally and vertically averaged transport of planetary angular momentum by v (1) from Equation 20vanishes in any case due to Equation 22.

MODEL DESCRIPTION
In this section we discuss the parameterization of diabatic and frictional effects, and for this purpose the redimensionalized variables are considered.

Diabatic processes
The diabatic processes described by the term S b in Equation 9are modelled with the simple relaxation ansatz and diffusion where  is a relaxation time-scale and  b is a diffusion constant.The prescribed buoyancy profile b eq is separated into zonally symmetric part and deviations from it: Here ⟨b eq ⟩  accounts for the meridional temperature differences in a radiative equilibrium atmosphere.We utilize a relaxation profile of Held and Hou (1980) form, where    (, ) denotes the associated Legendre polynomial, corresponding to the zonal wavenumber m, and the total wavenumber n.The constants  h = 100 K and  v = 40 K are measure for the meridional and vertical temperature gradients, respectively.1To the buoyancy profile ⟨b eq ⟩  corresponds a zonal jet in thermal wind balance with the form if zero surface wind ⟨u eq (z = 0)⟩  is assumed.
The zonally asymmetric part  *  from Equation 30 models the differential heating of the atmosphere due to the land-sea thermal contrast.Here we choose an idealized representation of this effect by prescribing a buoyancy anomaly which is the sum of two spherical harmonics with zonal wave number 2: The magnitude of the zonally asymmetric perturbation is set to  p = 5 K.  *  has its maximum at the surface at around 50 • latitude and decays in the vertical with an exponential decay length scale  −1 = 1 km.The relaxation profile from Equations 31 and 33 is shown in Figure 1.Note that this profile is characterized by super-rotation at the Equator.

Frictional effects
The frictional effects on the baroclinic flow are incorporated by including turbulent eddy diffusion in the momentum equation: where K is an eddy diffusivity constant.Without frictional effects the baroclinic wind from Equation 17 will diverge at the Equator, if the gradient of Φ ′ does not vanish faster than the Coriolis parameter f as  → 0. Thus, the inclusion of diffusion acts as a form of regularization for the PGEs.We have to stress, however, that those frictional effects are not accounted for in the present asymptotic derivation.However, when considering the equatorial region with small f , it is natural to  33expect that higher-order effects (such as eddy dissipation) will modify the geostrophic balance.Such effects are often taken into account in studies on the tropical circulation by adding Rayleigh drag to the geostrophic balance (Matsuno, 1966;Gill, 1980) or diffusion (Schneider and Lindzen, 1977).The value of K used in our model was set uniformly to 5 m 2 ⋅s −1 , a value used in other idealized studies of the atmospheric circulation (e.g., Held and Hou, 1980).
The baroclinic wind stress at the surface is parameterized using the drag coefficient C D : With the inclusion of vertical diffusion in Equation 34 even if In Equation 36, Δz s is the thickness of the upper layer.We refer to this layer as the stratosphere (note, however, that in Equation 31 no separate assumptions on the stratification within this layer are made) and z t marks the troposphere height.Similar boundary conditions are used in other PGE-based models; see, for example, Petoukhov et al. (2000) and eqs.15 and 22 therein.Equation 36 is a considerable limitation of the dynamics at z s , but as discussed in section 4.2 this has no pronounced effect on the major features of the tropospheric dynamics.
Consistent with the eddy diffusion closure (Equations 34 and 35), the frictional effects in the vorticity source term ⟨S  ⟩ z,  from Equation 25 are represented by Ekman friction: where (0) denotes  at the lowest model level.

Summary of the model equations and numerical implementation
The dimensional governing equations of the planetary geostrophic model take the form of two prognostic equations for the buoyancy and barotropic vorticity, and diagnostic relations for Φ ′ and the baroclinic/barotropic wind components In Equations 39-42 the primed variables indicate deviations from the vertical, as defined in Equations 10 and 11.We use (unless otherwise stated) horizontal spectral discretization with a triangular truncation of T21 and five equidistant layers in the vertical with layer thicknesses of Δz = 2 km.The variables b,  and − →  are defined at the centres of the layers with  at the interfaces.Centred differences are used for the  1.We integrate the model for 10 4 days.After 300 days a steady state is reached and we use data from this state only for the analysis.

Linear model
We linearize Equation 38 around a zonally symmetric stationary basic state given by and obtain (47) where an asterisk denotes the perturbations from the basic state.In Equation 47 we neglected the tendency of the basic state due to frictional and diabatic source terms.Note that without source terms the zonally symmetric basic state from Equations 44-46 is a stationary solution of the model equations.If source terms are omitted in the linearization of the barotropic vorticity equation, the tendency of ⟨ * ⟩ z,  disappears as well due to Equation 46.In this case Equation 47 describes the complete linear dynamics of the system.We look for normal mode solutions of the horizontally and vertically discretized version of Equation 47.All fields are represented as a truncated series of spherical harmonics, so for b* at vertical level z l this series reads with T denoting the spectral truncation,  the frequency and    (  ) the spectral coefficients.Thus, Equation 47 is written as a system of linear equations where the vector − →  has as entries the coefficients    (  ), L is a matrix dependent on the basic state and − →   describes the inhomogeneity of the equation.We refer to Equation 48 as the linear model associated with the PGE model from section 3.3.

Standard set-up
In the following, we refer to the set-up of the source terms described in the previous section as the standard set-up.The resulting stationary zonally averaged circulation in the model shows realistic key features of the temperature and wind fields as displayed in Figure 2a.The westerly jets are in thermal wind balance and have a jet maximum of 40 m⋅s −1 at around 30 • at the model top.At the surface, weak easterlies in the tropics and westerlies in the midlatitudes are visible.
The baroclinic wind u ′ (not shown) is easterly in the lower atmosphere and westerly in the upper atmosphere and has a linear vertical shear.In order to study the meridional overturning circulation, we introduce a stream function  satisfying The Hadley cell depicted in Figure 2b is too broad and extends into the high latitudes.Note that due to Equation 41 the meridional circulation is explained entirely in terms of viscous axisymmetric models (e.g., Schneider and Lindzen, 1977) and is missing the important advection of relative angular momentum (see Held and Hou, 1980 for a discussion).
The stationary zonal anomaly of the temperature field in the lower and upper atmosphere is displayed in Figure 3.Each individual extremum can be associated with an extremum in the forcing and there is no generation of wave trains as in the quasi-geostrophic dynamics.There is almost no vertical tilt of the disturbances but there is a phase shift of about 30 • to the east with respect to the forcing.The stationary wave amplitude is small and decays rapidly with height.The wave amplitude is considerably underestimated when compared with quasi-stationary waves in the atmosphere.This discrepancy might result from the Boussinesq approximation and from the vertical profile of the forcing considered here.
The time-averaged meridional momentum transport is analysed in Figure 4a.For that purpose the zonally and vertically averaged momentum flux ⟨ ′  ′ ⟩ , , entering the barotropic Equation 39, is separated into contributions from the zonally symmetric circulation (mean meridional overturning circulation) and from the zonally asymmetric part, defined as eddies.Figure 4a shows that the meridional momentum transport by the eddies is not significant.Whereas the magnitude of the momentum flux by the zonally symmetric circulation is realistic, the one by the eddies is considerably underestimated when compared with observations (Peixoto and Oort, 1992).
Next, we consider the budget in the barotropic momentum equation, which determines the zonally averaged surface zonal wind.By dividing Equation 28 by acos and substituting the source term S m corresponding to Equation 37, one obtains the following barotropic momentum equation for stationary motion: The term on the left-hand side describes the contribution from the momentum flux divergence and from the metric term (the fourth term on the left-hand side of Equation 20), while the two terms on the right-hand side account for Ekman friction by the barotropic and baroclinic flow.The contributions of the different terms are displayed in Figure 4b.The barotropic wind and the baroclinic surface wind have opposite signs everywhere outside the tropics and there is considerable cancellation between the two components.The Ekman friction by the full surface zonal wind ⟨u(0)⟩  = ⟨u⟩ z + ⟨u ′ (0)⟩  balances the momentum flux term on the left-hand side of Equation 49 (except at the Poles due to interpolation errors).
The PGEs generate transient disturbances by baroclinic instability (see section 5.5), which propagate with the mean flow and are concentrated in the subtropical and tropical regions.Those transients are damped in the standard model configuration by choosing sufficiently high diffusivity  b 2 .We also performed simulations allowing for baroclinic eddies with zonal wavenumber 10 and obtained qualitatively similar results (not shown) as described in this section.

4.2
The effect of the barotropic closure In order to assess the effect of the evolution equation for the barotropic flow on the circulation, we perform simulations where the closure (Equation 39) is omitted.Without closure, ⟨u⟩ z does not change from its initial value and the total zonal mean zonal wind in the model is too strong and shows super-rotation at the Equator; see Figure 5. Whereas the zonally symmetric baroclinic flow remains unchanged, the amplitude of the stationary zonal perturbations is reduced by a factor of 50% if the temperature fields are considered (not shown).The latter are affected by the barotropic wind through buoyancy advection.

Sensitivity with respect to the diffusion coefficient K
Since in the free atmosphere the effect due to eddy dissipation should be negligible, we perform simulations with non-uniform diffusion coefficient K in Equation 41.In particular, we choose for the meridional dependence of K a Gaussian profile centred at the Equator: 2 The value of  b considered here (see Table 1) corresponds to a damping time-scale of three quarters of a day of the spherical harmonics with total wavenumber n = 21.F I G U R E 5 Time mean zonal mean zonal wind (contours) and potential temperature (shading) in a simulation without the closure Equation 39with K 0 = 5m 2 s −1 , corresponding to the uniform K value in the standard set-up.By setting  K = 4 • , K decreases rapidly away from the Equator.Note that the equatorial region should have non-vanishing K to prevent the singularity of the PGEs discussed earlier.Almost no difference is visible in the zonally averaged temperature and wind fields from Figure 6a compared to the standard set-up (Figure 2a).There is some weakening of the Hadley cell observed; see Figures 6b and 2b.
We have also performed simulations with uniform K, but taking one tenth of the reference magnitude 5 m 2 ⋅s −1 .The resulting circulation (not shown) was almost unchanged compared to the standard set-up, and we conclude that the results are not sensitive with respect to the diffusion coefficient K. Simulations were also carried out by replacing the diffusion in Equation 41by Rayleigh friction and qualitatively similar results (not shown) were obtained.

Sensitivity with respect to resolution
The effect of the upper boundary condition (Equation 36) on the dynamics is studied by performing a model simulation with a doubled number of vertical levels (10 levels with Δz = 1 km), where the stratosphere is resolved by two layers instead of a single layer in the standard set-up.In the stratosphere the vertical structure of the baroclinic winds − →  ′ is set to a linear profile.The latter is determined by imposing ⟨ − →  ′ ⟩  = 0 and requiring continuity of − →  ′ at the tropopause.
In the special case of a single stratospheric layer this approach reduces to the one described by Equation 36.The convergence of the results with respect to the horizontal resolution

F I G U R E 6
The same as Figure 2, but for a simulation with the diffusion coefficient K confined in the tropics; see Equation 50was verified by considering a T42 spectral resolution.The zonally averaged circulation in the model is summarized in Figure 7.All the main features in the circulation remain the same as in the standard set-up.The maximum of the stream function  max from the meridional overturning circulation (Figure 7b) is slightly reduced from 2,172.2 m 2 ⋅s −1 (standard set-up) to 1,999.8 m 2 ⋅s −1 .No changes in the zonally asymmetric disturbances are observed (not shown).
In addition, a model run with the height of the troposphere extended to z t = 10 km (z a = 12 km) was performed where the total number of vertical levels was increased to 12 (Δz = 1 km).In the simulation the westerly jets increase further above 10 km and reach a maximum of about 52 m⋅s −1 (not shown).There is an intensification of the Hadley cell ( max = 2,306.8m 2 ⋅s −1 ), where the upper branch shifts to higher altitudes (not shown).

LINEAR ANALYSIS
Many aspects of the forced stationary waves in the atmosphere can be explained using linear theory within the quasi-geostrophic framework (e.g., Held, 1983;Pedlosky, 1987).For the PGEs, topographically and thermally forced stationary wave solutions were recently presented by Egger and Hoinka (2017).Here wave solutions of the linear PGE Boussinesq model are used to interpret the results from the numerical simulations in section 4.
In the following we consider the linearized equations with  b = K = C D = 0 and a basic state from Equations 44-46 with  = (),   = const and   = 0. Differentiating Equation 47 with respect to z, using the thermal wind relation and hydrostatic balance, yields where Φ* denotes the perturbation of Φ ′ from Φ ′ and  *  = ( *  −  * )∕.

Free waves
Looking for solutions of the form Φ * = Φ() exp{( +  − )} (with k the zonal wavenumber, m the meridional wavenumber and  the frequency) and setting  *  = 0, one obtains the dispersion relation This corresponds to the long wavelength limit of Rossby waves from quasi-geostrophic theory.The waves become stationary if  =   2  2   .Note that the left-hand side of Equation 51 does not involve any meridional derivatives of Φ* and the equation decouples in the meridional direction.As a consequence, the waves can have arbitrary meridional structure.

Forced stationary waves: General case
We consider forcing of the form with  *  = () cos( 0 ) − , which has the form of the zonally asymmetric forcing (Equation 33).The stationary form of Equation 51 divided by ∕ cos  is where the following definitions were introduced:

F I G U R E 7
The same as Figure 2, but for a simulation with T42 spectral resolution and 10 vertical levels The rigid lid boundary conditions for Φ* are given by the linearized version of Equation 47 under the assumption of   = 0 and  b = 0: A particular solution of the inhomogeneous Equation ( 54) reads where The homogeneous solution to Equation 54 takes the form where the real numbers  r ,  i satisfy  =  r + i i with  2 = −  1+ 2 (1 + ).The particular solution (Equation 59) alone does not fulfil the boundary conditions (Equation 58), but together with the homogeneous solution they can be satisfied by setting k = k 0 in Equation 62 and choosing appropriate constants a j .However, the explicit form of the constants soon becomes tedious and in the following section we introduce an approximation in order to simplify the analytical expressions.

Forced stationary waves: No-relaxation case
Neglecting the damping term in the buoyancy forcing, we consider here In this case we can set γ in Equation 54to zero and the particular solution has the form Again we have to add the corresponding homogeneous solution in order to satisfy the boundary conditions.The full solution takes the form

Forced stationary waves: Comparison with the nonlinear model
We compare the solutions of the linearized equations described in sections 5.2 and 5.3 with the full nonlinear stationary model response from section 4.1.We consider altogether three models describing linear dynamics.The first model is the analytical solution (Equation 65) evaluated by setting the basic state zonal wind () to the time-averaged zonal mean zonal wind at 3 km height from the nonlinear simulation and setting the buoyancy vertical gradient  48, with  = 0, where the basic state entering L is the same as for the analytical solution from Equation 65and  b = K = C D = 0. We refer to the resulting model as the linear inviscid model.Note that there is no vertical shear in the basic state zonal wind in this model.The third model is the stationary solution of the linear model from Equation 48but with a vertically varying basic state, where (, ⟨⟩  ) are set to the time-averaged profiles (⟨b⟩  , ⟨u⟩ z ) from the nonlinear simulation.In addition, effects due to friction and diffusion are taken into account in the linear model by using the values of  b , K and C D from the standard set-up.Due to numerical reasons the solution of the linear inviscid model is computed with T42 spectral resolution, whereas for the full linear model T21 resolution is sufficient.
The results for the different models are summarized in Figures 8 and 9.The linear model (Figures 8d and 9d) reproduces the meridional and vertical structures of the stationary waves in the nonlinear simulation (Figures 8a and 9a).Small deviations in the meridional structure are visible only equatorward of 30 • .The analytical solution (Equation 65) captures the large-scale structure of the stationary waves poleward of 30 • but produces spurious oscillations in the tropics (Figure 8b).Similar oscillations are observed in the inviscid linear model due to the neglected frictional effects (Figure 8c).
Figure 10 shows the stationary linear solutions when the relaxation time-scale  is reduced from the standard value of 15 days to 5 days.For  = 5 days the full linear model again accurately reproduces the nonlinear response (not shown).The magnitude of the analytical solution increases for smaller  in accordance with the nonlinear solution.Further inspection shows that there is a slight increase by about 5 • in the phase lag of the analytical solution with respect to the nonlinear solution.The inviscid linear model captures the nonlinear response to diabatic forcing for  = 5 days.The exact match of the magnitudes at 50 • is to some extent by chance (but the phase match is not).The reason for this is that the basic state wind in the inviscid linear model was set to the stationary wind field at 3 km height from the nonlinear simulation, which is an arbitrary level choice.From Figure 10 we conclude that for smaller  the difference between the linear inviscid model and the analytical solution increases.At the same time, the effect of the vertical shear (present in the linear model and not included in the linear inviscid model) becomes less important over the relaxation effect (included in both linear models).

5.5
Baroclinic instability within the PGEs on the sphere Until now we have considered disturbances generated by diabatic heating, but even in the absence of forcing the PGEs can produce exponentially growing disturbances by the mechanism of baroclinic instability.To study the latter process we utilize the linear model from Equation 48 without forcing and dissipation, that is,  b , K, C D and 1/ are set to zero.As a basic state we set (, ⟨⟩  ) to the time-averaged profiles (⟨b⟩  , ⟨u⟩ z ) from the nonlinear simulation.The results reported here are computed using T85 spectral resolution for convergence reasons.
Growth rates of the most unstable modes are shown in Figure 11a.The growth rates increase linearly with the zonal wavenumber without any bound.This is consistent with previous -plane analyses of the PGEs (Wiin-Nielsen, 1961;Colin de Verdiere, 1986).Since the PGEs are valid on very large spatial scales, only the results for the lowest wavenumbers should be relevant for the atmosphere.For example, for wavenumber 3 the growth rate corresponds to an e-folding time-scale of about 1 week, which is comparable with the time-scale of radiative processes.As shown in Figure 11b, the phase speed of the unstable modes nearly does not depend on the wavenumber and takes values between 1 and 4 m⋅s −1 .Due to the meridional decoupling of the inviscid eigenvalue problem, as discussed below in Equation 52, the horizontal structure of the unstable modes cannot be determined.The disturbances show a westward vertical tilt of about quarter of a wavelength in the lowest atmosphere; see Figure 12.

CONCLUSIONS
We present numerical simulations of the PGEs for Boussinesq fluid on the sphere supplemented by a novel evolution equation for the barotropic flow.The latter is affected by meridional momentum flux due to baroclinic flow and drag by the surface wind; see Equation 39.On the other hand, the barotropic wind affects the baroclinic flow through buoyancy advection.In order to remove the singularity of the PGEs at the Equator, the geostrophic balance is modified by including turbulent eddy diffusion.This is a different approach compared to other PGE-type models, where f is fixed at a constant value f (±15 • ) in the tropical region of each hemisphere (Petoukhov et al., 2000).The model is forced by relaxation towards a prescribed buoyancy profile.The model climatology shows westerly jets and surface tropical easterlies consistent with other Boussinesq simulations (e.g., Held and Hou, 1980).Due to the inclusion of turbulent eddy momentum diffusion, the model produces a viscous Hadley cell.This overturning circulation is responsible for the meridional momentum transport, whereas the flux due to eddies is negligible.The stationary zonally averaged surface zonal wind is determined entirely by the baroclinic meridional momentum flux ⟨u ′ v ′ ⟩ z,  (see Equation 49).There is considerable cancellation between the barotropic wind and the baroclinic surface zonal wind when time and zonal averages are considered (see Figure 4b).It is observed We study the response of the model to an idealized land-sea thermal forcing with exponential vertical decay.The stationary waves observed in the simulation are confined to the lower atmosphere and have no vertical tilt.It is shown that the response can be understood entirely in terms of linear dynamics.Forced stationary wave solutions within the PGEs were derived by Egger and Hoinka (2017).Here we consider analytical solutions but for different forcing profiles and under the Boussinesq assumption.It is shown that those solutions reproduce key features of the vertical and horizontal structures of the model response in midlatitudes.
The analysis of Wiin-Nielsen (1961) on baroclinic instability within the PGEs on a -plane is extended to the sphere by considering the growth rate, phase speed and vertical structure of the most unstable modes.The growth rates increase linearly with the wavenumber.This unbounded increase is due to the neglected relative vorticity advection in the PGEs (Wiin-Nielsen, 1961;Colin de Verdiere, 1986) and makes the numerical treatment of the equations challenging, since the highest resolved scales are the most unstable.In our model the inclusion of buoyancy diffusion introduces a cut-off in the growth rates.In the standard model configuration, the baroclinic eddies are suppressed using sufficiently high diffusion.Simulations with baroclinic eddies, however, indicate qualitatively similar results for the zonally averaged circulation.
Due to the Boussinesq assumption the wave disturbances (forced and baroclinic) in our model do not show an increase of amplitude with height as typically observed in the atmosphere.Consequently, momentum and temperature transport by the waves is underestimated.In the future we plan to relax the Boussinesq approximation to account for the missing effect.This requires further analysis to pose an appropriate upper boundary condition for the model.
Another important ingredient absent in the present model is the midlatitude synoptic-scale dynamics.In the case of small-amplitude eddies, Boljka and Shepherd (2018) and Boljka et al. (2018) provide a framework for studying interactions of the mean flow with planetary and synoptic scales.In the case of large-amplitude eddies, the two-scale model of Dolaptchiev and Klein (2013) would be the asymptotic consistent extension of the present planetary scale model to the synoptic scale.Interestingly, the planetary barotropic flow equation in that case provides the only feedback mechanism from the synoptic scale to the planetary scale.This highlights the importance of the barotropic closure equation considered here for the dynamics.
a) Zonally symmetric potential temperature (shading) and zonal wind (contours) corresponding to the relaxation profile from Equation 31.(b) Zonally asymmetric potential temperature distribution at a height of 1 km corresponding to the relaxation profile from Equation imposed in the model by setting the baroclinic horizontal velocity at the upper model level z s to Time mean zonal mean circulation in the PGE model: (a) zonal wind (contours) and potential temperature (shading); (b) stream function of the meridional overturning circulation Time mean zonally asymmetric potential temperature at heights of (a) 1 km and (b) 9 km a) Zonally and vertically averaged time mean meridional momentum flux ⟨u ′ v ′ ⟩ z,  by the mean meridional overturning circulation (MMC) and by the eddies (10 4 EDD), where the magnitude of the eddy flux has been multiplied by a factor of 10 4 to make it visible on the scale.(b) Contributions in the stationary barotropic momentum equation (Equation49): Ekman friction by the baroclinic (EKZ) and barotropic (EKB) surface winds, and contributions from the momentum flux divergence and metric term (MFD) and the residuum (RES).See the text for details

)
Time mean zonally asymmetric potential temperature at a height of 1 km: (a) nonlinear simulation; (b) analytical solution (Equation 65); (c) inviscid linear model and (d) linear model.See the text for explanations of the different models.Note also that in Figure 8b the ±2 contour line is drawn to indicate the large amplitudes around 15 • Longitude-height cross-section at 50 • N of time mean zonally asymmetric potential temperature: (a) nonlinear simulation; (b) analytical solution (Equation 65); (c) inviscid linear model and (d) linear model to ⟨  ⟩   =   ∕  .The second model is described by Equation 54, but instead of solving Equation 54 for Φ* we solve the equivalent equation for b*.This has the form of Equation Time mean zonally asymmetric potential temperature as a function of longitude at 50 • N and at a height of 1 km from the nonlinear model (solid line), the linear inviscid model (dotted line) and the analytical solution (dashed line), for a relaxation time-scale of (a11 (a) Growth rate (day −1 ) and (b) phase speed (m⋅s −1 ) corresponding to the most unstable mode as a function of the Longitude-height cross-section at 50 • N of the real part of the most unstable mode with the zonal wavenumber 3, in units of  0 that the barotropic wind affects only the zonally asymmetric part of the baroclinic flow (see section 4.2).