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

This paper introduces the Bergen Dynamical Model (Bedymo), an idealized atmospheric circulation model, which combines the quasi-geostrophic approximation and the hydrostatic primitive equations into 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 uses σ coordinates in the vertical. In addition to the atmospheric core, Bedymo also includes a slab ocean model with options for prescribed and wind-induced currents. Further, Bedymo has a graphical user interface, making it particularly useful for teaching.

Bedymo is evaluated for four atmosphere-only test cases and one coupled test case including the slab ocean component. The atmosphere-only test cases comprise the growth of a cyclonic disturbance in a baroclinic environment and the excitation of Rossby waves and inertia–gravity waves by isolated orography, as well as the simulation of a mid-latitude storm track, all in a mid-latitude channel. The atmosphere–ocean coupled test case is based on an equatorial channel and evaluates the coupled response to an isolated equatorial temperature anomaly in the ocean mixed layer. For all test cases, results agree well with expectations from theory and results obtained with more complex models.

1 Introduction

Since the 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 idealized 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 – into 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 QGPV (e.g. Charney and Phillips, 1953). In contrast, in PE models the horizontal winds are forecasted independently rather than diagnostically derived from geostrophy. This remains true even if many PE models forecast vorticity and divergence instead of the wind components directly (e.g. Smagorinsky, 1958). As our main goal is to combine the different approximations in one model, we devise a common approach based on prognostic equations for temperature and surface pressure. This approach allows us to solve the respective set of equations as similarly as possible. In contrast to previous studies (e.g. Whitaker, 1993; Rotunno et al., 2000), we can then intercompare QG with PE without using different models.

The simplifications from QG to PE have far-reaching consequences for the representation of the mid-latitude atmospheric dynamics. While in general the synoptic-scale dynamics are relatively well represented in QG (Charney and Phillips, 1953), there are several aspects of synoptic-scale weather systems that QG cannot represent. First and foremost, QG cyclones cannot develop fronts (cf. Hoskins, 1975), which are otherwise the most dynamically active parts of a cyclone. Further, the neglect of advection by the ageostrophic winds makes QG cyclones and anticyclones much more symmetrical than their real counterparts (Wolf and Wirth, 2015). Both fronts and cyclone–anticyclone asymmetries are well represented in a PE model.

From a more conceptual perspective, the difference between QG and PE is the difference between a model that solely represents a balanced component of the flow (QG) and a model that can support flow imbalances (PE). In QG, the model state is given by the distribution of a single variable (i.e. QGPV), all other variables are derived diagnostically from balance relations. In contrast, in PE the horizontal winds and temperature can evolve principally independently from each other. As a consequence, only a PE model can represent the unbalanced flow associated with inertia–gravity waves. The ability to switch between the QG and PE representations makes Bedymo thus ideally suited to address long-standing questions on the relation between the balanced and unbalanced flow components and the role of the unbalanced flow to keep the flow nearly balanced (Plougonven and Snyder, 2007; McIntyre, 2009; Plougonven and Zhang, 2014).

Another reason for starting the model development from scratch was to make use of comparatively new features of Fortran 95 and 2003. These features aid the modularity and readability (and hence maintainability) of the source code. In addition to the combined QG-PE atmospheric model, Bedymo also includes a slab ocean and atmospheric tracers as optional modules. Furthermore, the source code is organized in a way to make it easily accessible from Python. Besides allowing for flexible yet easy-to-read runscripts, the Python bindings provide the basis for a graphical user interface (GUI) that allows us to interactively run the model and watch the flow evolution “live” while the model is running. A user guide is available (Spensberger, 2022b). All these features make Bedymo an ideal tool not only for research but also for education and student research projects.

Table 1 Summary of the common approach to solve the QG and PE systems. Prognostic equations are marked in bold. In this table, u and v denote advecting wind velocities.

2 The model

2.1 Atmospheric dynamics

Our joint approach to solve the equations for QG and for hydrostatic PE is based on the thermodynamic equation, because it features only small modifications between the different approximations. The only difference is the wind velocity components used in the advection (Table 1). In conjunction with a lower boundary condition provided by the forecasted surface pressure, 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 they must be forecasted 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 three-dimensional QG-balance, and σ-vertical velocity is derived from there. 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.

Hoskins et al. (1978)

Table 2 List of symbols. Where ambiguous, superscripts “a” or “o” denote atmospheric or oceanic variables, respectively.

In the following, we summarize the model equations for QG and PE. We use a β-plane approximation for the Coriolis parameter $f=f_{\mathrm{0}}+\mathit{\beta }y$. In order to derive a self-consistent QG system, we need to approximate the specific volumes by time-invariant and horizontally homogeneous background values $\stackrel{\mathrm{‾}}{\mathit{\alpha }}=\stackrel{\mathrm{‾}}{\mathit{\alpha }}\left(\mathit{\sigma }\right)$. This approximation is optional for PE, and we will introduce and evaluate both variants. All symbols used in the following equations are summarized in Table 2.

2.1.1 Full primitive equations

The PE equations, as used in Bedymo, are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}-(f_{\mathrm{0}}+\mathit{\beta }y)v=-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}-\mathit{\alpha }\mathit{\sigma }{\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial x}}-ru+D{\mathrm{\nabla }}^{\mathrm{2}}u,\text{(2)}& {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}}+(f_{\mathrm{0}}+\mathit{\beta }y)u=-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial y}}-\mathit{\alpha }\mathit{\sigma }{\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial y}}-rv+D{\mathrm{\nabla }}^{\mathrm{2}}v,\text{(3)}& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{RT}{\mathit{\sigma }{p}_{\mathrm{s}}}}\text{(4)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\sigma }}}=-{\displaystyle \frac{RT}{\mathit{\sigma }}},\text{(5)}& {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}-{\displaystyle \frac{\mathit{\alpha }}{c_p}}\mathit{\omega }={\displaystyle \frac{J}{c_p}},\text{(6)}& {\displaystyle \frac{\mathrm{1}}{p_{\mathrm{s}}}}{\displaystyle \frac{\mathrm{d}{p}_{\mathrm{s}}}{\mathrm{d}t}}+{\mathrm{\nabla }}_{\mathit{\sigma }}\cdot \mathit{v}+{\displaystyle \frac{\partial \dot{\mathit{\sigma }}}{\partial \mathit{\sigma }}}=\mathrm{0}.\end{array}$$

In order to provide an energy sink for long and short waves, respectively, the momentum equations include non-scale-selective linear Ekman friction and a scale-selective damping term.

In this system, we infer geopotential ϕ by integrating the hydrostatic Eq. (4) upwards, starting from the time-invariant surface geopotential ϕ_s which represents the model orography. The local surface pressure tendency results from continuity,

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial t}}=-\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathrm{\nabla }\cdot \left(\mathit{v}{p}_{\mathrm{s}}\right)\mathrm{d}\mathit{\sigma }.\end{array}$$

And finally,

$$\begin{array}{}\text{(8)}& {\displaystyle }\dot{\mathit{\sigma }}& {\displaystyle }=-{\displaystyle \frac{\mathit{\sigma }}{p_{\mathrm{s}}}}{\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial t}}-{\displaystyle \frac{\mathrm{1}}{p_{\mathrm{s}}}}\underset{\mathrm{0}}{\overset{\mathit{\sigma }}{\int }}\mathrm{\nabla }\cdot \left(\mathit{v}{p}_{\mathrm{s}}\right)\mathrm{d}\mathit{\sigma }\text{, and}\text{(9)}& {\displaystyle }\mathit{\omega }& {\displaystyle }=-\underset{\mathrm{0}}{\overset{\mathit{\sigma }}{\int }}\mathrm{\nabla }\cdot \left(\mathit{v}{p}_{\mathrm{s}}\right)\mathrm{d}\mathit{\sigma }+\mathit{\sigma }\mathit{v}\cdot \mathrm{\nabla }{p}_{\mathrm{s}},\end{array}$$

which closes the system of Eqs. (1)–(6). All other parameters in the equations are either time-invariant fields, are constants, or represent an external forcing.

2.1.2 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 $\mathit{\alpha }\approx \stackrel{\mathrm{‾}}{\mathit{\alpha }}\left(\mathit{\sigma }\right)$ 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 $\stackrel{\mathrm{‾}}{\mathit{\alpha }}\mathit{\sigma }$ only varies with height. In fact, using the ideal gas law $\stackrel{\mathrm{‾}}{\mathit{\alpha }}\mathit{\sigma }=\frac{R\stackrel{\mathrm{‾}}{T}}{p_{\mathrm{s}}}$, the product $\stackrel{\mathrm{‾}}{\mathit{\alpha }}\mathit{\sigma }$ becomes constant with an isothermal background state. Further, vertical temperature advection is approximated by advection of the homogeneous background state only, $\dot{\mathit{\sigma }}\frac{\partial T}{\partial \mathit{\sigma }}\approx \dot{\mathit{\sigma }}\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}$ with $\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}=\frac{\partial \stackrel{\mathrm{‾}}{T}}{\partial \mathit{\sigma }}$.

In summary, the resulting momentum and thermodynamic equations are

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}-(f_{\mathrm{0}}+\mathit{\beta }y)v& {\displaystyle }=-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\mathit{\sigma }{\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial x}}-ru+D{\mathrm{\nabla }}^{\mathrm{2}}u,\text{(11)}& {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}}+(f_{\mathrm{0}}+\mathit{\beta }y)u& {\displaystyle }=-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial y}}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\mathit{\sigma }{\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial y}}-rv+D{\mathrm{\nabla }}^{\mathrm{2}}v,\text{(12)}& {\displaystyle \frac{d_{h}T}{\mathrm{d}t}}+\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}\dot{\mathit{\sigma }}-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}{c_p}}\mathit{\omega }& {\displaystyle }={\displaystyle \frac{J}{c_p}}.\end{array}$$

Continuity (Eq. 6), hydrostasy (Eq. 4), and Eqs. (7)–(9) that determine the surface pressure tendency and the vertical velocities all remain unchanged.

2.1.3 Quasi-geostrophy

In QG, the geostrophic wind follows from the geostrophic streamfunction ψ_g ,

$$\begin{array}{}\text{(13)}& {\mathit{v}}_{\mathrm{g}}=-\mathit{k}\times \mathrm{\nabla }{\mathit{\psi }}_{\mathrm{g}},\end{array}$$

which is in turn defined by

$$\begin{array}{}\text{(14)}& f_{\mathrm{0}}{\mathit{\psi }}_{\mathrm{g}}=\mathit{\varphi }+\stackrel{\mathrm{‾}}{\mathit{\alpha }}\mathit{\sigma }{p}_{\mathrm{s}}.\end{array}$$

The prognostic equations are

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{d_{\mathrm{g}}T}{\mathrm{d}t}}+\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}\dot{\mathit{\sigma }}-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}{c_p}}\mathit{\omega }={\displaystyle \frac{J}{c_p}}\text{and}\text{(16)}& {\displaystyle \frac{d_{\mathrm{g}}{p}_{\mathrm{s}}}{\mathrm{d}t}}={\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial t}}+{\mathit{v}}_{\mathrm{gs}}\cdot \mathrm{\nabla }{p}_{\mathrm{s}}={\mathit{\omega }}_{\mathrm{s}},\end{array}$$

in which the surface pressure tendency Eq. (16) is evaluated at the lower surface (σ=1), indicated by the subscript “s”.

The remaining unknown variables at this point are the vertical wind $\dot{\mathit{\sigma }}$ and the pressure vertical velocity ω. Starting from the following form of continuity

$$\begin{array}{}\text{(17)}& {\mathrm{\nabla }}_{\mathit{\sigma }}\cdot {\mathit{v}}_{\mathrm{a}}+{\displaystyle \frac{\mathrm{1}}{p_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathit{\omega }}{\partial \mathit{\sigma }}}=\mathrm{0},\end{array}$$

the pressure vertical velocity can be derived from the ω equation

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{s}{\mathrm{\nabla }}_{\mathit{\sigma }}^{\mathrm{2}}\mathit{\omega }+{\displaystyle \frac{f_{\mathrm{0}}^{\mathrm{2}}}{p_{\mathrm{0}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\omega }}{\partial {\mathit{\sigma }}^{\mathrm{2}}}}& =\mathrm{2}{\mathrm{\nabla }}_{\mathit{\sigma }}\cdot (Q_{\mathrm{1}},Q_{\mathrm{2}})+\mathit{\beta }{f}_{\mathrm{0}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial \mathit{\sigma }}}\\ & +f_{\mathrm{0}}{\displaystyle \frac{\partial r{\mathit{\zeta }}_{\mathrm{g}}}{\partial \mathit{\sigma }}}-f_{\mathrm{0}}{\displaystyle \frac{\partial D{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{g}}}{\partial \mathit{\sigma }}}-{\displaystyle \frac{\mathit{\kappa }}{\mathit{\sigma }}}{\mathrm{\nabla }}^{\mathrm{2}}J,\end{array}\end{array}$$

in which Q₁ and Q₂ represent the two components of the Hoskins et al. (1978) Q vector

$$\begin{array}{}\text{(19)}& {\displaystyle }{Q}_{\mathrm{1}}=+& {\displaystyle }{f}_{\mathrm{0}}{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{g}}}{\partial x}}\cdot {\mathrm{\nabla }}_{\mathit{\sigma }}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{g}}}{\partial \mathit{\sigma }}},\text{(20)}& {\displaystyle }{Q}_{\mathrm{2}}=-& {\displaystyle }{f}_{\mathrm{0}}{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{g}}}{\partial y}}\cdot {\mathrm{\nabla }}_{\mathit{\sigma }}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{g}}}{\partial \mathit{\sigma }}}.\end{array}$$

Further, $\stackrel{\mathrm{‾}}{s}=\stackrel{\mathrm{‾}}{s}\left(\mathit{\sigma }\right)$ is a horizontally homogeneous stability parameter

$$\begin{array}{}\text{(21)}& \stackrel{\mathrm{‾}}{s}=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}{\mathit{\sigma }}}{\displaystyle \frac{c_v}{c_p}}-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\alpha }}}{\partial \mathit{\sigma }}}.\end{array}$$

In order to solve the elliptic ω equation (Eq. 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

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{s}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\omega }}_{\mathrm{s}}-{\displaystyle \frac{f_{\mathrm{0}}^{\mathrm{2}}}{p_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathit{\omega }}{\partial \mathit{\sigma }}}{|}_{\mathrm{s}}& =-{\mathit{v}}_{\mathrm{gs}}\cdot \mathrm{\nabla }\left({\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{s}}\right)\\ & -f_{\mathrm{0}}\mathit{\beta }{v}_{\mathrm{gs}}-f_{\mathrm{0}}{r}_{\mathrm{s}}{\mathit{\zeta }}_{\mathrm{gs}}+f_{\mathrm{0}}{D}_{\mathrm{s}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{gs}},\end{array}\end{array}$$

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 $\dot{\mathit{\sigma }}$ from

$$\begin{array}{}\text{(23)}& \dot{\mathit{\sigma }}={\displaystyle \frac{\mathrm{1}}{p_{\mathrm{0}}}}(\mathit{\omega }-\mathit{\sigma }{\mathit{\omega }}_{\mathrm{s}}-\mathit{\sigma }({\mathit{v}}_{\mathrm{g}}-{\mathit{v}}_{\mathrm{gs}})\cdot \mathrm{\nabla }{p}_{\mathrm{s}}),\end{array}$$

which closes the system of equations.

QG potential vorticity q is not used anywhere in the model. Nevertheless, for reference, q takes the form

$$\begin{array}{}\text{(24)}& q={\mathrm{\nabla }}_{\mathit{\sigma }}^{\mathrm{2}}{\mathit{\psi }}_{\mathrm{g}}+\mathit{\beta }y+{\displaystyle \frac{f_{\mathrm{0}}^{\mathrm{2}}}{p_{\mathrm{0}}}}{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left[{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{s}}}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{g}}}{\partial \mathit{\sigma }}}\right],\end{array}$$

in this QG system.

2.2 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 timescales on which the internal ocean dynamics can be neglected. Nevertheless, as the oceanic heat transport might play a role even at these timescales, 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,

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\partial T^{\mathrm{o}}}{\partial t}}=-{\displaystyle \frac{F_{\mathrm{sh}}}{{\mathit{\rho }}^{\mathrm{o}}{c}_p^{\mathrm{o}}H}}-\text{OHT}-r_{\mathrm{T}}^{\mathrm{o}}(T^{\mathrm{o}}-T_{\mathrm{e}}^{\mathrm{o}}).\end{array}$$

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_{\mathrm{e}}^{\mathrm{o}}$, which may be used to crudely represent the neglected oceanic circulation (see Table 2 for other symbols).

The heat exchange is parameterized with a bulk flux formulation,

$$\begin{array}{}\text{(26)}& F_{\mathrm{sh}}={\mathit{\rho }}_{\mathrm{s}}^{\mathrm{a}}{c}_p^{\mathrm{a}}{C}_{\mathrm{sh}}\left|{\mathit{v}}_{\mathbf{s}}^{\mathbf{a}}\right|(T_{\mathrm{s}}^{\mathrm{a}}-T^{\mathrm{o}}).\end{array}$$

Here (and in the following equations) 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 Sect. 2.4.

The options for parametrizing the oceanic heat transports are

$$\begin{array}{}\text{(27)}& \text{OHT}=\left\{\begin{array}{ll}\mathrm{0}& \text{for a 0.5-layer model,}\\ {\mathit{v}}^{\mathrm{o}}\cdot {\mathrm{\nabla }}_h{T}^{\mathrm{o}}& \text{for a 1-layer model, and}\\ \mathrm{\nabla }\cdot \left({\mathit{u}}^{\mathrm{o}}{T}^{\mathrm{o}}\right)& \text{for a 1.25-layer model.}\end{array}\right.\end{array}$$

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):

$$w^{\mathrm{o}}={\mathrm{\nabla }}_h\cdot {\mathit{v}}^{\mathbf{o}},$$

that transports water from a lower layer of temperature $T_{\mathrm{2}}^{\mathrm{o}}$ into the mixed layer,

$$\mathrm{\nabla }\cdot \left({\mathit{u}}^{\mathrm{o}}{T}^{\mathrm{o}}\right)={\mathrm{\nabla }}_h\cdot \left({\mathit{v}}^{\mathrm{o}}{T}^{\mathrm{o}}\right)-\left\{\begin{array}{ll}{w}^{\mathrm{o}}{T}^{\mathrm{o}}& \text{if downwelling,}{w}^{\mathrm{o}}<\mathrm{0},\\ w^{\mathrm{o}}{T}_{\mathrm{2}}^{\mathrm{o}}& \text{if upwelling,}{w}^{\mathrm{o}}\ge \mathrm{0}.\end{array}\right.$$

The temperature $T_{\mathrm{2}}^{\mathrm{o}}$ might vary spatially but is kept constant over time. By not letting $T_{\mathrm{2}}^{\mathrm{o}}$ vary with time, we implicitly assume the deeper layer to be motionless and to have infinite heat capacity.

Figure 1 Screenshot of the graphical user interface for Bedymo taken while the model is running.

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.25-layer 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_{\mathrm{2}}^{\mathrm{o}}=T^{\mathrm{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 ${\mathit{v}}^{\mathbf{o}}={\mathit{v}}_{\mathbf{e}}+{\mathit{v}}_{\mathrm{p}}$ is a combination of a user-prescribed flow pattern v_p and the wind-driven Ekman flow v_e. Our formulation of v_e also includes a small parameter ϵ in addition to the standard Ekman flow,

$$\begin{array}{}\text{(28)}& {\mathit{v}}_{\mathrm{e}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}^{\mathrm{a}}{C}_D\left|{\mathit{v}}^{{\mathrm{a}}_{\mathrm{s}}}\right|}{{\mathit{\rho }}^{\mathrm{o}}H({\mathit{ϵ}}^{\mathrm{2}}+f^{\mathrm{2}})}}(\mathit{ϵ}{\mathit{v}}^{{\mathrm{a}}_{\mathrm{s}}}-f\mathit{k}\times {\mathit{v}}^{{\mathrm{a}}_{\mathrm{s}}}),\end{array}$$

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.

2.3 Numerics and model infrastructure

2.3.1 Coordinate system and discretization

Bedymo uses Cartesian coordinates in the horizontal and a terrain-following pressure coordinate σ in the vertical, which is defined by

$$\begin{array}{}\text{(29)}& \mathit{\sigma }={\displaystyle \frac{p}{p_{\mathrm{s}}}}.\end{array}$$

Variations in surface pressure thus affect the coordinate system throughout the atmospheric column up to the model top at p=0.

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.

Table 3 Summary of pertinent model parameters for the test setups used to evaluate the model. The first four columns represent the four mid-latitude test cases discussed in Sect. 3.1, with IGW being short for inertia–gravity wave. The fifth column represents the coupled test case, which is the Matsuno–Gill-like response to equatorial heating discussed in Sect. 3.2.

Bedymo uses the third-order Runge–Kutta time integration scheme in conjunction with a third-order upwind-biased interpolation of the advected quantities to determine the advective fluxes at the staggered locations of the wind speed components, following the derivations of Smolarkiewicz (1982) and Tremback et al. (1987). A first-order upwind-biased interpolation is only used for vertical advection in the uppermost and lowermost level. We tested higher- and lower-order advection schemes for the interior of the model domain, but this choice turned out to be largely inconsequential. The explicit scale-selective damping active for most of the test setups dominates the numerical diffusion (see Table 3). The remainder of the model is discretized using second-order centred differences and interpolations.

2.3.2 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 stabilized (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.

2.4 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, i.e. $\dot{\mathit{\sigma }}=\mathrm{0}$ and ω=0. At the surface $\dot{\mathit{\sigma }}=\mathrm{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 boundary-normal fluxes.

2.5 Python bindings and graphical user interface

While the model is written in Fortran, the code is structured such to make it easily accessible from python using F2PY (NumPy Developers, 2022). In particular, the model can be run from python and model variables are accessible from python during runtime.

These python bindings are the basis for a rudimentary python run script illustrating the basic use of the model through python, as well as a graphical user interface (GUI) based on python and PyQt (Fig. 1). The GUI allows us to start, stop, and single-step the model and watch the evolution of all prognostic and pertinent diagnostic variables interactively while the model is running. Both the run script and the GUI are included in the Bedymo source code repository.

Figure 2 Baroclinic development in a three-layer channel model on an β plane for (a, d, g) QG, (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 d, (d–f) 3 d, and (g–i) 5 d lead time, respectively. The shading shows temperature in Kelvin, and the barbs show the winds (in m s⁻¹), both representing the lowest of the three layers at $\mathit{\sigma }=\mathrm{5}/\mathrm{6}$. The contours show surface pressure with a contour interval of 10 hPa centred on 1000 hPa, with contours below 1000 hPa shown as dashed lines. Note that the entire model domain is not shown in the y direction.

3 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 four test cases focus on important aspects of the mid-latitude atmospheric dynamics in isolation. In a 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 Table 3.

3.1 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 zonal-mean statistics of the storm track simulated by Bedymo. Finally, test cases three and four evaluate the development towards the stationary wave responses to isolated orography, with a focus on Rossby waves and inertia–gravity waves, respectively.

3.1.1 Baroclinic cyclogenesis on an f and β-plane

The baroclinic channel for the cyclogenesis test case is 16 000 km long and periodic in the zonal direction. In the meridional direction it is 12 000 km wide and bounded by impermeable walls. The horizontal resolution is 100 km, and we use three levels in the vertical. We initialize a baroclinic zone with a temperature contrast of in total about 32 K distributed over about 3000 km. Temperature stratification is determined by $\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}=\mathrm{0.8}{\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}}_{\mathrm{d}}$, with ${\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}}_{\mathrm{d}}$ the dry adiabatic lapse rate in σ coordinates. This approximately corresponds to an overall Brunt–Väisälä frequency $N^{\mathrm{2}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻². The baroclinic zone is meridionally centred in the channel. There is no initial surface pressure gradient, and thus there are 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⁻¹ in the uppermost layer ($\mathit{\sigma }=\mathrm{1}/\mathrm{6}$).

To localize the baroclinic development, we perturb this initial state by a warm anomaly of 2 K located at the centre 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 or 70 m s⁻¹, respectively.

The baroclinic development for the control setup and the three model variants is summarized in Fig. 2. For the first 3 d, the evolution in the three model variants is quite similar (Fig. 2a–f), and only at day 5 do 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.

Figure 3 Growth rate of eddy kinetic energy averaged over 24 h periods for QG (orange) and the two variants of PE (blue, light blue). Line weight indicates the initial and jet strength with maximum winds of 30 m s⁻¹ (thin), 50 m s⁻¹ (medium), and 70 m s⁻¹ (bold).

Comparing the two PE variants, only minor differences occur up to day 5 (Fig. 2h and 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 are comparatively cooler than their counterparts in full PE. Given that the only differences between the two PE variants is linearized vertical advection and a linearized pressure gradient force, it seems implausible that these differences in temperature originate from the linearized 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 h (QG), 19.0 h (full PE), and 21.1 h (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

$${\mathit{\sigma }}_{\mathrm{E}}=\mathrm{0.31}{\displaystyle \frac{u_{\max }}{L_{\mathrm{d}}}}\approx \mathrm{1.66}\times {\mathrm{10}}^{-\mathrm{5}}{\text{s}}^{-\mathrm{1}},$$

for the control jet intensity of 50 m s⁻¹. This corresponds to an e-folding timescale of τ_E=16.7 h, which fits quite well, particularly 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⁻¹, the Eady growth timescale reduces to 11.9 h. This fits almost perfectly to the simulated maximum intensification with a timescale of 11.7 h 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⁻¹, the correspondence between theory and simulation is less tight, with an Eady growth timescale of 27.8 h and a QG maximum intensification with a time scale of 35.6 h. This might at least partly be due to the β effect and downstream development playing a comparatively larger role with weaker baroclinicity. But this cannot be the full explanation because for full PE the reduction in the peak growth rate (about 44 %) is only slightly larger than the reduction in baroclinicity (40 %).

Figure 4 Sensitivity of the downstream development to the initial jet speed in a three-layer channel model on a β plane using jet speeds of (a–c) 30 m s⁻¹, (d–f) 50 m s⁻¹, and (g–i) 70 m s⁻¹. The presented lead times are scaled linearly with the initialized jet speed, such that the rows show the development after (a–c) 8.25 d, (d–f) 5 d, and (g–i) 3.5 d lead time, respectively. Columns show the different versions of the model: panels (a, d, g) show QG, panels (b, e, h) show homogeneous-density PE, and panels (c, f, i) show full PE. Barbs and contours are the same as in Fig. 2, and panels (d–f) are identical to panels (g–i) of Fig. 2 except for the temperature scale.

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 normalized 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 and 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 centre 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 idealized cyclones with more complex models (e.g. Schemm et al., 2013; Terpstra and Spengler, 2015, and the cyclones in the colder environments of Tierney et al., 2018). This consistency is not surprising given the much earlier results of Simmons and Hoskins (1978), for instance, which used a model similar in complexity to Bedymo to study cyclogenesis. However, 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 parameterizations.

3.1.2 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 parameterizations for three physical processes to the model setup. First, we enable a temperature relaxation towards the initial state with a timescale of 10⁶ s≈11.6 d throughout the model domain to represent all processes that replenish baroclinicity in real storm tracks. Second, we enable linear Ekman friction in the lowest model layer with r=10⁶ s. Finally, we enable bi-harmonic diffusion that damps 2Δx waves with a timescale of 10⁵ s≈28 h. This timescale 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 d each, of which the first 2 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.

Figure 5 Statistical equilibrium state of the QG and full PE storm tracks, respectively. Line colour and weight are consistent throughout all the panels and identical for panels (a), (b), and (d). Throughout this figure, blue lines represent full PE, and orange represent lines QG. Line weight indicates the vertical level, with the boldest lines representing the uppermost level.

Because of the absence of any zonal asymmetries, the resulting storm track 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 $\frac{\partial }{\partial y}\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}<\mathrm{0}$ indicates a convergence of eddy momentum fluxes in the time mean, for example, analogous to $\frac{\partial }{\partial y}\stackrel{\mathrm{‾}}{v^{\prime }{T}^{\prime }}<\mathrm{0}$ and $\frac{\partial }{\partial \mathit{\sigma }}\stackrel{\mathrm{‾}}{{\dot{\mathit{\sigma }}}^{\prime }{T}^{\prime }}<\mathrm{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). Maximum zonal winds in the uppermost level are just below 50 m s⁻¹, which would be the wind speed required to thermally balance the initial and relaxation temperature gradient in the absence of near-surface winds. Consistent with the mean zonal winds, the mean convergence of momentum is symmetric around the imposed baroclinic zone with maximum momentum convergence at the jet core, which is also strongest at the jet level (Fig. 5b). In contrast, heat transports are strongest in the lowest level (Fig. 5d and e), consistent with the largest deviation of mean temperature from the relaxation state (Fig. 5c).

The described QG storm track is qualitatively similar to the storm tracks produced by other models (e.g. Held and Suarez, 1994; Vallis et al., 2004; Frierson et al., 2006; Voigt and Shaw, 2016) and the ones actually observed (e.g. Figs. 7.7 and 11.8 of Peixoto and Oort, 1992). While these qualitative similarities are encouraging, the evaluation of Bedymo through comparison with existing models is somewhat limited by important differences between these models. Aquaplanet and Held–Suarez-type model simulations are global and thus both include a Hadley circulation and take into account the spherical geometry of the Earth. Both of these features are missing in Bedymo but will certainly affect the simulated storm track in these other models. For example, by construction, the storm track and jets in Bedymo are purely eddy-driven, while Held–Suarez-type simulations also include a thermally driven subtropical jet.

In contrast to QG, the simulated PE storm track is not symmetric around the imposed baroclinic zone. The core of the jet 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 are 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, with spherical geometry and the Hadley circulation again likely affecting the dynamics. In fact, some setups of the spherical one-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 do not, however, affect the PE storm track, as the meridional profiles of all parameters remain nearly unchanged in a simulation with a meridionally wider channel.

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 colour and weight is consistent across panels.

Finally, the simulated storm track remains largely consistent across a range of temperature relaxation timescales (Fig. 6). In addition to the default relaxation with a timescale of 11.6 d, we conducted simulations with full PE using timescales of 3.47 and 34.7 d. These simulations are called “strong” and “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 and stronger jets (Fig. 6).

3.1.3 Rossby waves excited by orography

The third 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 initialized by homogeneous zonal winds of 10 m s⁻¹, which are balanced by a surface pressure gradient. In order to better resolve the vertical structure of the Rossby wave pattern, we use 10 vertical levels compared to 3 for the storm track test cases. Orography is introduced as an isolated elliptical cone-shaped mountain, meridionally centred in the domain, with a meridional scale of 1500 km, a zonal scale of 1000 km, and a height of 100 m² s⁻².

Finally, we increase the scale-selective damping parameter D to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$m² s⁻¹ to be able to damp out gravity waves from slight initial imbalances in the PE solutions. These imbalances are zonally symmetric and thus appear as zonally symmetric anomalies in the meridional wind on top of the Rossby wave pattern. Such anomalies are visible at day 1 and 5 (Fig. 7b, c, e, f) but have largely disappeared at day 9 (Fig. 7h and i).

Figure 7 Developing stationary Rossby wave response to orographic forcing on a β plane for (a, d, g) QG, (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 d, (d–f) 5 d, and (g–i) 9 d lead time. The shading shows the meridional wind at σ=0.75 in m s⁻¹. The contours show the orography by the 4 and 8 m contours. Note that the entire model domain is not shown in the y direction.

Figure 7 shows the developing Rossby wave train excited by the mountain for the three model variants. As expected from theory, the QG response is meridionally symmetric. In contrast, the two variants of the PE again show slight meridional asymmetries. These asymmetries become increasingly more pronounced with increasing mountain height (not shown) and are thus most likely due to non-linear interactions of the perturbation flow with the mountain and itself. But despite the slight asymmetries in the PE solutions, the response is very consistent in both shape and amplitude across the three model variants.

Figure 8 Linear quasi-geostrophic solution of the stationary Rossby wave response to isolated orography. Shading and contours are the same as in Fig. 7 except that the meridional wind is shown at z=2500 m.

Further, the response fits qualitatively well to the response expected from linear wave theory (Fig. 8; Hoskins and Karoly, 1981; Held et al., 2002). The horizontal wave pattern in particular compares very well between the analytic model and Bedymo, whereas there is about a factor 2 and 4 difference in amplitude compared with the PE and the QG solutions, respectively (Figs. 7 and 8). For reasonably low values of the damping parameter (corresponding to r in Bedymo), the amplitude in the analytic model is first and foremost set by the lower boundary condition, which is markedly different from both the QG and the PE modes in Bedymo. In the analytic model, the lower boundary is given by a temperature anomaly at z=0, whereas in Bedymo it is given by a surface geopotential at σ=0. The surface geopotential enters the QG system by modifying the surface vorticity field (cf. Eq. 22), whereas it enters the PE system through the continuity equation. Thus, only in the PE mode of Bedymo does the mountain actually block any volume. Some discrepancies with respect to the amplitude of the wave pattern are thus to be expected.

3.1.4 Inertia–gravity waves excited by orography

The fourth and final mid-latitude test case evaluates the developing stationary inertia–gravity wave (IGW) pattern response to isolated orography. In comparison to the previous Rossby wave test case, all horizontal scales, including the extent of the mountain and the grid spacing, are decreased by an order of magnitude. We further increased the number of vertical levels to 30 in order to properly resolve the vertical propagation of the IGWs. To isolate the IGW from the Rossby wave response, we approximate the Coriolis effect by an f plane. QG cannot represent IGWs, and the two variants of PE yield nearly indistinguishable results. We thus focus our presentation on the results from full PE.

Figure 9 Developing stationary inertia–gravity wave response to orographic forcing on an f plane for full PE. The left column (a, c) shows the horizontal wave pattern at σ=0.7. The right column (b, d) shows the vertical wave propagation in the x–z plane at y≈0. From top to bottom, the rows show the development after (a, b) 12 h and (c, d) 48 h. The shading shows the pressure vertical wind (in Pa s⁻¹), the contours show the orography by the 4 and 8 m contours. Note that the entire model domain is not shown in the y direction.

Figure 10 Linear stationary inertia–gravity wave response to the forcing from isolated orography. Shading and contours are the same as in each row of Fig. 9 except that the shading shows the Cartesian vertical wind (in mm s⁻¹).

Again we compare the transient solutions from Bedymo (Fig. 9) with the linear stationary solution (Fig. 10). Both the horizontal and the vertical wave propagation patterns qualitatively fit very well. Only the wave amplitude is again about a factor of 2 larger in the Bedymo-PE solution compared to the linear stationary one. Further, 48 h into the simulation, the vertical wave pattern shows first signs of interference with waves reflected from the model top (Fig. 9d), whereas the horizontal propagation of the inertia–gravity waves is still largely undisturbed by the model lateral boundaries (Fig. 9c).

3.2 Coupled tests

In addition to the mid-latitude test case, we subject Bedymo to a test case centred on tropical air–sea interaction. The model domain remains a zonally periodic channel with a meridional width of 10 000 km. Zonally, the channel is extended to cover the actual circumference of the Earth of 40 000 km and the Coriolis parameter is adapted to represent an equatorial β plane. As the model domain contains the Equator, we use only the two PE variants of Bedymo for this test case.

Initial surface air temperatures and sea surface temperatures (SSTs) are homogeneous at 288 K except for a +5 K SST anomaly with a radius of 1500 km centred on the Equator. Below this warm mixed layer is a 15 K colder deeper-ocean layer without heat anomaly. In case of the 1.25-layer ocean model, this cold deep water can be upwelled to affect the ocean mixed-layer temperature and thus the SST. The atmosphere is initialized with homogeneous and barotropic easterlies with a wind speed of 10 m s⁻¹ that is balanced as far as possible and necessary by a meridional surface pressure gradient. Temperature stratification in the atmosphere is as in the cyclogenesis test case, corresponding approximately to $N^{\mathrm{2}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻². Neither ocean nor atmospheric temperatures are relaxed.

Figure 11 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 d, (b) 3 d, (c) 5 d, and (d) 7 d lead time, respectively. The panel setup is the same as in Fig. 7 but showing wind in the lowest of three levels ($\mathit{\sigma }=\mathrm{5}/\mathrm{6}$) and with the contour interval for the surface pressure anomaly decreased to 1 hPa. The additional orange and red contours show zonal wind anomalies relative to the initialization with a contour interval of 1 m s⁻¹ centred around zero and the ±2.5 m s⁻¹ contour highlighted.

The initial transient atmospheric response to the surface heating is shown in Fig. 11. As the ocean heat transport only has a minor influence within the initial 7 d shown in Fig. 11 and homogeneous-density PE and full PE yield very similar results (cf. Figs. 11 and A2), here we focus on the combination of full PE with the 0.5-layer ocean that ignores ocean heat transport (Eq. 27).

Although the forcing from the +5 K SST anomaly is clearly finite in amplitude, the initial atmospheric response to equatorial surface heating is consistent with the linear Matsuno–Gill response derived in Matsuno (1966) and Gill (1980). 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.

Figure 12 Near-stationary wave response to a coupled equatorial SST anomaly on a β plane for full PE after 60 d of lead time. From top to bottom, the rows show the development of (a, b) the 0.5-layer model, (c, d) the 1-layer model, and (e, f) the 1.25-layer model. The panel setup for the left column is the same as in Fig. 11, but to accommodate the much larger sea level pressure anomalies in (e), black contours are shown at −1.5 through 1.5 hPa in steps of 1 hPa and at 5 and 10 hPa, and zonal wind contours are omitted beyond 2.5 m s⁻¹. The maximum zonal wind anomaly on either side of the Equator is about 7 m s⁻¹. The right column shows ocean mixed-layer temperatures in Kelvin, with the numbers in the top-right corners indicating the domain-integral mixed-layer heat anomaly relative to 288 K.

After 60 d 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 parametrization 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 and A3).

With the 0.5-layer ocean, the shape of the initial SST anomaly is entirely intact after 60 d (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 and 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 parameterized diffusion.

With the 1.0-layer ocean, direct wind-induced currents and associated ocean heat transports considerably deformed the SST anomaly over 60 d (Fig. 12d). With this parametrization 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 ocean compared to the 0.5-layer ocean (compare Fig. 12d and b). 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 and 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 d. 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 symmetric (Fig. 12e and f). The colder SSTs lead to a marked increase in surface pressure over the Equator and a zonally averaged 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 parameterizations of the ocean heat transport demonstrates that both ocean internal dynamics and air–sea exchange have a profound influence on the coupled solution on timescales longer than a week or so. The atmosphere–ocean interactions occurring in the 1.0-layer and 1.25-layer ocean are however beyond what could be captured by analytic solutions, such that the work of Matsuno (1966) and Gill (1980) can only provide limited guidance. A further caveat with this conclusion is that Bedymo so far cannot represent dynamical balances in the ocean, and thus cannot represent the ocean gyre circulations in the mid-latitudes, for example.

4 Conclusions

We introduced a joint approach to consistently solve the quasi-geostrophic (QG) equation 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.

We implemented this approach in the Bergen Dynamic Model (Bedymo) and demonstrated the feasibility of the approach and the performance of Bedymo on the basis of five test cases. These cases are (a) the baroclinic development of a cyclonic disturbance, (b) the representation of mid-latitude storm tracks, (c, d) the excitation of Rossby and inertia–gravity waves by isolated orography, and (e) the coupled response of the PE models with a slab ocean to an equatorial temperature anomaly in the ocean mixed layer. In all cases, the model results agree well with either an analytical solution for the corresponding linearized problem or conceptual models.

By successfully combining QG and PE into one consistent model, Bedymo considerably simplifies the comparison of the dynamical differences between QG and PE because it eliminates all error sources associated with comparing two different models. The ability to simply switch between the approximations is especially valuable in cases where the formal validity of QG becomes questionable. One such example is the treatment of orography because the assumptions underlying QG formally require orographic slopes to be negligibly small.

5 Outlook

The development of Bedymo will not stop with this publication. We foresee two main avenues for future development. First, we plan to include a parametrization of moist diabatic effects. Bedymo already contains the technical basis for such a parametrization in the form of a passive tracer module and the model infrastructure allowing the tracer module to exert diabatic forcing on the dry dynamical core (not evaluated here). Second, we plan to complement the current Cartesian geometry in the horizontal using an option that takes the spherical geometry of the Earth into account. In a longer-term perspective, we envision including an option for semi-geostrophy (Hoskins, 1975) as an intermediate step between QG and PE, as well as an option for a more dynamically active ocean component.

However, already in its current state at the time of publication, Bedymo represents a unique research tool. It represents an idealized complement to general circulation models run in Held–Suarez or Aquaplanet simulations. Compared to these models, Bedymo omits moist diabatic effects, spherical geometry, and the Hadley circulation, which is advantageous for studying all phenomena for which these aspects of the mid-latitude dynamics should not play a role. Further, due to the ease of intercomparisons between the QG and PE approximations, Bedymo is ideally suited to assess the still-unclear relation between the balanced and the unbalanced components of mid-latitude flow (cf. Plougonven and Zhang, 2014).

Appendix A

Figure A1 The same as Fig. 2 but for simulations on an f plane instead of a β plane.

Figure A2 The same as Fig. 11 but for homogeneous-density PE instead of full PE.

Figure A3 The same as Fig. 12 but for homogeneous-density PE instead of full PE. As in Fig. 12, the maximum zonal wind anomaly on either side of the Equator is about 7 m s⁻¹.