Bedymo: a combined quasi-geostrophic and primitive equation model in sigma coordinates

. This paper introduces the idealised atmospheric circulation model Bedymo, which combines the quasi-geostrophic approximation and the hydrostatic primitive equations in one modelling framework. The model is designed such that the two systems of equations are solved as similarly as possible, such that differences can be unambiguously attributed to the different approximations, rather than the model formulation or the numerics. As a consequence, but in contrast to most other quasi-geostrophic models, Bedymo is using sigma-coordinates in the vertical. In addition to the atmospheric core, Bedymo 5 also includes a slab ocean model and passive tracer module that will provide the basis for future idealised parametrisation of moisture and latent heat release (cid:58)(cid:58)(cid:58) with (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) options (cid:58)(cid:58) for (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) prescribed (cid:58)(cid:58)(cid:58)(cid:58) and (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) wind-induced (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) currents. Further, Bedymo has a graphical user interface, making it particularly useful for teaching.


Introduction
Since the 1950ies :::: 1950s, several studies introduced different quasi-geostrophic (QG; e.g.Charney and Phillips, 1953) and hydrostatic primitive-equation models (PE; e.g.Smagorinsky, 1958).Table 1 of Schär and Wernli (1993) lists a more recent selection of QG and PE models, and the development of idealised models is still ongoing (e.g.Hogg et al., 2003;Maze et al., 2006).Given this wealth of existing models, why should one develop another model?The main reason for developing the BErgen DYnamical MOdel (Bedymo) is to combine two approximations -quasi-geostrophy and the dry hydrostatic primitive equations -in one modelling framework.To our knowledge, there is no other model that combines several approximations in one numerical framework.Hence, Bedymo provides the only dynamical core that incorporates two levels in the Held (2005) hierarchy of models.
Typically, the approach to solve the underlying set of equations is different for each approximation.For example, QG models usually forecast the QG potential vorticity (QGPV) and diagnose the geostrophic wind and temperature field by inverting the temperature distribution fully determines the atmospheric state in the QG system via hydrostatic and geostrophic balance.In PE, the horizontal wind velocity components evolve independently and hence must be forecasted as well :::::::: separately.
In both systems, pressure and σ vertical velocity is required to integrate the thermodynamic equation forward in time.In QG, the pressure vertical velocity follows from an inversion of the omega-equation, which implicitly establishes the threedimensional QG-balance, and σ-vertical velocity is derived therefrom.In PE, both vertical velocities are derived from the continuity equation, using the divergence of the forecasted horizontal flow and the local surface pressure tendency.In PE, the local surface pressure is given by the column-integrated mass flux divergence and thus also follows from continuity.In contrast, in QG the local surface pressure tendency is derived from the pressure vertical velocity at the lower surface, and thus by the lower boundary condition used for the ω-inversion.This boundary condition is given by the vorticity equation evaluated at the lower surface such that in QG surface pressure evolves following QG vorticity dynamics.
In the following, we summarise the model equations for QG and PE.We use a β-plane approximation for the Coriolis parameter f = f 0 + βy.In order to derive a self-consistent QG system, we need to approximate the specific volumes by timeinvariant and horizontally homogeneous background values α = α(σ).This approximation is optional for PE, and we will introduce and evaluate both variants.All symbols used in the following equations are summarised in Tab. 2.

Full primitive equations
The PE equations, as used in Bedymo, are In order to provide an energy sink for long and short waves, respectively, the momentum equations include scale-inselective linear Ekman friction and a scale-selective damping term.
In this system, we infer geopotential φ by integrating the hydrostatic equation ( 4) upwards, starting from the time-invariant surface geopotential φ s which represents the model orography.The local surface pressure tendency results from continuity, And, finally, which closes the system of equations ( 1)-( 6).All other parameters in the equations are either time-invariant fields, constants, or represent an external forcing.

Homogeneous density approximation
As a first step towards quasi-geostrophy, we first assume horizontally homogeneous density within the PE system.In mathematical form, the approximation α ≈ α(σ) resembles the anelastic approximation in Cartesian z-coordinates.However, while the continuity equation for the anelastic approximation reduces to that for incompressible flow, the continuity equation in σcoordinates is unchanged by this approximation.As a result of this difference, (barotropic) acoustic waves are still supported, while they are not with the anelastic approximation.
As a result of the homogeneous density approximation, the pressure gradient terms become linear, because ασ only varies with height.In fact, using the ideal gas law ασ = RT ps , the product ασ becomes constant with an isothermal background state.Further, vertical temperature advection is approximated by advection of the homogeneous background state only, σ ∂T ∂σ ≈ σΓ with Γ = ∂T ∂σ .
In summary, the resulting momentum and thermodynamic equations are Continuity (eq.6), hydrostasy (eq.4) and eqs.( 7)-( 9) that determine the surface pressure tendency and the vertical velocities all remain unchanged.

Quasi-geostrophy
In QG, the geostrophic wind follows from the geostrophic streamfunction ψ g , which is in turn defined by The prognostic equations are in which the surface pressure tendency equation ( 16) is evaluated at the lower surface (σ = 1), indicated by the subscripts s.
The remaining unknown variables at this point are the vertical wind σ and the pressure vertical velocity ω.Starting from the following form of continuity the pressure vertical velocity can be derived from the ω-equation in which Q 1 and Q 2 represent the two components of the Hoskins et al. (1978) Q-vector Further, s = s(σ) is a horizontally homogeneous stability parameter In order to solve the elliptic ω-equation (18), we require a condition for ω at the lower surface.We obtain ω s from the vorticity tendency equation evaluated at the lower surface.The resulting condition is in which the subscript s again indicates variables that are evaluated at the surface.Finally, after using ω s to obtain ω for the entire model domain, we derive σ from which closes the system of equations.
QG potential vorticity q is not used anywhere in the model.Nevertheless, for reference, q takes the form in this QG system.

Ocean dynamics
In addition to the atmospheric component, Bedymo also includes a slab ocean.The slab ocean is intended to represent an oceanic mixed layer that interacts with the atmosphere on time scales on which the internal ocean dynamics can be neglected.
Nevertheless, as the oceanic heat transport might play a role even at these time scales, we provide several options to provide oceanic flow and heat transport.
The only prognostic variable in the slab ocean model is the mixed layer temperature T o , It can change due to sensible heat exchange F sh between the ocean and the atmosphere and oceanic heat transport OHT .
In addition, the model includes a relaxation term towards a prescribed climatological temperature T o e , which may be used to crudely represent the neglected oceanic circulation (see Tab. 2 for other symbols).
The heat exchange is parameterised with a bulk flux formulation, Here and in the following, atmospheric variables are denoted by a superscript a, and the index s denotes values at the interface between the atmosphere and ocean.Details on how these surface variables are defined are given in section 2.4.
The options for parametrising the oceanic heat transports are for a 1-layer model, and for a 1.25-layer model. ( The first option represents a motionless ocean, while the other two options differ in the way divergence in the flow field is treated, with u o denoting the 3D-flow field and v o the horizontal part of the flow.Whereas divergence does not influence the mixed layer temperature in the 1-layer model, the 1.25-layer model assumes a compensating vertical flow (positive upwards) that transports water from a lower layer of temperature T o 2 into the mixed layer, The temperature T o 2 might vary spatially, but is kept constant over time.By not letting T o 2 vary with time, we implicitly assume the deeper layer to be motionless and to have infinite heat capacity.
Formally, both the 1-layer model and the 1.25-layer model for the OHT do not conserve energy, because the average mixed layer temperature can change without temperature relaxation or heat exchange with the atmosphere.Upwelling in the 1.25layer model lets the average mixed layer temperature approach the deep layer temperature, without the latter changing due to the downwelling required by continuity.The 1.25-layer model reduces to the 1.0-layer model if T o 2 = T o , such that the average mixed layer temperature will over time approach the local mixed layer temperature in locations of divergent flow.
The total transporting flow v o = v e + v p is a combination of a user-prescribed flow pattern v p and the wind-driven Ekmanflow v e .Our formulation of v e also includes a small parameter in addition to the standard Ekman flow, which controls the magnitude of the non-rotating flow excited by wind stresses.This addition results in finite Ekman velocities at the equator.Codron (2012) refers to as the inverse damping timescale of oceanic currents.

Coordinate system and discretisation
Bedymo uses Cartesian coordinates in the horizontal and a terrain-following pressure coordinate σ in the vertical, which is defined by Variations in surface pressure thus affect the coordinate system throughout the atmospheric column up to the model top at In both the horizontal and the vertical directions, the velocity components are staggered with respect to the main grid points, following the C-grid setup of Arakawa and Lamb (1977).Temperature, geopotential, specific volume, and surface pressure are all defined on the unstaggered grid.

Elliptical solver
In QG-mode, an elliptic equation has to be solved to determine the vertical velocity.Through extensive testing, we found the Full Multigrid Method (e.g.Saad, 2003) in conjunction with the Bi-Conjugate gradient Stabilised (BiCGstab) Method of van der Vorst (1992) to be the most effective and stable configuration for Bedymo.We use BiCGstab both to solve the elliptic equation on the coarsest grid and to iteratively refine the solution on finer grids.

Boundary conditions
The boundaries of the model domain are located at staggered grid points.At the model top, p = 0 and also both vertical velocities vanish, σ = 0 and ω = 0.At the surface, σ = 0, but ω s = 0 as surface pressure can change both locally and in a Lagrangian reference frame.In QG, ω s is determined through a condition derived from the surface vorticity tendency (eq.22); in PE, it is derived from continuity (eq.9).
There are two available options for boundary conditions along the lateral boundaries.The first option represents periodic boundaries in the respective direction, and the second option features an impermeable free-slip "wall" with zero boundarynormal fluxes.

Evaluation
We subject Bedymo to several test cases to evaluate the performance of the model.We chose different test cases to isolate pertinent aspects of the atmospheric and coupled dynamics.Furthermore, we chose test cases that have been studied comprehensively, such that the expected results are well established.The first three ::: four : test cases focus on important aspects of the mid-latitude atmospheric dynamics in isolation.In a fourth :::: fifth test case, we evaluate the PE model in conjunction with the slab ocean by simulating the coupled response to a temperature anomaly in the ocean mixed layer located at the equator.:: A ::::::: summary ::: of ::::::: pertinent :::::: model ::::::::: parameters ::: for :: all :::: test :::: cases :: is ::::: given :: in :::: Tab.:: 3.

Atmospheric test cases
The first atmosphere-only test case evaluates the representation of baroclinic cyclogenesis in a mid-latitude channel on the f and β-plane, as well as the sensitivity of the baroclinic development to the magnitude of the initial baroclinicity.Second, we extend the baroclinic channel setup by including temperature relaxation and friction and evaluate the long-term zonalmean statistics of the storm track simulated by Bedymo.Finally, test case three evaluates :::: cases ::::: three ::: and :::: four :::::::: evaluate the development towards the stationary barotropic Rossby wave response :::: wave :::::::: responses : to isolated orography : , :::: with : a ::::: focus ::: on
This corresponds approximately to an overall Brunt Vaisala frequency N 2 = 1.2•10 −4 s −2 .The baroclinic zone is meridionally centred in the channel.There is no initial surface pressure gradient and thus no initial surface winds.With Coriolis parameters corresponding to a latitude of 50 • N, the baroclinic zone is initially balanced by a jet of maximum intensity of about 50 m s −1 in the uppermost layer (σ = 1/6).
To localise the baroclinic development, we perturb this initial state by a warm anomaly of 2 K located at the center of the baroclinic zone.The temperature perturbation is invariant with height and has a zonal and meridional extent corresponding to the width of the baroclinic zone (approx.3000 km).The perturbation is initially balanced by a perturbation thermal wind.
In our control setup, we use the β-plane approximation, as it is this setup we will later extend to multi-decadal storm track simulations.We will compare this control setup to simulations on an f -plane and to simulations in which we vary the initial baroclinicity to yield a balanced initial jet of either about 30 m s −1 or 70 m s −1 , respectively.The baroclinic development for the control setup and the three model variants is summarised in Figure 2.For the first three days, the evolution in the three model variants is quite similar (Fig. 2a-f), and only at day 5 structural differences become apparent between QG and the two PE variants (Fig. 2g-i).In QG, cyclones are considerably larger in scale than in the PE variants and largely symmetric in size and structure compared to the anticyclones.In contrast, PE cyclones are smaller, more circular, and surrounded by steeper pressure gradients than PE anticyclones.Further, only PE produces qualitatively realistic fronts with near-discontinuities in the temperature field.This is expected from theory, as Hoskins (1975) showed that advection by ageostrophic winds must be taken into account in order to realistically capture the evolution of fronts.
Comparing the two PE variants, only minor differences occur up to day 5 (Fig. 2h-i).The two most apparent differences are a slightly slower cyclone intensification with homogeneous density PE, and, potentially related, changes in the temperature structure.In the homogeneous density PE, cyclone cores are comparatively warmer, and anticyclones comparatively cooler than their counterparts in full PE.Given that the only differences between the two PE variants is linearised vertical advection and a linearised pressure gradient force, it seems implausible that these differences in temperature originate from the linearised vertical advection.Warmer cyclones and colder anticyclones yield an overall somewhat reduced baroclinicity in homogeneous density PE compared to full PE, which is consistent with the slightly slower intensification in homogeneous density PE.
Considering the growth rate of the eddy kinetic energy, both variants of PE develop somewhat slower and later than QG (Fig. 3), with maximum growth rates of 17.7 hours (QG), 19.0 hours (full PE), and 21.1 hours (homo.PE).For our channel setup, the Rossby deformation radius is approximately 930 km (e.g.Vallis, 2006), which implies a most unstable wave length of about 3600 km, and a maximum Eady growth rate of for the control jet intensity of 50 m s −1 .This corresponds to an e-folding time scale of τ E = 16.7 hours, which fits quite well, in particular with the QG variant.
Further, the simulated evolution generally follows the expectations from theory for varying magnitudes of the baroclinicity.
For a stronger jet of 70 m s −1 , the Eady growth time scale reduces to 11.9 hours.This fits almost perfectly to the simulated maximum intensification with a time scale of 11.7 hours in the QG variant (darkest blue line in Fig. 3).The two PE variants are again somewhat slower and later in their development (lighter blue lines in Fig. 3).For the weaker jet of 30 m s −1 , the correspondence between theory and simulation is less tight, with an Eady growth time scale of 27.8 hours and a QG maximum intensification with a time scale of 35.6 hours.This might, at least partly, be due to the beta-effect and downstream development   playing a comparatively larger role with weaker baroclinicity.But this cannot be the full explanation, because for full PE the 265 reduction in the peak growth rate (about 44%) is only slightly larger than the reduction in baroclinicity (40%).
Despite these differences in growth rate, structurally the evolution remains very consistent across the three tested magnitudes of the baroclinicity (Fig. 4).Considering lead times normalised by the different jet intensities, the wave structures are nearly identical within each model variant (columns in Fig. 4), showing that downstream propagation of wave energy scales largely linearly with the jet intensity despite the considerable amplitude of the waves.
The evolution of the cyclone is qualitatively similar on the f and β-plane for all three model variants (compare Figs. 2 & A1).
In particular the differences in the size and shape of cyclones and anticyclones translate from the β to the f-plane.However, at day 5 it becomes apparent that meridional movement is much less constrained on the f-plane compared to the β-plane.The meridional scale of the synoptic systems is markedly larger, and soon after day 5 the large cold sectors in the zonal center of the domain start to interact with the southern boundary of the domain.
Overall, the structure and sensitivities of the simulated cyclone development is very much in line :::: with previous simulations of idealised cyclones with more complex models (e.g., Schemm et al., 2013, Terpstra andSpengler, 2015, and the cyclones in the colder environments of Tierney et al., 2018).This consistency is not surprising given the much earlier results of for example Simmons and Hoskins (1978), :::::::::::::::::::::::: Simmons and Hoskins (1978), ::: for :::::::: instance, which used a model similar in complexity to Bedymo to study cyclogenesis.But the consistency remains worth noting, because these more recent studies typically employ a factor of 5-10 higher resolution in both the horizontal and vertical and use a full suite of physics parametrisations.

Mid-latitude storm track
The control setup for the baroclinic instability test case also serves as the basis for long-term simulations to evaluate the representation of a mid-latitude storm track in Bedymo.In order to achieve a statistically stationary storm track, we follow Held and Suarez (1994) and add simple parametrisations for three physical processes to the model setup.First, we enable a temperature relaxation towards the initial state with a time scale of 10 6 s ≈ 11.6 days throughout the model domain to represent all process that replenish baroclinicity in real storm tracks.Second, we enable linear Ekman friction in the lowest model layer with r = 10 6 s.Finally, we enable bi-harmonic diffusion that damps 2∆x-waves with a time scale of 10 5 s ≈ 28 hours.This time scale increases quadratically with increasing wave length and thus predominantly affects the shortest waves.The storm track simulations cover a time period of 32 years of 360 days each, of which the first two years are discarded as spin-up.After 1-1.5 years there is no discernible trend anymore in the distribution of sea-level pressure in any of the simulations, which we take as an indication that a statistical equilibrium has been reached.
Because of the absence of any zonal asymmetries, the resulting storm track is ::::::: statistics ::: are zonally symmetric.Figure 5 thus presents the zonal and time average state for the last 30 years of the simulations.In addition to the average temperature and zonal winds, we also show the eddy momentum and heat fluxes, in which the temperature and wind perturbations (indicated by primes) are taken to be the deviations from the zonal-and-time average.These fluxes are interesting, because, for example ∂ ∂y u v < 0 indicates a convergence of eddy momentum fluxes in the time mean, and analogously for ∂ ∂y v T < 0 and ∂ ∂σ σ T < 0 that indicate a convergence of heat in the meridional and vertical, respectively.To keep the line plots legible, we restrict our presentation and discussion to QG and full PE.
In QG, the resulting storm track is symmetric and centred on the imposed baroclinic zone.In the angular momentum budget, weak surface easterlies on either side of the baroclinic zone balance the near-surface westerlies under the jet (Fig. 5a).
In contrast to QG, the simulated PE storm track is not symmetric around the imposed baroclinic zone.The core of the jet 315 is shifted poleward at all levels, and near-surface easterlies mainly occur on the equatorward side of the imposed baroclinic zone.This asymmetry is likely related to the asymmetry between cyclones and anticyclones observed in the cyclogenesis test case (Fig. 2).As discussed there, the strongest gradients in both sea-level pressure and temperature occur around cyclones, and thus on average somewhat poleward of the imposed baroclinic zone (Fig. 2).It thus seems plausible that the time mean reflects this asymmetry.If this interpretation is correct, the asymmetric storm track would thus be another consequence of taking into account advection by the ageostrophic winds (cf., Hoskins, 1975;Wolf and Wirth, 2015).
The storm tracks observed on Earth display a similar meridional asymmetry with regards to the surface winds.Near-surface easterlies are much more pronounced in the subtropics, on the equatorward side of the baroclinic zone, compared to the polar regions.Despite the similarity, the mechanisms leading to these easterlies might nevertheless be different.For real storm tracks, both , ::::: with ::::::: spherical ::::::::: geometry ::: and : the Hadley circulation and the spherical geometry of the Earth certainly introduce meridional asymmetries that cannot be captured in a Cartesian mid-latitude channel ::::: again ::::: likely ::::::: affecting ::: the :::::::: dynamics.In fact, some setups of the spherical 1-layer QG model of Vallis et al. (2004) yield meridional asymmetries similar to the one seen here for PE.
In addition to the asymmetries, the PE storm track features larger heat and momentum fluxes at all levels compared to QG.
Nevertheless, the meridional profiles of these fluxes are generally consistent across QG and PE at all levels.The only exception is near-surface eddy momentum transport, but here amplitudes are quite small both in QG and PE.Consistent with the more vigorous heat transport in PE, the average temperature deviations from the relaxation state are larger in PE than in QG.On the equatorward side of the mean jet, weak but noticeable temperature deviations extend all the way to the southern boundary.
Interactions with the boundary dohowever : , ::::::: however, : not affect the PE storm track, as the meridional profiles of all parameters remain nearly unchanged in a simulation with a meridionally wider channel.
Finally, the simulated storm track remains largely consistent across a range of temperature relaxation time scales (Fig. 6).In addition to the default relaxation with a time scale of 11.6 days, we conducted simulations with full PE using time scales of 3.47 and 34.7 days, respectively.These simulations are called "weak ::::: strong" and "strong ::::: weak" in Fig. 6, respectively.Within this range of parameters, stronger relaxation yields a more vigorous storm track with more intense heat and momentum transports as well as stronger jets (Fig. 6).

Rossby waves excited by orography
The third and final mid-latitude test case evaluates the developing stationary wave response to isolated orography.Orographically forced stationary waves are one of the main ingredients determining the zonal asymmetries in the Northern Hemisphere storm tracks (Held et al., 2002).The test case considered here uses the same mid-latitude channel used in the previous test cases, but the model is initialised by homogeneous zonal winds of 10 m s −1 which are balanced by a surface pressure gradient.
The solution consists of a pair of cyclonic vortices, symmetric on either side of the equator, on the western side of the heating.
These vortices represent the Rossby wave-component of the response as derived in Matsuno (1966).These Rossby waves are accompanied by a (baroclinic) Kelvin wave propagating eastwards along the equator Gill (1980).When comparing our transient solution to those in Gill (1980), it is important to remember that Gill (1980) includes a large amount of damping in order to arrive at a stationary solution.We only apply weak scale-selective damping as described in the storm track test case, and the excited Rossby and Kelvin waves can thus propagate away from the wave source.
After 60 days lead time, the atmosphere achieved a near-equilibrium, which is only slowly evolving in tandem with the evolving ocean.At this time, the coupled solution depends strongly on the chosen parametrisation for the ocean heat transport (Fig. 12).We here again focus the discussion on the simulations using full PE, as the homogeneous density PE solutions are very consistent (cf.Figs. 12 & A3).
With the 0.5-layer ocean, the shape of the initial SST anomaly is entirely intact after 60 days (Fig. 12b), because the ocean itself does not internally redistribute heat.Nevertheless, the ocean mixed layer lost about a quarter of its initial heat anomaly to the atmosphere and the SST anomaly decreased by about 1.2 K.The atmospheric response is much weaker in amplitude than the initial shock-like response to the heating (cf.Figs. 11 & 12a), but it is now following the stationary solution in Gill (1980) more closely, as the transient response has largely dissipated due to numerical and parameterised diffusion.
With the 1.0-layer ocean, direct wind-induced currents and associated ocean heat transports considerably deformed the SST anomaly over 60 days (Fig. 12d).With this parametrisation of the ocean-heat transport, diverging currents do not change the mixed-layer temperature.Atmospheric easterlies induce divergent ocean currents along the equator, thereby increasing the area covered by the warm anomaly without decreasing its amplitude.Overall this leads to a slightly larger ocean heat content with the 1.0-layer compared to the 0.5-layer ocean (compare Fig. 12d and Fig. 12b).Consistent with the slightly larger ocean heat anomaly, the atmospheric response in the 1.0-layer ocean setup is slightly larger in amplitude than for the 0.5-layer ocean (Fig. 12a, c).
Finally, with the 1.25-layer ocean, divergent currents along the equator lead to the upwelling of considerably colder water, leading to large negative SST anomalies along the equator (Fig. 12f).This upwelling is larger in amplitude than the initially positive SST anomaly, which thus largely disappeared after 60 days.In fact, the coldest SSTs appear just eastward of the location of the initial SST anomaly, because the initial response intensified the surface winds in this region.With the warm anomaly largely disappeared, both the SST anomaly field and the atmospheric response is to a first approximation zonally The panel setup for the left column is as in Fig. 11, but to accommodate the much larger sea-level pressure anomalies in (e), black contours are shown at -1.5 hPa through 1.5 hPa in steps of 1 hPa as well as at 5 and 10 hPa, and zonal wind contours are omitted beyond 2.5 m s −1 .
symmetric (Fig. 12e, f).The colder SSTs lead to a marked increase in surface pressure over the equator and a zonally-average meridional flow component away from the equator on either side.
In summary, Bedymo successfully captures both the transient (cf., Matsuno, 1966) and (near-)stationary response (cf., Gill, 1980) to equatorial heating.Beyond these theoretical expectations, the comparison of the different parametrisations of the ocean heat transport demonstrates that both ocean-internal dynamics and air-sea exchange have a profound influence on the We introduced a joint approach to consistently solve the quasi-geostrophic (QG) and two variants of the primitive equations (PE).In all systems, we forecast temperature and surface pressure.In PE, the horizontal wind velocity components also need to be forecasted, whereas all other variables follow diagnostically in QG.Table 1.Summary of the common approach to solve the QG and PE systems.Prognostic equations are marked bold.In this table, u and v denote advecting wind velocities.

Figure 2 .
Figure 2. Baroclinic development in a three-layer channel model on an β-plane for (a,d,g) QG and (b,e,h) homogeneous density PE, and (c,f,i) full PE.From top to bottom, the rows show the development after (a-c) 1 days, (d-f) 3 days, and (g-i) 5 days lead time, respectively.The shading shows temperature in Kelvin and barbs the winds in m s −1 , both representing the lowest of the three layers at σ = 5/6.The contours show surface pressure with a contour interval of 10 hPa centred on 1000 hPa, which contours below 1000 hPa dashed.Note that the not the entire model domain is shown in the y-direction.

Figure 6 .
Figure 6.Sensitivity of the full PE storm track to the strength of temperature relaxation.The panels (a-c) here correspond to panels (a), (b) and (d) in Fig. 5, with the control simulation here being identical to the PE simulation there.Line color :::: colour : and opacity ::::: weight : is consistent across panels.

Figure 7 .Figure 8 .
Figure 7. Developing stationary Rossby wave response to orographic forcing on a β-plane for (a,d,g) QG and (b,e,h) homogeneous density PE, and (c,f,i) full PE.From top to bottom, the rows show the development after (a-c) 1 days, (d-f) 5 days, and (g-i) 9 days lead time, respectively.The shading shows the (barotropic) meridional wind : at :::::::: σ = 0.75 in m s −1 .The contours show surface pressure anomaly with respect to the balanced initial conditions with a contour interval of 4 hPa, showing only the surface pressure anomaly associated with the orography ::by :: the ::: 4 m ::: and :::::::::: 8 m-contours.Note that not the entire model domain is shown in the y-direction.

Figure 11 .
Figure11.Developing wave response in response to an equatorial SST anomaly on a β-plane for full PE.From top to bottom, the rows show the development after (a) 1 day, (b) 3 days, (c) 5 days, and (d) 7 days lead time, respectively.The panel setup is as in Fig.7, but showing wind in the lowest of three levels (σ = 5/6) and with the contour interval for the surface pressure anomaly decreased to 1 hPa.The additional orange-red contours show zonal wind anomalies relative to the initialisation with a contour interval of 1 m s −1 centered around zero and the ±2.5 m s −1 -contour highlighted.

Figure 12 .
Figure 12.Near-stationary wave response to a coupled equatorial SST anomaly on a β-plane for full PE after 60 days of lead time.From top to bottom, the rows show the development of (a : ,b) the 0.5-layer model, (b ::: c,d) the 1-layer model, and (c :: e,f) the 1.25-layer model, respectively.

Figure A1 .Figure A2 .Figure A3 .
Figure A1.As Figure 2, but for simulations in on an f -plane instead of a β-plane.