: Rotating flows, atmospheric dynamics, Discontinous Galerkin (DG) method, Non-conforming mesh adaptivity, IMEX Runge-Kutta method, compressible Euler equations.

1 Introduction↩︎

Atmospheric flows are characterized by the interaction of multiple spatial and temporal scales, ranging from large-scale balanced circulations to fast acoustic and gravity waves [1], [2]. In fully compressible atmospheric models, the coexistence of slow geostrophic modes and advective dynamics as well as fast acoustic and gravity modes introduces severe stiffness, strongly constraining the time step size used in fully explicit time integration methods. Efficient numerical discretisations must therefore accurately capture the meteorologically relevant dynamics while reducing the impact of the acoustic stability constraint.

To address the stiffness arising from fast acoustic waves, implicit-explicit (IMEX) time integration methods are particularly attractive for compressible atmospheric flows [3]–[5]. In IMEX formulations, the stiff acoustic contributions are treated implicitly, while the slower advective dynamics are advanced explicitly [6], [7], allowing the time step to be determined by the meteorologically relevant flow scales rather than by the speed of sound, while still retaining the fully compressible formulation. For the spatial discretisation, high-order discontinuous Galerkin (DG) methods provide several advantages for atmospheric modelling, including high-order accuracy, local conservation properties, geometric flexibility, and suitability for non-conforming adaptive meshes [6], [8], [9]. Moreover, DG discretisations are particularly well suited to modern high-performance computing architectures because of their element-local structure [10] and efficient matrix-free implementation strategies [11].

A high-order IMEX-DG solver for compressible atmospheric flows has recently been developed within the deal.II finite element library [12]–[14]. The framework combines IMEX time integration with matrix-free DG operator evaluation, adaptive mesh refinement, and distributed-memory parallelism, enabling scalable simulations of compressible atmospheric dynamics on non-conforming meshes [15], [16]. However, the current formulation is restricted to non-rotating Cartesian configurations.

For atmospheric applications and numerical weather prediction (NWP), planetary rotation is a leading-order physical effect and a key ingredient of geostrophic balance, i.e. the equilibrium between horizontal pressure-gradient and the Coriolis force that governs much of the large-scale atmospheric circulation [1], [2]. Extending IMEX-DG atmospheric solvers to rotating flows is therefore essential both for standard benchmark configurations involving rotating stratified flows and to proceed towards more realistic atmospheric applications.

In this work, we extend the existing IMEX-DG framework to include the Coriolis force within the fully compressible Euler equations. Two alternative treatments of the rotational term within the scheme’s solution strategy, based on a nonlinear fixed-point formulation, are introduced and compared in terms of accuracy, robustness, and computational efficiency.

The proposed rotating formulation is validated through two- and three-dimensional benchmark problems. In particular, inertia-gravity wave tests are used to assess convergence properties and preservation of balanced dynamics, while fully three-dimensional flow over orography is considered to investigate performance in a rotating stratified atmosphere using established benchmarks and newly design large-scale tests. Simulations are performed on both uniform and, to the authors’ knowledge for the first time, non-conforming adaptive meshes. The results show that the inclusion of rotation does not adversely affect the accuracy, robustness, and scalability properties of the original IMEX-DG framework while correctly capturing the expected rotating-flow dynamics.

The remainder of the paper is organized as follows. The governing equations are presented in Section 2. The temporal and spatial discretisation methods are described in Section 3. Numerical results in two and three dimensions, including flow over orography and the comparison between the two rotational treatments, are presented in Section 4. Some conclusions and considerations about future work are presented in Section 5.

2 The mathematical model↩︎

We consider a dry atmosphere modeled as an inviscid fluid governed by the compressible Euler equations on a bounded domain \(\Omega \subset \mathbb{R}^d\), with \(2 \le d \le 3\) [2], [17]. The spatial coordinate is denoted by \(\boldsymbol{x}\in\Omega\) and time by \(t \in [0,T_f]\). The system is supplemented by suitable initial and boundary conditions, which will be specified according to the test cases considered.

The conservation laws of mass, momentum, and total energy read \[\begin{align} \partial_t \rho + \nabla\cdot(\rho\boldsymbol{u}) &= 0, \tag{1} \\ \partial_t(\rho\boldsymbol{u}) + \nabla\cdot(\rho\boldsymbol{u}\otimes\boldsymbol{u}) + \nabla p &= \rho\boldsymbol{a}_{\text{ext}}, \tag{2} \\ \partial_t(\rho E) + \nabla\cdot\!\big((\rho E+p)\boldsymbol{u}\big) &= \rho\boldsymbol{a}_{\text{ext}}\cdot\boldsymbol{u}, \tag{3} \end{align}\] where \(\rho\) is the density, \(\boldsymbol{u}=(u,v,w)\) is the velocity field, \(p\) is the pressure, and \(E\) is the total specific energy defined as \(E = e + \frac{1}{2}|\boldsymbol{u}|^2,\) with \(e\) the internal energy per unit mass. Furthermore, \(\partial_t\) denotes the time derivative, while \(\nabla\cdot\big(\cdot\big)\) and \(\nabla\big(\cdot\big)\) denote the three-dimensional divergence and gradient operators, respectively. The system is closed by assuming a calorically perfect ideal gas, for which \[p = \rho R T = (\gamma - 1)\rho e,\] with \(\gamma = c_{p}/c_{v}\) the constant ratio of specific heats at constant pressure and constant volume and \(R = c_{p} - c_{v}\) the specific gas constant [18].

We adopt a local Cartesian tangent-plane coordinate system \((x,y,z)\) aligned with the eastward, northward, and upward directions. When expressing the equations in a rotating non-inertial reference frame, additional apparent forces arise [2]. Under the \(f\)-plane approximation, the Coriolis parameter is constant, \[f = f_0.\]

The centrifugal contribution is absorbed into the definition of the gravitational acceleration, so that the external body force becomes \[\boldsymbol{a}_{\text{ext}} = - g \boldsymbol{k} - f\,\boldsymbol{k}\times\boldsymbol{u},\] where \(\boldsymbol{k}\) denotes the unit vector in the vertical direction and \(g\) is the gravitational acceleration. The term \(-f\,\boldsymbol{k}\times\boldsymbol{u}\) represents the Coriolis acceleration. Since the Coriolis force is orthogonal to the velocity field, i.e. \((\boldsymbol{k} \times \boldsymbol{u})\cdot \boldsymbol{u}= 0,\) it does not contribute to the work of external forces in the energy equation.

Using the identity \((\rho E + p)\boldsymbol{u}= (h + k)\rho \boldsymbol{u}\), where \(k = |\boldsymbol{u}|^2/2\) and \(h\) is the specific enthalpy defined as \(h = e + p/\rho\), following the formulation adopted in [12], [19], the system can be equivalently rewritten as \[\label{eq:governing95equations952} \begin{align} \partial_t \rho \;+\; \nabla\!\cdot(\rho \boldsymbol{u}) &= 0, \\[4pt] \partial_t(\rho \boldsymbol{u}) \;+\; \nabla\!\cdot(\rho \boldsymbol{u}\otimes\boldsymbol{u}) \;+\; \nabla p &= -\,\rho g\,\boldsymbol{k} \;-\; \rho f\,\boldsymbol{k}\times\boldsymbol{u}, \\[4pt] \partial_t(\rho E) \;+\; \nabla\!\cdot\!\big[(h + k)\,\rho \boldsymbol{u}\big] &= -\,\rho g\,\boldsymbol{k}\!\cdot\!\boldsymbol{u}. \end{align}\tag{4}\]

For two-dimensional configurations, we employ a 2D vertical slice approximation that restricts the dynamics to the \(x-z\) plane, while still allowing for a nonzero meridional velocity component \(v\). In practice, this assumption is implemented by imposing \(\partial _y (\cdot)=0\) in 4 , while the transverse velocity \(v\) remains an active prognostic variable. This approach is often used in atmospheric dynamics to reduce the full three-dimensional problem to a two-dimensional slice that still captures the essential balance between advection, pressure gradients, buoyancy, and rotation [5], [7]. Using this approximation, the equations, written in scalar form, become as follows

\[\begin{align} \partial_t \rho + \partial_x (\rho u) + \partial_z (\rho w)& = 0,\\ \partial_t (\rho u) + \partial_x (\rho u^2) + \partial_z (\rho u w) + \partial_x p &=f \rho v,\\ \partial_t (\rho v) + \partial_x (\rho u v) + \partial_z (\rho v w) &= - f \rho u,\\ \partial_t (\rho w) + \partial_x (\rho u w) + \partial_z (\rho w^2) + \partial_z p &= -g\rho,\\ \partial_t (\rho E) + \partial_x \big[(h+k)\rho u \big] + \partial_z \big[(h+k)\rho w \big] &= - g \rho w. \label{eq:vertical95slice} \end{align}\tag{5}\]

3 The numerical framework↩︎

The numerical approximation of the governing equations 4 is based on a high-order Discontinuous Galerkin (DG) spatial discretization coupled with an implicit-explicit (IMEX) time integration strategy. This combination is well suited for atmospheric flows, where slow advective dynamics coexist with fast acoustic waves. In the low Mach number regime, terms proportional to \(1/\text{Ma}^2\) in the non-dimensional formulation of the equations, such as the pressure gradient, introduce stiffness, leading to severe stability restrictions for fully explicit schemes [9], [20], [21].

Following [12], [19], time discretization is performed using an IMEX additive Runge–Kutta (ARK) method, in which stiff contributions are treated implicitly and non-stiff terms explicitly [22]–[24]. In particular, pressure-gradient terms, together with gravity and rotation, are included in the implicit part, while advective and nonlinear transport terms are handled explicitly. This choice is motivated by the low-Mach nature of atmospheric flows, where fast acoustic modes only carry a small fraction of the total energy: treating them implicitly relaxes stability constraints while retaining accuracy and efficiency [12], [25]. Although the Coriolis term and the gravity term are not intrinsically stiff, they are commonly included in the implicit treatment to better preserve the hydrostatic and geostrophic balances.

In this work, we employ a second-order IMEX-ARK scheme, defined by two sets of coefficients, associated with the explicit and implicit parts, respectively. More specifically, we adopt the TR-BDF2 scheme [26]. The corresponding Butcher tableaux can be written as

with parameter \(\chi=2-\sqrt{2}\), which ensures \(L\)-stability and thus provides effective damping of fast acoustic modes in stiff low-Mach regimes [26]–[28].

For the explicit part, we employ a second-order RK scheme [6]. The conditions \(b=\tilde{b}\) and \(c=\tilde{c}\) are enforced, where \(b,\widetilde{b}\) are the final update weights and \(c,\widetilde{c}\) the stage nodes of the explicit and implicit Runge-Kutta tableaux, respectively [23]. The condition \(b=\tilde{b}\) guarantees preservation of linear invariants [6], while \(c=\tilde{c}\) ensures consistency between the explicit and implicit stages. The remaining free parameter in the explicit tableau is set to \(a_{32}=0.5\), following [19], so as to improve monotonicity properties. The resulting scheme belongs to the Explicit Singly Diagonal Implicit Runge-Kutta (ESDIRK) class. Thus, for a ODE of the form \(\mathbf{y}'=\mathbf{f}_{\mathrm{S}}(\mathbf{y},t)+\mathbf{f}_{\mathrm{NS}}(\mathbf{y},t)\), where NS denotes the non-stiff term and S denotes the stiff one, the \(l\)-th stage of the solution is computed as \[\begin{align} \mathbf{y}^{(n,l)}=\mathbf{y}^{n}+&\Delta t \sum_{m=1}^{l-1}[a_{lm}\mathbf{f}_{\mathrm{NS}}(\mathbf{y}^{(n,m)},t^n+c_m\Delta t)+\widetilde{a}_{lm}\mathbf{f}_{\mathrm{S}}(\mathbf{y}^{(n,m)},t^n+c_m\Delta t)] \\ +& \Delta t \, \widetilde{a}_{ll} \, \mathbf{f}_{\mathrm{S}}(\mathbf{y}^{(n,l)},t^n+c_l\Delta t), \end{align}\] for \(l=1,2,3\). The update is then given by \[\mathbf{y}^{n+1} = \mathbf{y}^n + \Delta t \sum_{l=1}^{3} b_l \left[ \mathbf{f}_{\mathrm{NS}}(\mathbf{y}^{(n,l)}, t^n + c_l \Delta t) + \mathbf{f}_{\mathrm{S}}(\mathbf{y}^{(n,l)}, t^n + c_l \Delta t) \right].\]

Applying the IMEX-ARK method presented above, the first stage is purely formal and simply reduces to \(\rho^{(n,1)}=\rho^n\) for the density and analogously for the other prognostic variables. From the second stage onward, it is possible to write the following general formulae \[\begin{align} \rho^{(n,l)} &= \rho^n - \sum_{m=1}^{l-1}a_{lm}\Delta t \, \nabla \cdot (\rho^{(n,m)} \boldsymbol{u}^{(n,m)}), \\ \rho^{(n,l)} \boldsymbol{u}^{(n,l)} & + \tilde{a}_{ll}\Delta t \nabla p^{(n,l)} +g \,\tilde{a}_{ll}\Delta t \rho^{(n,l)} \boldsymbol{k}+ f\tilde{a}_{ll} \Delta t \rho^{(n,l)}\boldsymbol{k}\times \boldsymbol{u}^{(n,l)} = \mathbf{m}^{(n,l)}, \\ \rho^{(n,l)} E^{(n,l)} &+ \tilde{a}_{ll}\Delta t \, \nabla \cdot (h^{(n,l)} \rho^{(n,l)} \boldsymbol{u}^{(n,l)}) + g\,\tilde{a}_{ll}\rho^{(n,l)} \boldsymbol{k} \cdot \boldsymbol{u}^{(n,l)}=\hat{e}^{(n,l)}, \\ \end{align} \label{eq:IMEXsystem95sstage}\tag{6}\]

where \[\begin{align} \mathbf{m}^{(n,l)} = \rho^n \boldsymbol{u}^n - \sum_{m=1}^{l-1} &[a_{lm}\Delta t \, \nabla \cdot (\rho^{(n,m)} \boldsymbol{u}^{(n,m)} \otimes \boldsymbol{u}^{(n,m)})+ \tilde{a}_{lm}\Delta t \nabla p^{(n,m)} \\ &\quad + g \tilde{a}_{lm}\Delta t \rho^{(n,m)} \boldsymbol{k} +f\tilde{a}_{lm}\Delta t \rho^{(n,m)} \boldsymbol{k}\times \boldsymbol{u}^{(n,m)}], \end{align}\] \[\begin{align} \hat{e}^{(n,l)} = &\rho^n E^n - \sum_{m=1}^{l-1} \Bigg[\tilde{a}_{lm}\Delta t \, \nabla \cdot (h^{(n,m)} \rho^{(n,m)} \boldsymbol{u}^{(n,m)})\\& \quad + a_{lm}\Delta t \,\nabla \cdot (\kappa^{(n,m)} \rho^{(n,m)} \boldsymbol{u}^{(n,m)}) +g\Delta t \,\tilde{a}_{lm}\rho^{(n,m)} \boldsymbol{k} \cdot \boldsymbol{u}^{(n,m)}\Bigg]. \end{align}\]
The spatial discretization of the governing equations is based on a high-order Discontinuous Galerkin (DG) method on a mesh \(\mathcal{T}_h\), composed of quadrilateral elements in 2D and hexaedral elements in 3D. This choice combines high-order accuracy with local conservation and geometric flexibility, making it well suited for compressible atmospheric flows [29].

For an integer \(r \geq 0\), we define the polynomial space \(\mathbb{Q}_r(K)\) as the space of polynomials of degree

\(\leq r\) in each coordinate direction on the element \(K\). The global DG finite element space is then \[Q_r = \{ v \in L^2(\Omega) \;: \;v|_K \in \mathbb{Q}_r(K), \;\forall K \in \mathcal{T}_h \}.\] Scalar quantities such as density and pressure are approximated in \(Q_r\), while vector-valued quantities such as velocity in \([Q_r]^d= \{ v = (v_1,\dots,v_d): v_i \in Q_r \},\) where \(d\) denotes the number of the physical dimensions. We employ a nodal DG formulation on tensor-product elements. On each element \(K\), the discrete solution is \[\mathbf{q}_h\big|_K(\boldsymbol{x},t) = \sum_{j=1}^{N_p} \mathbf{q}_j(t) \ell_j(\boldsymbol{x}),\] where \(\{\mathbf{q}_j \approx \mathbf{q}(\boldsymbol{x}_j)\}_{j=1}^{N_p}\), with \(N_p=(r+1)^d\), are nodal values at Gauss-Legendre points \(\{\boldsymbol{x}_j\}\), and \(\{\ell_j\}\) are the associated Lagrange basis functions [30]. This tensor-product nodal formulation is particularly suitable for matrix-free implementations, since the discrete operators can be evaluated directly at quadrature points without explicit assembly of the global matrices. In addition, the tensor-product structure enables efficient sum-factorized quadrature and operator evaluation on quadrilateral and hexahedral meshes [11], [29].

Following the standard DG formulation, the equations are multiplied by test functions and integrated by parts element-wise, leading to surface flux contributions that couple neighbouring elements. A Rusanov-type flux is used at interfaces, with stabilization parameter \(\lambda = |\mathbf{u}\cdot\mathbf{n}|\) [19]. This avoids excessive dissipation proportional to the sound speed. For the pressure-gradient contribution, centered fluxes are employed in order to simplify the implicit treatment within the IMEX formulation, as implicit time integration inherently ensures stability without the need for upwind stabilization.

Denoting \(\{\boldsymbol{\varphi}_i\}_{i=1}^{\dim([Q_r]^d)}\) the basis functions of the DG velocity space \([Q_r]^d\) and \(\{\psi_j\}_{j=1}^{\dim(Q_r)}\) the basis functions of the DG pressure space \({Q}_r\), velocity and pressure can be rewritten in a discrete form as \[\begin{align} \boldsymbol{u}(\boldsymbol{x},t)\approx\sum_{j=1}^{\text{dim}([{Q}_r]^d)} u_j(t) \boldsymbol{\varphi}_j(\boldsymbol{x}), & \qquad p(\boldsymbol{x},t)\approx \sum_{j=1}^{\text{dim}(Q_r)} p_j(t) {\psi}_j(\boldsymbol{x}). \end{align}\]

Given these definitions, it is possible to derive the weak formulation of the momentum equation at stage \(l\) in matrix form:

\[\mathbf{A}^{(n,l)}\,\mathbf{u}^{(n,l)} + \mathbf{B}^{(n,l)} \mathbf{p}^{(n,l)} + \mathbf{R}^{(n,l)}\mathbf{u}^{(n,l)}= \mathbf{f}^{(n,l)}- \mathbf{f}_g^{(n,l)}, \label{eq:momentum95matrix}\tag{7}\] where \(\mathbf{u}^{(n,l)}\) and \(\mathbf{p}^{(n,l)}\) are the vectors of velocity and pressure degrees of freedom, \[A^{(n,l)}_{ij} = \sum_K \int_K \rho^{(n,l)} \boldsymbol{\varphi}_j \cdot \boldsymbol{\varphi}_i \, \mathrm{d} \Omega, \label{eq:A}\tag{8}\]

\[B^{(n,l)}_{ij} = \sum_K \int_K -\tilde{a}_{ll}\Delta t \, \nabla \cdot \boldsymbol{\varphi}_i \, \psi_j \, \mathrm{d} \Omega + \sum_{\Gamma \in E} \int_\Gamma \tilde{a}_{ll}\Delta t \, \{\!\!\{{ \psi_j }\}\!\!\} \, \llbracket{\boldsymbol{\varphi}_i}\rrbracket \, \mathrm{d} \Sigma,\]

\[R^{(n,l)}_{ij}= \tilde{a}_{ll}\Delta t f \sum_{K\in\mathcal{T}_h}\!\int_K \rho^{(n,l)} ( \boldsymbol{k} \times \boldsymbol{\varphi}_j ) \cdot \boldsymbol{\varphi}_i\, \mathrm{d} \Omega\]

\[\mathbf{f}_{g_{i}}^{(n,l)}= \tilde{a}_{ll}\Delta t g \sum_{K\in\mathcal{T}_h}\!\int_K \rho^{(n,l)}\boldsymbol{k} \cdot \boldsymbol{\varphi}_{i}\, \mathrm{d} \Omega,\]

\[\begin{align} \mathbf{f}^{(n,l)} =& \sum_{K \in \mathcal{T}_h} \int_K \rho^n \boldsymbol{u}^n \cdot \boldsymbol{\varphi}_i \, \mathrm{d} \Omega + \sum_{m=1}^{l-1}\Bigg[ \sum_{K \in \mathcal{T}_h} \int_K a_{lm}\Delta t \, \left(\rho^{(n,m)} \boldsymbol{u}^{(n,m)} \otimes \boldsymbol{u}^{(n,m)}\right) : \nabla \boldsymbol{\varphi}_i \, \mathrm{d}\Omega \\ +& \sum_{K \in \mathcal{T}_h} \int_K \tilde{a}_{lm}\Delta t p^{(n,m)} \nabla\!\cdot \boldsymbol{\varphi}_i \mathrm{d} \Omega - \sum_{\Gamma \in \mathcal{E}} \int_\Gamma a_{lm}\Delta t \, \{\!\!\{{ \rho^{(n,m)} \boldsymbol{u}^{(n,m)} \otimes \boldsymbol{u}^{(n,m)}}\}\!\!\} : \langle\!\langle {\boldsymbol{\varphi}_i} \rangle\!\rangle \, \mathrm{d} \Sigma \\ -& \sum_{\Gamma \in \mathcal{E}} \int_\Gamma \tilde{a}_{lm}\Delta t \, \{\!\!\{{ p^{(n,m)}}\}\!\!\} \, \llbracket{\boldsymbol{\varphi}_i}\rrbracket \, \mathrm{d} \Sigma -\sum_{\Gamma \in \mathcal{E}} \int_\Gamma a_{lm}\Delta t \, \frac{\lambda^{(n,m)}}{2} \, \llbracket{ \rho^{(n,m)} \boldsymbol{u}^{(n,m)} }\rrbracket : \langle\!\langle {\boldsymbol{\varphi}_i} \rangle\!\rangle \, \mathrm{d} \Sigma \\ -& \sum_{K \in \mathcal{T}_h} \int_K g \tilde{a}_{lm}\Delta t \, \rho^{(n,m)} \boldsymbol{k} \cdot \boldsymbol{\varphi}_i \, \mathrm{d} \Omega - \sum_{K \in \mathcal{T}_h} \int_K f \tilde{a}_{lm}\Delta t \rho^{(n,m)}(\boldsymbol{k}\times \boldsymbol{u}^{(n,m)})\cdot\boldsymbol{\varphi}_i \mathrm{d} \Omega \Bigg]. \end{align}\] with \(\lambda^{(n,m)} = \max \left( | \boldsymbol{u}^{(n,m),+}\cdot \mathbf{n}^+|, \; | \boldsymbol{u}^{(n,m),-}\cdot \mathbf{n}^- | \right),\) and the symbol \(:\) indicates the scalar product for second-order tensors. Taking into account that \[\rho^{(n,l)} E^{(n,l)} = \rho^{(n,l)} e^{(n,l)}(p^{(n,l)},\rho^{(n,l)}) + \rho^{(n,l)} \kappa^{(n,l)},\] it is also possible to derive the algebraic form of the energy equation as \[\mathbf{C}^{(n,l)} \mathbf{u}^{(n,l)} + \mathbf{D}^{(n,l)} \mathbf{p}^{(n,l)}+\mathbf{k}^{(n,l)}+\mathbf{M}^{(n,l)}_g \mathbf{u}^{(n,l)}= \mathbf{g}^{(n,l)},\] where \[C^{(n,l)}_{ij} = \sum_{K\in \mathcal{T}_h} \int_K -\tilde{a}_{ll}\Delta t \, h^{(n,l)} \rho^{(n,l)} \boldsymbol{\varphi}_j \cdot \nabla \psi_i \, \mathrm{d} \Omega + \sum_{\Gamma \in E} \int_\Gamma \tilde{a}_{ll}\Delta t \, \{\!\! \{ h^{(n,l)} \rho^{(n,l)} \boldsymbol{\varphi}_j \}\!\! \} \cdot \llbracket{\psi_i}\rrbracket \, \mathrm{d} \Sigma,\] \[D^{(n,l)}_{ij}=\sum_{K\in \mathcal{T}_h}\int_K \frac{1}{(\gamma -1)}(\psi_i)\psi_j \mathrm{d}\Omega ,\]

\[{k}^{(n,l)}_i=\sum_{K \in \mathcal{T}_h}\int_K \rho^{(n,l)} k^{(n,l)}\mathbf{\psi}_i\mathrm{d}\Omega,\]

\[[\mathbf{M}_{g}^{(n,l)}]_{ij}=\sum_{K\in\mathcal{T}_h}\int_K g\tilde{a}_{ll}\Delta t\,\rho^{(n,l)} \boldsymbol{k}\cdot\varphi_i\,\psi_j\,\mathrm{d}\Omega\] and the expression of \(\mathbf{g}\) is given by

\[\begin{align} \mathbf{g}_i^{(n,l)} &= \sum_{K\in\mathcal{T}_h}\int_K \rho^n E^n \psi_i \,\mathrm{d}\Omega +\sum_{m=1}^{l-1}\Bigg[ \sum_{K\in\mathcal{T}_h}\int_K \Delta t a_{lm}\,\left(\kappa^{(n,m)} \rho^{(n,m)}\boldsymbol{u}^{(n,m)}\right)\cdot\nabla\psi_i\,\mathrm{d}\Omega\\ & \quad +\sum_{K\in\mathcal{T}_h}\int_K \Delta t\tilde{a}_{lm}\left(h^{(n,m)}\rho^{(n,m)}\boldsymbol{u}^{(n,m)}\right) \cdot\nabla\psi_i\,\mathrm{d}\Omega \\ &\quad -\sum_{\Gamma\in\mathcal{E}}\int_\Gamma \Delta t\Big( a_{lm}\,\{\!\!\{{\kappa^{(n,m)}\rho^{(n,m)}\boldsymbol{u}^{(n,m)}}\}\!\!\} +\tilde{a}_{lm}\{\!\!\{{h^{(n,m)}\rho^{(n,m)}\boldsymbol{u}^{(n,m)}}\}\!\!\} \Big)\cdot\llbracket{\psi_i}\rrbracket\,\mathrm{d}\Sigma \\ &\quad -\sum_{\Gamma\in\mathcal{E}}\int_\Gamma \Delta t\,\frac{\lambda^{(n,m)}}{2} \Big(a_{lm}\llbracket{\rho^{(n,m)}\kappa^{(n,m)}}\rrbracket +\tilde{a}_{lm}\llbracket{\rho^{(n,m)}e^{(n,m)}}\rrbracket\Big)\cdot\llbracket{\psi_i}\rrbracket\,\mathrm{d}\Sigma \\ &\quad -\sum_{K\in\mathcal{T}_h}\int_K g \,\tilde{a}_{lm}\Delta t\,\rho^{(n,m)}\boldsymbol{k}\cdot\boldsymbol{u}^{(n,m)}\psi_i \,\mathrm{d}\Omega - \frac{\lambda^{(n,m)}\tilde{a}_{lm}\Delta t}{2}\sum_{\Gamma \in \mathcal{E}}\int_\Gamma \llbracket{\rho^{(n,m)}E^{(n,m)}}\rrbracket\cdot\llbracket{\psi_i}\rrbracket\mathrm{d}\Sigma \Bigg]\\ &\quad -\frac{\lambda^{(n,l)}\tilde{a}_{ll}\Delta t}{2}\sum_{\Gamma \in \mathcal{E}}\int_\Gamma\llbracket{\rho^{(n,l)}e^{(n,l)}\psi_i}\rrbracket\cdot\llbracket{\psi_j}\rrbracket\mathrm{d}\Sigma. \end{align}\] with \[\lambda^{(n,m)} = \max \left( | \boldsymbol{u}^{(n,m),+}\cdot \mathbf{n}^+|, \; | \boldsymbol{u}^{(n,m),-}\cdot \mathbf{n}^- | \right).\] Then, we end up with the following coupled system of equations \[\begin{align} &\mathbf{A}^{(n,l)}\,\mathbf{u}^{(n,l)}+\mathbf{R}^{(n,l)}\mathbf{u}^{(n,l)} + \mathbf{B}^{(n,l)} \mathbf{p}^{(n,l)} + \mathbf{f}_g^{(n,l)}= \mathbf{f}^{(n,l)},\\ &\mathbf{C}^{(n,l)} \mathbf{u}^{(n,l)} +\mathbf{M}^{(n,l)}_g \mathbf{u}^{(n,l)} + \mathbf{D}^{(n,l)} \mathbf{p}^{(n,l)}+\mathbf{k}^{(n,l)}= \mathbf{g}^{(n,l)}. \label{eq:matrix95form} \end{align}\tag{9}\]

The algebraic system for \((\mathbf{u}^{(n,l)},\mathbf{p}^{(n,l)})\) is nonlinear and coupled, as matrices and right-hand sides depend on the stage unknowns through thermodynamic quantities. Following the strategy commonly adopted in semi-implicit atmospheric solvers [12], [19], [31], the nonlinear system is solved through a fixed-point (Picard) iteration[30]. This approach avoids the solution of a fully nonlinear system at each IMEX stage while preserving the implicit coupling between pressure and velocity. At each iteration, selected terms are evaluated using the previous iterate, resulting in a sequence of linear systems. The iteration is initialized using the solution obtained at the previous IMEX stage and proceeds until the relative variation between two consecutive iterates falls below a prescribed tolerance. Different lagging strategies for the nonlinear terms result in separate formulations:

• In the first approach (R1), the term \(\mathbf{R}^{(n,l)}\mathbf{u}^{(n,l)}\) in the momentum equation is evaluated at the previous Picard iterate \(k\), while the remaining terms are treated at the current iteration \(k+1\). The same strategy is adopted for the gravitational contribution in the energy equation, where \(\mathbf{M}_g^{(n,l)}\mathbf{u}^{(n,l)}\) and \(\mathbf{k}^{(n,l)}\) are lagged at iteration \(k\). This choice moves the rotational and gravitational contributions to the right-hand side, where they are evaluated using the previous Picard iterate. The resulting system reads \[\begin{align} &\mathbf{A}^{(n,l)}\,\mathbf{u}^{(n,l,k+1)}+ \mathbf{B}^{(n,l)} \mathbf{p}^{(n,l,k+1)} = \mathbf{f}^{(n,l)}- \mathbf{f}_g^{(n,l)}- \mathbf{R}^{(n,l)}\mathbf{u}^{(n,l,k)} ,\\ &\mathbf{C}^{(n,l,k)}\mathbf{u}^{(n,l,k+1)}+\mathbf{D}^{(n,l)}\mathbf{p}^{(n,l,k+1)}=\mathbf{g}^{(n,l,k)}-\mathbf{k}^{(n,l,k)}-\mathbf{M}_g^{(n,l)}\mathbf{u}^{(n,l,k)}. \label{sbhpfndc} \end{align}\tag{10}\]

  This coupled system is solved using a Schur complement approach. Eliminating the velocity from the momentum equation yields a reduced system for the pressure, \[\left[\mathbf{D}^{(n,l)} - \mathbf{C}^{(n,l,k)} \big(\mathbf{A}^{(n,l)}\big)^{-1} \mathbf{B}^{(n,l)}\right]\, \mathbf{p}^{(n,l,k+1)} = \mathbf{RHS}^{(n,l,k)}_{R1},\] where \[\begin{align} \mathbf{RHS}^{(n,l,k)}_{R1} =&\;\mathbf{g}^{(n,l,k)} - \mathbf{M}_g^{(n,l)}\mathbf{u}^{(n,l,k)} - \mathbf{k}^{(n,l,k)} \\ &\; - \mathbf{C}^{(n,l)} \big(\mathbf{A}^{(n,l)}\big)^{-1} \Big(\mathbf{f}^{(n,l)} - \mathbf{R}^{(n,l)}\mathbf{u}^{(n,l,k)} - \mathbf{f}_g^{(n,l)}\Big). \end{align}\]

• In the second strategy (R2), both the rotational term \(\mathbf{R}^{(n,l)}\mathbf{u}^{(n,l)}\) and the gravitational contribution \(\mathbf{M}_g^{(n,l)}\mathbf{u}^{(n,l)}\) are evaluated at the current Picard iteration \(k+1\) and retained on the left-hand side. These contributions are therefore treated implicitly within the coupled momentum-energy system.

  Proceeding as before, the resulting linear system for the pressure reads \[\big(\mathbf{D}^{(n,l)} - (\mathbf{C}^{(n,l,k)}+\mathbf{M}_g^{(n,l)}) \big(\mathbf{A}^{(n,l)}+\mathbf{R}^{(n,l)}\big)^{-1} \mathbf{B}^{(n,l)}\big)\, \mathbf{p}^{(n,l,k+1)} = \mathbf{RHS}^{(n,l,k)}_{R2},\] with \[\begin{align} \mathbf{RHS}^{(n,l,k)}_{R2} =&\;\mathbf{g}^{(n,l,k)} - \mathbf{k}^{(n,l,k)} \\ &\; - (\mathbf{C}^{(n,l,k)}+\mathbf{M}_g^{(n,l)}) \big(\mathbf{A}^{(n,l)}+\mathbf{R}^{(n,l)}\big)^{-1} \Big(\mathbf{f}^{(n,l)} - \mathbf{f}_g^{(n,l)}\Big). \end{align}\] Once \(\mathbf{p}^{(n,l,k+1)}\) has been computed, the velocity \(\mathbf{u}^{(n,l,k+1)}\) is recovered from the momentum equation.

In both approaches, the pressure system is solved using a GMRES iterative method, with a stopping criterion based on the relative residual not exceeding a prescribed tolerance.

3.1 Matrix-free application of \((\mathbf{A}+\mathbf{R})^{-1}\)↩︎

The formulations discussed above require repeated applications of block-diagonal operators such as \(\mathbf{A}^{-1}\) and \((\mathbf{A}+\mathbf{R})^{-1}\). These operators are evaluated in a matrix-free fashion by exploiting the tensor-product structure of the DG discretisation on quadrilateral/hexahedral meshes [11], [32].

The density-weighted mass matrix can be expressed element-wise as \[\mathbf{A}_K = \mathbf{S}^T \mathbf{W}_\rho \mathbf{S},\] where \(\mathbf{S}\) is the interpolation matrix between nodal and quadrature bases, and \(\mathbf{W}_\rho\) is a diagonal matrix containing quadrature weights, geometric factors, and density [11]. This factorisation enables an efficient application of \(\mathbf{A}_K^{-1}\) through tensor-product operations.

For the R2 strategy, the inverse of the modified operator \((\mathbf{A}+\mathbf{R})\) is required. Rather than treating \(\mathbf{R}\) as a generic matrix contribution, we exploit the algebraic structure of the Coriolis operator. Introducing the linear mapping \(\mathcal{J}\mathbf{v} = \mathbf{k}\times\mathbf{v}\), one can show that \[\mathbf{R} = \beta\, \mathbf{A}\mathbf{J}, \qquad \beta = \tilde{a}_{ll}\Delta t f,\] so that \[\mathbf{A}+\mathbf{R} = \mathbf{A}(\mathbf{I}+\beta\mathbf{J}).\] This immediately yields the factorisation \[(\mathbf{A}+\mathbf{R})^{-1} = (\mathbf{I}+\beta\mathbf{J})^{-1}\mathbf{A}^{-1}.\]

Assuming a constant vertical rotation axis \(\mathbf{k}=(0,0,1)^T\), the operator \(\mathcal{J}\) acts as a planar rotation on the horizontal velocity components, satisfying \(\mathcal{J}^2=-I\). Hence, \[(\mathbf{I}+\beta\mathbf{J})^{-1} = \frac{1}{1+\beta^2} \begin{pmatrix} 1 & \beta & 0 \\ -\beta & 1 & 0 \\ 0 & 0 & 1+\beta^2 \end{pmatrix},\] which corresponds to a local transformation applied pointwise to the velocity field. Therefore, the application of \((\mathbf{A}+\mathbf{R})^{-1}\) reduces to a standard matrix-free application of \(\mathbf{A}^{-1}\), followed by a local algebraic transformation involving \((\mathbf{I}+\beta\mathbf{J})^{-1}\). This avoids the need to solve linear systems with the full coupled operator and preserves the efficiency of the matrix-free implementation.

4 Numerical results↩︎

The numerical method presented in the previous sections is validated through a series of two- and three-dimensional test cases for rotating atmospheric flows. The considered benchmarks are designed to evaluate the ability of the proposed IMEX-DG formulation to accurately capture wave propagation and balanced rotating dynamics while preserving long-time stability. The numerical experiments include both standard rotating inertia-gravity wave benchmarks and fully three-dimensional simulations of stratified flow over orography. The latter provide a more challenging setting to assess the behaviour of the method in realistic large-scale atmospheric regimes. Simulations are performed on both uniform and non-conforming meshes in order to evaluate the robustness of the formulation in the presence of local mesh refinement.

Particular attention is devoted to verifying that the inclusion of rotational effects preserves the accuracy, robustness, and scalability properties of the original non-rotating IMEX-DG solver implemented within the matrix-free distributed-memory framework of the deal.II finite element library [33].

For the numerical experiments, we define the minimum mesh element diameter \(\mathcal{H}\) as \[\mathcal{H}=\min{\{\text{diam}(K)|K \in\mathcal{T}_h}\},\] together with the acoustic and advective Courant numbers, defined respectively as \[C=\frac{rc \Delta t}{\mathcal{H}}\sqrt{d}, \qquad C_{\mathrm{adv}}=\frac{r U \Delta t}{\mathcal{H}}\sqrt{d},\] where \(c\) is the reference magnitude of speed of sound, \(U\) is the reference magnitude of the flow velocity, \(r\) is the polynomial degree of the DG method and \(\mathcal{H}\) is interpreted as an effective spatial resolution for a DG discretisation [12]. In all the simulations, we set \(\gamma=1.4\), \(g=9.81\text{m s}^{-2}\), and \(R=287 \text{ J Kg}^{-1}\text{K}^{-1}\).

Simulations have been performed on computing clusters at the Technical University of Denmark on 28 Lenovo ThinkSystem SD530 compute nodes equipped with two Intel Xeon Gold 6226R processors (16 physical cores, 2.90 GHz), [34] and on the MeluXina HPCF, using up to 1024 2x AMD EPYC Rome 7H12 64c 2.6GHz CPUs [35], compiled with GCC 12.3, targeting 256-bit vectorization.

4.1 Two-dimensional inertia-gravity wave↩︎

For the two-dimensional tests, a vertical-slice approximation is adopted to retain the effect of rotation on the vertical structure of the flow while reducing the problem to a two-dimensional physical domain. The main benchmark considered is the propagation of inertia-gravity waves in a rectangular channel with no background flow and non-zero rotation, following the standard test cases proposed in [7], [36]. The computational domain of the first case is a two-dimensional channel of size \(\SI[parse-numbers=false]{(0,6000)\times(0,10)}{\kilo\meter\squared}\). Periodic boundary conditions are prescribed in the horizontal direction, while wall boundary conditions are imposed at the top and bottom boundaries. The flow is initialized with a stationary, hydrostatic, horizontally homogeneous, isothermal background state: \[\begin{align} T_0, \quad p_0(z)=p_s e^{-\delta z}, \quad \rho_0=\rho_s e^{-\delta z}, \quad \delta=\frac{g}{RT_0},\quad p_s=\frac{g \rho_s}{\delta}, \quad \boldsymbol{u}_0=\boldsymbol{0}. \label{eq:baldauf95background} \end{align}\tag{11}\] Following [7], [36], the thermodynamic perturbations are initialized through Bretherton-transformed variables \[\begin{align} & T'(\mathbf{r},t=0)=e^{\tfrac{1}{2}\delta z} \cdot T_b(\mathbf{r},t=0),\\ p'(\mathbf{r},t=0)&=0, \qquad \rho'(\mathbf{r},t=0)=e^{-\tfrac{1}{2}\delta z} \cdot \rho_b(\mathbf{r},t=0), \label{eq:baldauf95bubble} \end{align}\tag{12}\] where the Bretherton-transformed temperature and density are given by \[\begin{align} T_b(\mathbf{r},t=0)=\Delta T \cdot e^{- \frac{(x-x_c)^2}{a^2}}\cdot \sin \pi \frac{z}{H}, \qquad \rho_b(\mathbf{r},t=0)=\rho_s \left( - \dfrac{T_b}{T_0} \right). \end{align}\] The perturbation bubble is centered at \(x_c=\SI{3000}{\kilo\meter}\), with \(a=\SI{100}{\kilo\meter}\). The velocity perturbation is initially set to zero. Rotation is included through a constant Coriolis parameter \(f\approx \SI[parse-numbers=false]{1.03126 \cdot 10^{-4}}{\per\second}\), corresponding to a latitude of \(\phi=\ang{45} °\). The final time is set to \(T_f=8\) h. To reproduce the effective resolution used in the DUNE simulations reported in [36], we consider a mesh composed by \((300 \times 20)\) elements and a polynomial degree \(r=4\), yielding an effective resolution of \(\SI{5}{\kilo\meter}\) in \(x\) and \(\SI{125}{\meter}\) in \(z\). The time step is taken equal to \(\SI{0.5}{\second}\), corresponding to an advective Courant number \(C_{\mathrm{adv}}\approx2.28 \cdot 10^{-6}\) and an acoustic Courant number \(C\approx 0.044\). As discussed in [36], a sufficiently small time step is necessary to properly resolve the acoustic–gravity wave dynamics and avoid spurious fast-wave activity contaminating the balanced solution. Following [36], the analytical reference solution is obtained by linearising the fully compressible equations around an isothermal hydrostatic state, applying a Bretherton transformation, and evolving the resulting inertial gravity wave modes. This yields an exact solution of the linearised system used as a benchmark for fully nonlinear simulations. Using a reference implementation to evaluate the analytical solution for this configuration, Figures  1 (a), 1 (b) and  1 (c) compare the numerical solution against the analytical reference fields at the final time. The results are in line with those of [36], and confirm that the present solver is able to accurately capture the propagation and structure of the inertia–gravity waves emanating from the centre of the domain, where a geostrophic mode remains present until the end of the simulation.

Convergence test and comparison between approaches R1 and R2↩︎

A convergence study is performed to assess the accuracy of the numerical scheme and to compare the two approaches (R1 and R2) introduced in Section 3. Both strategies are validated against the analytical reference solution at the final time \(T_f = \SI{8}{\hour}\), following [36].

The error is evaluated for the perturbation variables \(w\), \(T'\), and \(p'\) using the discrete \(l^2\) and \(l^\infty\) norms \[\begin{align} l^2(\psi) = \left(\frac{1}{N_x N_z}\sum_{i=1}^{N_x}\sum_{k=1}^{N_z}|\psi(x_i,z_k)-\psi_{num}(x_i,z_k)|^2\right)^{1/2},\, l^\infty(\psi) = \max_{i,k} |\psi(x_i,z_k)-\psi_{num}(x_i,z_k)|, \end{align}\] and the experimental order of convergence (EOC) is defined as \[EOC_p = \frac{\log\left(l^p(\Delta x_2)/l^p(\Delta x_1)\right)}{\log\left(\Delta x_2/\Delta x_1\right)},\]

The reference analytical solution is evaluated on a fixed grid with \(N_x = 6000\) and \(N_z = 400\), while the numerical solution is interpolated onto this grid for error computation. The initial discretization consists of \(150 \times 10\) elements, which corresponds to a resolution of \(\SI{40}{\kilo\meter}\) in the \(x\) direction and of \(\SI{1}{\kilo\meter}\) in the \(z\) direction, with \(\Delta t = 1\) s and polynomial degree \(r=4\). The standard hyperbolic scaling of space and time intervals is adopted, in which the mesh resolution and the time step are both concurrently halved. The expected convergence rate is \(\min\{2, r+1\}\).

Acceptable values are obtained for the convergence rates both with the R1 and the R2 approach in the three variables (Tables 1 and 2). More in detail, for \(w\) both approaches achieve near-optimal convergence rates consistent with the expected theoretical order. For \(p'\), a reduced but consistent convergence trend is observed. For some other variables, the convergence rate decreases at the finest refinement levels, suggesting error saturation associated with the evaluation of the analytical reference solution rather than with the numerical discretisation itself. This interpretation is supported by additional self-convergence tests, which recover the expected asymptotic behaviour.

Overall, these results indicate that the two rotational treatments are numerically equivalent for the considered test case, both in terms of convergence properties and solution accuracy.

4.2 Planetary-scale inertia-gravity wave↩︎

To further assess the robustness of the proposed formulation in a more demanding regime, the inertia-gravity wave benchmark is repeated on a larger computational domain. The horizontal extent of the domain is increased by one order of magnitude, yielding a domain of size \(\SI[parse-numbers=false]{(0,60000) \times (0,10)}{\kilo\meter\squared}\). Consistently with the reference configuration considered in the previous test, the horizontal semi-amplitude of the initial perturbation is also increased by a factor of 10, corresponding to \(a=\SI{1000}{\kilo\meter}\). All the remaining physical and numerical parameters are kept unchanged, including the background stratification, the Coriolis parameter, the number of elements of the mesh, the polynomial degree, and the boundary conditions. The simulation is run up to a longer final time \(T_f=\SI{96}{\hour}\) using a time step \(\Delta t = \SI{10}{\second}\). The corresponding acoustic and advective Courant numbers are approximately \(C=0.09\) and \(C_{\text{adv}}=2.81\times 10^{-6}\), respectively.

The larger horizontal scale leads to inertia-gravity waves with longer wavelengths and slower propagation speeds, while rotational effects become comparatively more significant in the overall flow evolution. This configuration therefore provides a more stringent test for the preservation of balanced dynamics and for the long-time robustness of the IMEX-DG formulation.

Figure 2 shows the evolution of the temperature perturbation field at different times during the simulation. Initially, the perturbation is localized around the center of the domain and exhibits the expected vertically elongated structure associated with the imposed stratification. As the simulation evolves, the perturbation generates pairs of inertia-gravity waves propagating horizontally away from the central region. At \(t=\SI{48}{\hour}\), the wave packets are clearly visible on both sides of the initial disturbance while the central balanced structure remains well preserved. At the final time \(t=\SI{96}{\hour}\), the waves have propagated farther across the domain without exhibiting noticeable numerical distortion or spurious oscillations.

Overall, the numerical solution remains stable over the entire simulation and preserves the expected symmetry and large-scale wave structure. These results confirm that the rotating IMEX-DG formulation remains robust and accurate also in large-domain configurations characterized by long integration times and stronger rotational influence.

To further assess the ability of the scheme to preserve balanced dynamics over long integration times, we also provide a quantitative measure of the departure from geostrophic equilibrium during the simulation. Since the present formulation is not exactly well-balanced at the discrete level, the pressure-gradient and Coriolis contributions do not cancel to machine precision, and small imbalance errors may progressively develop during the time integration.

The discrete residual of the horizontal momentum balance is computed, following [37], as

\[\label{eq:egeo} \mathcal{E}_{\mathrm{geo}} = \frac{\int_0^{z_{\text{max}}}\int_0^{x_{\text{max}}}|\frac{\partial p}{\partial x}\frac{1}{\rho}-f v| \mathrm{d}x \mathrm{d}z}{z_{\text{max}}x_{\text{max}}},\tag{13}\] where \(x_{\text{max}}\) and \(z_{\text{max}}\) denote the the maximum horizontal and vertical coordinates of the domain, respectively.

Although a slow growth of \(\mathcal{E}_{\mathrm{geo}}\) is observed over the four-day simulation time, the magnitude of the error remains lower than \(10^{-7}\) (Figure 3), which is compatible with the truncation error of the numerical discretisation. Moreover, the error growth remains bounded over the full integration interval and does not trigger spurious oscillations or visible distortions of the large-scale wave structures. This behaviour indicates that, despite the lack of a geostrophically and hydrostatically well-balanced discretisation, the proposed high-order IMEX-DG formulation preserves geostrophic equilibrium with satisfactory accuracy and robustness also over long simulation times and in large-scale rotating configurations (M. Baldauf, personal communication).

4.3 3D flow over a bell-shaped hill↩︎

In this section, we consider fully three-dimensional atmospheric dynamics simulations with the IMEX-DG scheme. These experiments enable the assessment of the interaction between stratification, orographic forcing, and rotation in more realistic and representative dynamical configurations. Two test cases of flow over orography are considered: an existing medium-steep benchmark [5], [13], [15], [16], used to analyse the parallel performance of the solver, and a new, larger-scale atmospheric configuration designed to investigate the behaviour of the method in a regime where the effect of rotation is more prominent. Both configurations are simulated with and without rotation in order to assess the impact of the Coriolis term and verify the robustness of the solver in fully three-dimensional flows. For both tests, the computational domain is a hexahedral box where the lower boundary is deformed by a smooth bell-shaped mountain centered at \((x_c,y_c)\) km in the horizontal plane, with height \[h(x,y)=h_0\left(1+\frac{(x-x_c)^2+(y-y_c)^2}{a_c^2}\right)^{-1}, \label{eq:bell95mountain}\tag{14}\] where \(h_0\) denotes the maximum mountain height and \(a_c\) controls the horizontal half-width [13]. The mesh is generated on the physical domain and mapped to the terrain-following coordinates through the Gal-Chen transformation [38]. Periodic boundary conditions are imposed in the horizontal directions, while solid-wall conditions are applied at the bottom boundary. At the upper boundary, non-reflecting conditions are employed in order to reduce artificial wave reflections. As in [13], [14], to further suppress spurious reflections of outgoing gravity waves, Rayleigh damping (sponge) layers are introduced near the upper boundary and in lateral buffer regions. The damping coefficient is set to zero below a prescribed height \(z_B\) and increases smoothly toward the top boundary \(z_T\) according to \[\begin{align} &\lambda(z)=\bar{\lambda}\sin^2\left[\frac{\pi}{2}\left(\frac{z-z_B}{z_T-z_B}\right)\right] \qquad &&\text{for } z\ge z_B, \\ &\lambda(z)=0 \qquad &&\text{for } z<z_B . \end{align} \label{eq:damping}\tag{15}\] An analogous formulation is employed in the lateral damping regions.

Benchmark test↩︎

The first domain considered is a \(\Omega=\SI[parse-numbers=false]{(0,60) \times (0,40) \times (0,16)}{\kilo\meter\cubed}\) box. The center of the mountain is at \((x_c,y_c)=\SI[parse-numbers=false]{(30,20)}{\kilo\meter}\) and \(h_0=\SI{400}{\meter}\), \(a_c=\SI{1000}{\meter}\). Damping is applied in the first and last \(\SI{20}{\kilo\meter}\) in the \(x\) direction, in the first and last \(\SI{10}{\kilo\meter}\) in the \(y\) direction, and in the uppermost \(\SI{6}{\kilo\meter}\) in the \(z\) direction, with \(\bar{\lambda} \Delta t=1.2\). The flow is initialized from a stationary hydrostatic background state. The expression for the background potential temperature is given by \(\bar \theta=\theta_{\mathrm{ref}}\exp(\frac{N^2}{g}z)\). The buoyancy frequency is set to \(N=\SI{0.01}{\per\second}\), corresponding to a nonhydrostatic regime, and \(\theta_{\mathrm{ref}}=\SI{293.15}{\kelvin}\). When rotation is active, the background is completed with a geostrophically balanced horizontal pressure gradient. The background pressure is decomposed into a height-dependent hydrostatic contribution and a geostrophic correction accounting for rotation: \[\label{eq:init95cond} \begin{align} &p_0=p_{0,h}(z)+p_{0,g}(x,y), \\ & p_{0,h}= p_{\text{ref}}(\pi)^{\frac{\gamma}{\gamma-1}} ,\quad \pi =1+ \frac{g^2}{c_p T_{\text{ref}}N^2} \left[\exp{\left(-\frac{N^2}{g} z\right)}-1\right],\\ & p_{0,g}(x,y)= (x-x_0)\,\rho_0 f\, U\sin\alpha-(y-y_0)\,\rho_0 f\, U\cos\alpha, \end{align}\tag{16}\] where \(\pi\) is the Exner pressure. In the non-rotating case (\(f=0\)), the geostrophic contribution vanishes and \(p_0\) reduces to the hydrostatic profile. A uniform horizontal wind is imposed, \(\mathbf{u}_0=U_0(\cos{\alpha},\sin{\alpha},0),\) with \(\alpha=0\) (flow aligned with the \(x\) direction) and \(U_0=\SI{10}{\meter\per\second}\). We use polynomial degree \(r=4\) on a uniform mesh of \(30 \times 20 \times 8 = 4800\) elements, yielding an effective resolution of \(\SI{500}{\meter}\) in all directions, and \(6 \times 10^5\) degrees of freedom for each scalar variable. The time step is \(\SI{8}{\second}\) and the final time is \(T_f=8\) h. The corresponding acoustic and advective Courant numbers are approximately \(C=0.21\) and \(C_{\text{adv}}=5.5\), respectively. We consider two simulations: a non-rotating run (\(f=0\)) and a rotating run with \(f=\SI[parse-numbers=false]{10^{-4}}{\per\second}\).

In both configurations, the orographic forcing generates a wake and a downstream gravity-wave pattern with comparable peak amplitudes, indicating that the inclusion of the Coriolis terms does not introduce unphysical effects (Figure 4). However, in the rotating case, the downstream structure loses the axial symmetry of the non-rotating case, consistently with the effect of rotation. This is evident in the horizontal slices of the \(v\) component (panel (d) in Figure 4), and even more clearly in the vertical sections (third and fourth row in Figure 4): at the same downstream distance from the mountain, the non-rotating case exhibits contours that are mirror images with opposite sign across the wake, whereas the rotating case clearly shows a difference in the contour shapes, revealing the induced asymmetry. Since the computational domain is relatively limited in horizontal extent, the rotational effects remain moderate. Overall, this benchmark provides a first validation of the implementation of the rotating terms in the IMEX-DG scheme in a fully three-dimensional setting.

To obtain cleaner diagnostics while keeping the computational cost manageable, two additional simulations are performed using a non-conforming mesh. The non-conforming discretization allows for local refinement near the orography, where the strongest dynamics and wave generation occur, while maintaining a coarser resolution elsewhere [13].

The base mesh has an effective horizontal resolution of \(1\) km and is locally refined by two additional levels around the mountain, reaching a finest resolution of approximately \(250\) m. Refinement is applied in the region \(x\in[20,40]\) km and \(y\in[10,30]\) km; vertically, the first refinement level is used for \(z\le7\) km and the second for \(z\le2\) km (Figure 5), following the same approach as in [13], [14].

With this non-conforming mesh, the total number of active cells is 1594 (about one third of the conforming mesh), and the number of DOF per scalar variable is 199 250, i.e. also about one third of the uniform-grid configuration.

The simulation using the non-conforming mesh reproduces the same large-scale wave dynamics observed with the uniform grid, while providing smoother and better-resolved structures near the mountain and in the near-wake region, where the strongest gradients occur (Figure 6). The comparison between rotating and non-rotating cases leads to the same qualitative conclusions as before, with rotation producing only a weak asymmetry in the downstream wave pattern. Overall, the non-conforming discretization improves the local accuracy around the orography at a significantly reduced computational cost, while also enabling targeted vertical refinement near the terrain.

Since this benchmark coincides with the configuration considered in [15], [16] for the non-rotating case, we repeat the strong-scaling analysis in the presence of rotation in order to assess the impact of the rotational term on the parallel performance of the solver and to directly compare the results with the previously published results. We consider both a uniform mesh composed of \(120\times80\times32 = 307200\) elements with polynomial degree \(r=4\), and a non-conforming mesh consisting of \(204816\) elements. The strong-scaling analysis evaluates the wall-clock time required to perform the simulations while keeping the computational workload, namely the spatial resolution, fixed and progressively increasing the available computational resources. In particular, the simulations are executed using 1, 2, 4, 8, and 16 full MeluXina CPU nodes, each providing 128 cores. Figure 7 reports the comparison between the previously published strong-scaling results obtained in the absence of rotation and the corresponding results including the rotational term.

Overall, the introduction of rotation does not significantly deteriorate the parallel efficiency of the solver. The speedup curves obtained with and without rotation exhibit a very similar trend for both the uniform and non-conforming meshes, confirming that the additional rotational contribution introduces only a limited overhead in the parallel implementation. In particular, the solver maintains a close-to-ideal scaling behaviour up to 1024 cores, while a moderate degradation is observed at 2048 cores, where communication costs and load-balancing effects become more relevant. Moreover, the non-conforming configuration shows a scaling behaviour comparable to the uniform mesh case, demonstrating the robustness of the AMR framework also in the presence of rotational effects.

Larger-scale test↩︎

In order to better assess the impact of rotation on the flow dynamics and to verify the correct preservation of geostrophic balance by the IMEX-DG solver, we designed a new test of 3D flow over orography on a larger computational domain. The enlarged domain allows rotational effects to develop over larger spatial scales, reducing the influence of the lateral damping regions and making the differences between the rotating and non-rotating configurations more evident. The computational domain is a \(\Omega=\SI[parse-numbers=false]{(0,240)\times(0,240)\times(0,30)}{\kilo\meter\cubed}\) box. The mountain is centered at \((x_c,y_c)=\SI[parse-numbers=false]{(120,120)}{\kilo\meter}\) and is defined by 14 with \(h_0=\SI{100}{\meter}\) and \(a_c=\SI{10}{\kilo\meter}\). Damping layers are applied in the first and last \(\SI{60}{\kilo\meter}\) of both horizontal directions and in the uppermost \(\SI{14}{\kilo\meter}\) of the domain, using expressions 15 with \(\bar{\lambda}\Delta t = 1.2\).
The flow is initialized from a stationary hydrostatic background state. In this configuration, the Exner pressure profile is prescribed as \[\pi=\left(\frac{p_0}{\bar{p}}\right)^{\frac{\gamma-1}{\gamma}} =\exp\left(-\frac{g}{c_p\bar{T}}z\right),\] corresponding to an isothermal hydrostatic atmosphere. As in the previous test case, a uniform horizontal background wind is imposed, \(\mathbf{u}_0=U_0(\cos\alpha,\sin\alpha,0)\), with \(\alpha=0\) and \(U_0=\SI{10}{\meter\per\second}\). The reference potential temperature is set to \(\theta_{\mathrm{ref}}=\SI{250}{\kelvin}\). The buoyancy frequency is equal to \(N=\SI{0.02}{\per\second}\). The resulting non-dimensional parameter \(\frac{Na_c}{U_0}=20\gg1\), so, according to the classification in [39], this configuration falls in the hydrostatic regime.

The simulations are performed using polynomial degree \(r=4\) on uniform mesh composed of \(50 \times 50 \times 50\) elements. The time step is \(\Delta t=\SI{4}{\second}\) and the final time is set to \(T_f=\SI{8}{\hour}\). The corresponding acoustic and advective Courant numbers are approximately \(C=1.29\) and \(C_{\text{adv}}=0.04\), respectively.

As in the smaller-domain configuration, the interaction between the stratified background flow and the orography generates a wake and a downstream gravity-wave pattern with comparable amplitudes in both the rotating and non-rotating simulations, confirming that the inclusion of the Coriolis term does not introduce spurious or non-physical effects (Figure 8). However, in the present larger-scale configuration the influence of rotation is considerably more pronounced. While the non-rotating solution maintains an approximately symmetric structure with respect to the wake axis (left panels), the rotating case develops a clear asymmetry in both the horizontal and vertical velocity fields (right panels). This behaviour is particularly evident in the horizontal slices of the \(v\) component shown in panels (c) and (d) of Figure 8, where the rotating solution exhibits a visible lateral distortion of the wave pattern. The effect becomes even clearer in the vertical \((x-z)\) sections reported in the third and fourth rows of Figure 8. These slices are taken at fixed \(y\) locations and show that, at comparable downstream distances from the mountain, the non-rotating case retains nearly symmetric contour structures with opposite signs across the wake, whereas the rotating simulation displays tilted and distorted wave fronts induced by the Coriolis acceleration. The larger horizontal extent of the domain allows the rotational effects to develop over longer spatial and temporal scales, making the associated asymmetry significantly more visible than in the smaller-domain test. At the same time, the numerical solution remains stable and free from noticeable spurious oscillations, indicating that the proposed rotating IMEX-DG formulation preserves robustness also in larger fully three-dimensional atmospheric configurations.

5 Conclusions↩︎

We have presented an extension of the IMEX-DG framework for compressible atmospheric flows proposed in [12], [13] to rotating configurations including the Coriolis force. The rotational contribution has been incorporated consistently within the IMEX decomposition and matrix-free DG discretisation and implemented in the deal.II framework. Two alternative treatments of the rotational term within the nonlinear fixed-point formulation have been introduced and compared. A relevant aspect of the proposed work is the efficient matrix-free treatment of the inverse of the operator \((\mathbf{A}+\mathbf{R})\), obtained by exploiting the discrete structure of the rotational operator.

Simulations of standard and newly introduced rotating benchmarks demonstrate that the proposed formulation accurately reproduces inertia-gravity wave dynamics and rotating flow over orography in both two- and three-dimensional configurations. The numerical results show that the method is able to capture the asymmetry induced by rotation without introducing spurious oscillations or loss of stability, even in long integration time and large-scale configurations. Convergence studies confirm that the inclusion of rotation preserves the expected accuracy properties of the original framework, while the comparison between the two rotational treatments indicates that they lead to nearly equivalent numerical behaviour for the considered test cases.

The numerical experiments further show that the rotating extension remains robust on both uniform and non-conforming adaptive meshes and does not significantly affect the parallel scalability of the original solver, as measured in strong scaling experiments on up to 2048 cores. Future developments will focus on planetary-scale three dimensional configurations, on detailed comparisons with exactly well-balanced discretisations for longer lead times, and on extensions to spherical geometries.

Acknowledgments↩︎

This manuscript expands on the first author’s Master’s thesis for a double degree program with Politecnico di Milano and Technical University of Denmark [40]. G.O. is part of the INdAM-GNCS National Research Group. We gratefully acknowledge Prof. Antonella Abbà at Politecnico di Milano for useful comments on this work. We thank Dr. Michael Baldauf for kindly making available the code to compute the analytical solution in [36]. We acknowledge the EuroHPC Joint Undertaking for awarding this project access to the EuroHPC supercomputer MeluXina through EuroHPC Benchmark Access calls #EHPC-BEN-2026B01-008 and #EHPC-BEN2026B04-091. Generative AI tools have been use exclusively to improve language clarity, without contributing to the scientific content of the work.

References↩︎