# Impact of curved elements for flows over orography with a Discontinuous Galerkin scheme

Giuseppe Orlando, Tommaso Benacchio, Luca Bonaventura

► To cite this version:

Giuseppe Orlando, Tommaso Benacchio, Luca Bonaventura. Impact of curved elements for flows over orography with a Discontinuous Galerkin scheme. 2024. hal-04547681

HAL Id: hal-04547681

https://hal.science/hal-04547681

Preprint submitted on 16 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Impact of curved elements for flows over orography with a Discontinuous Galerkin scheme

Giuseppe Orlando, Tommaso Benacchio, Luca Bonaventura

► To cite this version:

Giuseppe Orlando, Tommaso Benacchio, Luca Bonaventura. Impact of curved elements for flows over orography with a Discontinuous Galerkin scheme. Journal of Computational Physics, 2024, 519, pp.113445. 10.1016/j.jcp.2024.113445. hal-04547681v3

HAL Id: hal-04547681

https://hal.science/hal-04547681v3

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Impact of curved elements for flows over orography with
a Discontinuous Galerkin scheme

Giuseppe Orlando$^{(1)}$,
Tommaso Benacchio$^{(2)}$, Luca Bonaventura$^{(3)}$

(1) CMAP, CNRS, École Polytechnique, Institute Polytechnique de Paris
Route de Saclay, 91120 Palaiseau, France
giuseppe.orlando@polytechnique.edu

(2) Weather Research, Danish Meteorological Institute
Sankt Kjelds Plads 11, 2100 Copenhagen, Denmark
tbo@dmi.dk

(3) Dipartimento di Matematica, Politecnico di Milano
Piazza Leonardo da Vinci 32, 20133 Milano, Italy
luca.bonaventura@polimi.it

**Keywords:** Numerical Weather Prediction, Discontinuous Galerkin methods, High-order mapping, Flows over orography, Curved boundary.
Impact of curved elements for flows over orography with
a Discontinuous Galerkin scheme

Giuseppe Orlando$^{(1)}$,
Tommaso Benacchio$^{(2)}$, Luca Bonaventura$^{(3)}$

(1) CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris
Route de Saclay, 91120 Palaiseau, France
giuseppe.orlando@polytechnique.edu

(2) Weather Research, Danish Meteorological Institute
Sankt Kjelds Plads 11, 2100 Copenhagen, Denmark
tbo@dmi.dk

(3) Dipartimento di Matematica, Politecnico di Milano
Piazza Leonardo da Vinci 32, 20133 Milano, Italy
luca.bonaventura@polimi.it

**Keywords:** Numerical Weather Prediction, Discontinuous Galerkin methods, High-order mapping, Flows over orography, Curved boundary.
Abstract

We present a quantitative assessment of the impact of high-order mappings on the simulation of flows over complex orography. Curved boundaries were not used in earlier numerical methods, whereas they are employed nowadays to an increasing extent in combination with high-order methods, such as the Finite Element Method (FEM) and the Spectral Element Method (SEM). We consider here a specific Discontinuous Galerkin (DG) method implemented in the framework of the *deal.II* library, which natively supports high-order mappings. A number of numerical experiments based on classical benchmarks over idealized orographic profiles demonstrate the positive impact of curved boundaries on the accuracy of the results. These findings are also supported by results of the application of this approach to non-smooth and realistic orographic profiles.
Abstract

We present a quantitative assessment of the impact of high-order mappings on the simulation of flows over complex orography. Curved boundaries were not used in early numerical methods, whereas they are employed to an increasing extent in state of the art computational fluid dynamics codes, in combination with high-order methods, such as the Finite Element Method and the Spectral Element Method. Here we consider a specific Discontinuous Galerkin (DG) method implemented in the framework of the deal.II library, which natively supports high-order mappings. A number of numerical experiments based on classical benchmarks over idealized orographic profiles demonstrate the positive impact of curved boundaries on the accuracy of the results, with no significantly adverse effect on the computational cost of the simulation. These findings are also supported by results of the application of this approach to non-smooth and realistic orographic profiles.
# 1 Introduction

High-order numerical methods are employed to an increasing extent for many relevant physical problems, in particular for applications in Computational Fluid Dynamics (CFD), see e.g. [6, 16, 20] among many others. Such methods must be supplemented by high-order approximations of curved geometries to maintain their high-order accuracy [28] and achieve optimal convergence rates [11]. For instance, considering Dirichlet boundary conditions, any finite element method on curved domains is at most second order accurate [48, 49], unless a curved boundary element is adopted [11]. Thus, the use of curved meshes takes on considerable importance for many realistic applications and the treatment of boundary conditions for curved boundaries is an active research area, see e.g. the recent contributions [9, 10]. While designing numerical schemes in combination with curved boundaries posed limitations for early numerical methods, a number of approaches to deal naturally with curved geometries have been proposed for high-order methods over the last fifty years. The most popular approach to build a curvilinear element method relies on the iso-parametric approximation of the curved boundary, which, since the seminal contributions [24, 57], has widely been employed in the literature, see among many others [15, 18, 39, 50, 51]. The most advanced developments use rational B-spline or NURBS approximations in the so-called iso-geometric analysis (IGA) framework [29]. An alternative to curvilinear elements is to improve the treatment of the boundary conditions, taking into account the features of the “real” curved geometry on a straight faced mesh [34, 53]. However, this method can be formulated only for slip-wall conditions and for 2D geometries. We also mention the approach described e.g in [1], in which a sub-tessellation is considered to evaluate integrals on straight faced meshes. Finally, recent developments employ a reconstruction on-site data (ROD) approach [12, 13], using a least-squares method to handle several constraints from scattered mean values associated to the elements, while yet another approach employs the so-called Shifted Boundary Method (SBM), see [9], to which we refer for all the details.

We focus here on the Discontinuous Galerkin (DG) method, which combines high-order accuracy and flexibility [22]. The iso-parametric approximation and IGA have been widely used in conjunction with DG schemes [6, 27, 33]. However, these approaches are particularly expensive, both in terms of computational cost and memory overhead. A simpler treatment of curved elements in the context of DG methods uses instead high-order mappings from the straight-sided reference element to each curved element [35, 52], leading to a non-trivial expression for the determinant of the mapping Jacobian (see Section 3).

While the accuracy loss due to the use of low-order mappings can be clearly identified and analyzed when dealing with smooth boundaries, the use and impact of high-order mappings in presence of irregular and non-smooth boundaries is less straightforward. This aspect is especially important when dealing with applications to numerical weather prediction (NWP) and climate simulations, in which the lower boundary is defined by irregular orographic data. In NWP and climate models it is often necessary to filter the orographic data to avoid or reduce spurious numerical effects [54]. Subgrid-scale orographic drag parametrization is then employed to compensate the insufficient resolution of orographic features [38, 44]. On the other hand, increasing the resolution without parametrization is beneficial to improve the description of atmospheric pro-
# 1 Introduction

High-order numerical methods are employed to an increasing extent for many relevant physical problems, in particular for applications in Computational Fluid Dynamics (CFD), see e.g. [6, 19, 24], among many others. Moreover, in recent years, a general effort to design numerical methods allowing for more general element shapes, such as polygonal/polyhedral elements, is witnessed, see, e.g., [9, 14, 23], and the references therein. However, such methods must be supplemented by high-order approximations of curved geometries to maintain their high-order accuracy [33] and achieve optimal convergence rates [13]. For instance, considering Dirichlet boundary conditions, any finite element method on curved domains is at most second order accurate [54, 55], unless a curved boundary element is adopted [13]. Thus, the use of curved meshes takes on considerable importance for many realistic applications, and the treatment of boundary conditions for curved boundaries is an active research area, see e.g. the recent contributions [11, 12]. While the use of numerical schemes in combination with curved boundaries posed limitations for early numerical methods, a number of approaches to deal naturally with curved geometries have been proposed for high-order methods over the last fifty years.

A very popular approach to build a curvilinear element method relies on the isoparametric approximation of the curved boundary which, since the seminal contributions [29, 63], has widely been employed in the literature, see among many others [18, 21, 44, 56, 57]. Advanced developments use rational B-spline or NURBS approximations in the so-called iso-geometric analysis (IGA) framework [34]. An alternative to curvilinear elements is to improve the treatment of the boundary conditions, taking into account the features of the real curved geometry on a straight faced mesh [39, 59]. However, this method can be formulated only for slip-wall conditions and for 2D geometries. We also mention the approach described, e.g, in [1], in which a sub-tessellation is considered to evaluate integrals on straight faced meshes. Finally, recent developments employ a reconstruction off-site data approach [15, 16], using a least-squares method to handle several constraints imposed by scattered mean values associated to the elements, while yet another approach employs the so-called Shifted Boundary Method, see [11], to which we refer for all the details.

We focus here on the Discontinuous Galerkin (DG) method, which combines high-order accuracy and flexibility [27]. A theoretical analysis on the applicability of the DG method on essentially arbitrary shaped-elements, including curved polygonal/polyhedral elements, was presented in [10]. The iso-parametric approximation and IGA have been widely used in conjunction with DG schemes [6, 32, 38]. However, these approaches are particularly expensive, both in terms of computational cost and memory overhead. A simpler treatment of curved elements in the context of DG methods uses instead high-order mappings from the straight-sided reference element to each curved element [40, 58], leading to a non-trivial expression for the determinant of the mapping Jacobian (see Section 3).

While the accuracy loss due to the use of low-order mappings can be clearly identified and analyzed when dealing with smooth boundaries, the use and impact of high-order mappings in presence of irregular and non-smooth boundaries is less straightforward. This aspect is especially important when dealing with applications to numerical weather prediction (NWP) and climate simulations, in which the lower boundary is
cesses over complex orography and forecast skill [19, 30]. Moreover, the positive impact of parameterizations decreases as the model resolution increases. Hence, it is important to assess the performance of high-order mappings also in the case of non-smooth orography and to analyse the interplay of those mappings with standard filtering procedures. Atmospheric motion over idealized orography is the focus of a number of popular benchmarks proposed and analyzed in the NWP literature [7, 31, 32, 46]. The orography is typically described by non-linear smooth analytical profiles and therefore the computational domain is characterized by a curved boundary. Nevertheless, to the best of our knowledge, curved elements have not been employed for such benchmarks and no survey of their impact on numerical results is available.

In this work, we assess the impact of high-order mappings on the accuracy of numerical simulations of atmospheric flow over both idealized and real orography. For the purpose of this study, we use the IMEX-DG scheme proposed in [40, 41], to which we refer for a complete analysis and description of the method. The method was validated for the simulations of weakly compressible atmospheric flows in [42, 43]. The solver is implemented in the framework of the open-source library *deal.II* [2, 5], which natively supports high-order mappings. We show that, for smooth profiles, using high-order mappings leads to more accurate results than using linear mapping. In addition, high-order mappings are found to provide analogous results to those obtained with linear mappings at higher resolution, meaning that the use of high-order mappings acts as a sort of sub-tessellation. Next, we modify the customary benchmarks by adding a non-differentiable perturbation. We compare the results of very high resolution simulations performed with linear mappings that resolve well the irregular profile, to those of lower resolution simulations carried out with both low- and high-order mappings. High-order mappings are found to provide more accurate results also in this case. We then assess the impact of two common filtering procedures on the results. Finally, we consider the use of high-order mappings for a realistic complex orography described by a set of point data.

The paper is structured as follows. The model equations are briefly presented in Section 2. A short review of the use of high-order mappings for DG schemes is reported in Section 3. Relevant numerical simulations showing the impact of high-order mappings and curved elements are presented in Section 4. Finally, some conclusions are reported in Section 5.

## 2 The model equations

The mathematical model consists of the fully compressible Euler equations of gas dynamics in conservation form under the influence of gravity, completed with the ideal gas law [42, 43], in a vertical $(x, z)$ slice domain. Let $\Omega \subset \mathbb{R}^d, d = 2$, be a connected open bounded set with a possibly curved boundary $\partial\Omega$, with $\mathbf{x}$ denoting the spatial coordinates and $t$ the temporal coordinate. We consider the case $d = 2$ in this work, but the analysis and the considerations outlined in the following sections can be straightforwardly extended to $d = 3$. We consider the unsteady compressible Euler
defined by irregular orographic data. Indeed, in NWP and climate models it is often
necessary to filter the orographic data to avoid or reduce spurious numerical effects
[60]. Subgrid-scale orographic drag parametrizations are then employed to compensate
the insufficient resolution of orographic features [43, 49]. On the other hand, increasing
the resolution without parametrization is beneficial to improve the description of atmo-
spheric processes over complex orography and forecast skill [22, 35]. The importance of
an accurate description of orography is also highlighted in studies like [50], where the
impact of improved orographic resolution is assessed via idealized climate simulations.
Moreover, the positive impact of parametrizations decreases as the model resolution
increases. Hence, it is important to assess the performance of high-order mappings also
in the case of non-smooth orography and to analyse the interplay of those mappings
with standard filtering procedures. Atmospheric motion over idealized orography is
the focus of a number of popular benchmarks proposed and analyzed in the NWP
literature [7, 36, 37, 52]. The orography is typically described by non-linear smooth
analytical profiles and therefore the computational domain is characterized by a curved
boundary. Nevertheless, to the best of our knowledge, curved elements have not been
employed for such benchmarks and no survey of their impact on numerical results is
available.

In this work, we assess the impact of high-order mappings on the accuracy of
numerical simulations of atmospheric flow over both idealized and real orography. For
the purpose of this study, we use the IMEX-DG scheme proposed in [45, 46], to which
we refer for a complete analysis and description of the method. The method was
validated for the simulations of weakly compressible atmospheric flows in [47, 48].
The solver is implemented in the framework of the open-source library deal.II [2, 5],
which natively supports high-order mappings. We show that, for smooth orography
profiles, using high-order mappings leads to more accurate results than using linear
mapping. In addition, high-order mappings are found to provide analogous results
to those obtained with linear mappings at higher resolution, meaning that the use of
high-order mappings acts as a sort of sub-tessellation. Next, we modify the customary
smooth benchmarks by adding a non-differentiable perturbation. We compare the
results of very high resolution simulations using linear mappings that resolve well the
irregular profile to those of lower resolution simulations carried out with both low- and
high-order mappings. High-order mappings are found to provide more accurate results
also in this case. We then assess the impact of two common filtering procedures on
the results. Finally, we consider the use of high-order mappings for a realistic complex
orography described by a set of point data.

The paper is structured as follows. The model equations are briefly presented
in Section 2. A short review of the use of high-order mappings for DG schemes is
reported in Section 3. Relevant numerical simulations showing the impact of high-order
mappings and curved elements are presented in Section 4. Finally, some conclusions
are reported in Section 5.

2 The model equations

The mathematical model consists of the fully compressible Euler equations of gas
dynamics in conservation form under the influence of gravity, completed with the
equations, written in conservation form as

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) &= 0 \\
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p &= \rho \mathbf{g} \\
\frac{\partial (\rho E)}{\partial t} + \nabla \cdot [(\rho E + p) \mathbf{u}] &= \rho \mathbf{g} \cdot \mathbf{u},
\end{aligned}
\quad (1) $$

for $\mathbf{x} \in \Omega$, $t \in (0, T_f]$, endowed with suitable initial and boundary conditions. Here $T_f$ is the final time, $\rho$ is the density, $\mathbf{u}$ is the fluid velocity, $p$ is the pressure, and $\otimes$ denotes the tensor product. Moreover, $\mathbf{g} = -g\mathbf{k}$ represents the acceleration of gravity, with $g = 9.81 \text{ m s}^{-2}$ and $\mathbf{k}$ denoting the upward pointing unit vector in the standard Cartesian frame of reference. The total energy $\rho E$ can be rewritten as $\rho E = \rho e + \rho k$, where $e$ is the internal energy and $k = \frac{1}{2}|\mathbf{u}|^2$ is the kinetic energy. We also introduce the specific enthalpy $h = e + \frac{p}{\rho}$ and we notice that one can rewrite the energy flux as

$$ (\rho E + p) \mathbf{u} = \left( e + k + \frac{p}{\rho} \right) \rho \mathbf{u} = (h + k) \rho \mathbf{u}. \quad (2) $$

The above equations are complemented by the equation of state for ideal gases, given by

$$ p = \rho RT, \quad (3) $$

with $R$ being the specific gas constant. For later reference, we define the Exner pressure $\Pi$ as

$$ \Pi = \left( \frac{p_0}{p} \right)^{\frac{\gamma-1}{\gamma}}, \quad (4) $$

with $p_0 = 10^5 \text{ Pa}$ being a reference a pressure and $\gamma$ denoting the specific heats ratio. We take the specific heats ratio $\gamma = 1.4$ and the specific gas constant $R = 287 \text{ J kg}^{-1} \text{ K}^{-1}$.

## 3 High-order mappings for DG schemes

We consider a decomposition of the domain $\Omega$ into a family of quadrilaterals $\mathcal{T}_h$ and denote each element by $K$. The most classical approach for the treatment of curved boundaries in the finite element method is an iso-parametric approximation, in which both the geometry and the solution are approximated by some high-order polynomial functions [24, 57]. In particular, rational B-spline or NURBS approximations lead to the iso-geometric analysis (IGA) [29]. Even though its effectiveness has been shown in several contributions, see e.g. [39, 45, 50, 51], this approach can be expensive in terms of computational cost and memory overhead, because of the integration on curvilinear elements, especially when the boundary geometry is represented by very high-order polynomials [55, 56].

We consider here a simpler approach based on a local polynomial interpolation of the curved geometry. We define as mappings the transformations from the reference element to elements in the physical space. Let $\hat{\mathbf{x}}$ be a point in the reference element, namely the unit square in 2D. Hence, each mapping is a function $\mathbf{F}_K$ such that
ideal gas law [47, 48], in a vertical $(x, z)$ slice domain. Let $\Omega \subset \mathbb{R}^d, d = 2$, be a connected open bounded set with a possibly curved boundary $\partial\Omega$, with $\mathbf{x}$ denoting the spatial coordinates and $t$ the temporal coordinate. We consider the case $d = 2$ in this work, but the analysis and the considerations outlined in the following sections straightforward extend to $d = 3$. No ad hoc two-dimensional numerical procedure or code was employed for the numerical results presented in this work. We consider the unsteady compressible Euler equations, written in conservation form as

$$ 
\begin{aligned}
& \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \\
& \frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p = \rho \mathbf{g} \\
& \frac{\partial (\rho E)}{\partial t} + \nabla \cdot [(\rho E + p) \mathbf{u}] = \rho \mathbf{g} \cdot \mathbf{u},
\end{aligned}
\quad (1) $$

for $\mathbf{x} \in \Omega, t \in (0, T_f]$, endowed with suitable initial and boundary conditions. Here $T_f$ is the final time, $\rho$ is the density, $\mathbf{u}$ is the fluid velocity, $p$ is the pressure, and $\otimes$ denotes the tensor product. Moreover, $\mathbf{g} = -g\mathbf{k}$ represents the acceleration of gravity, with $g = 9.81 \text{ m s}^{-2}$ and $\mathbf{k}$ denoting the upward pointing unit vector in the standard Cartesian frame of reference. The total energy $\rho E$ can be rewritten as $\rho E = \rho e + \rho k$, where $e$ is the internal energy and $k = \frac{1}{2}|\mathbf{u}|^2$ is the kinetic energy. We also introduce the specific enthalpy $h = e + \frac{p}{\rho}$ and we notice that one can rewrite the energy flux as

$$ (\rho E + p) \mathbf{u} = \left( e + k + \frac{p}{\rho} \right) \rho \mathbf{u} = (h + k) \rho \mathbf{u}. \quad (2) $$

The above equations are complemented by the equation of state for ideal gases, given by

$$ p = \rho RT, \quad (3) $$

with $R$ being the specific gas constant. For later reference, we define the Exner pressure $\Pi$ as

$$ \Pi = \left( \frac{p_0}{p} \right)^{\frac{\gamma-1}{\gamma}}, \quad (4) $$

with $p_0 = 10^5 \text{ Pa}$ being a reference pressure and $\gamma$ denoting the specific heats ratio. We take the specific heats ratio $\gamma = 1.4$ and the specific gas constant $R = 287 \text{ J kg}^{-1} \text{ K}^{-1}$.

## 3 High-order mappings for DG schemes

We consider a tessellation of the domain $\Omega$ into a family of quadrilaterals $\mathcal{T}_h$ and denote each element by $K$. The most classical approach for the treatment of curved boundaries in the finite element method is the iso-parametric approximation, in which both the geometry and the solution are approximated by some high-order polynomial functions [29, 63]. In particular, rational B-spline or NURBS approximations lead to the iso-geometric analysis (IGA) [34]. Even though its effectiveness has been shown in several contributions, see e.g. [44, 51, 56, 57], this approach can be expensive in terms of computational cost and memory overhead, because of the integration on curvilinear
$\mathbf{x} = \mathbf{F}_K(\hat{\mathbf{x}})$. We assume that the function $\mathbf{F}_K$ is invertible, meaning that each point of $\hat{\mathbf{x}}$ of the reference element corresponds to one and only one point $\mathbf{x}$ in the real element $K$ [26]. A relevant quantity is the Jacobian of the transformation $\mathbf{J}_K(\hat{\mathbf{x}}) = \nabla \mathbf{F}_K(\hat{\mathbf{x}})$. Requiring that $\mathbf{F}_K$ is invertible is equivalent to assuming $\mathbf{J}_K(\hat{\mathbf{x}})$ has nonzero determinant [26]. For a generic function $\phi(\mathbf{x})$, which can be real-valued, vector-valued or tensor-valued, the mapping acts therefore as follows:

$$ \phi(\mathbf{x}) = \phi(\mathbf{F}_K(\hat{\mathbf{x}})) = \hat{\phi}(\hat{\mathbf{x}}). \qquad (5) $$

The DG method is characterized by integrals over elements and over faces that share two elements (see e.g. [22] for a general presentation of the method). Using a simple change of variables, it is possible to express integrals over a cell $K$ as integral over the reference element $\hat{K}$. More specifically, the following relation holds [4, 26]:

$$ \int_K \phi(\mathbf{x}) dx = \int_{\hat{K}} \hat{\phi}(\hat{\mathbf{x}}) |\det \mathbf{J}_K(\hat{\mathbf{x}})| d\hat{\mathbf{x}}, \qquad (6) $$

where $\det \mathbf{J}_K(\hat{\mathbf{x}})$ is the determinant of the Jacobian and $d\hat{\mathbf{x}}$ is the volume form in the reference element. Analogous considerations hold for face integrals. Note also that the transformation of differential forms, such as gradients of scalar functions, denoted by $\mathbf{v}$, and gradients of vector fields, denoted by $\mathbf{T}$, follows the general rule

$$ \mathbf{v}(\mathbf{x}) = \mathbf{A}(\hat{\mathbf{x}}) \hat{\mathbf{v}}(\hat{\mathbf{x}}) \qquad \mathbf{T}(\mathbf{x}) = \mathbf{A}(\hat{\mathbf{x}}) \hat{\mathbf{T}}(\hat{\mathbf{x}}) \mathbf{B}(\hat{\mathbf{x}}). \qquad (7) $$

The differential forms $\mathbf{A}$ and $\mathbf{B}$ are determined by the kind of object being transformed, as discussed in [4] (see also [3] for an analysis of finite element differential forms on curvilinear cubic meshes). For the sake of clarity, we report the transformation for the gradient of a shape function $\phi$, which reads as follows:

$$ \nabla \phi(\mathbf{x}) = \mathbf{J}_K^{-1}(\hat{\mathbf{x}}) \hat{\nabla} \hat{\phi}(\hat{\mathbf{x}}). \qquad (8) $$

The integral in (6) is then approximated by a quadrature rule, so as to obtain

$$ \int_K \phi(\mathbf{x}) dx \approx \sum_q \hat{\phi}(\hat{\mathbf{x}}_q) |\det \mathbf{J}_K(\hat{\mathbf{x}}_q)| w_q, \qquad (9) $$

where $w_q$ are the weights of the chosen quadrature formula. A visual impression of orography representation using high-order and low-order mapping is presented in Figure 1. Note that the transformation depends on the location of the physical element, therefore it is different for each element.

# 4 Numerical results

The technique outlined in Section 3 is now employed in a number of two-dimensional benchmarks concerning flows over orography. Several analytical profiles have been considered in literature for the bottom boundary. We focus here on two of them:
elements, especially when the boundary geometry is represented by very high-order polynomials [61, 62].

We consider here a simpler approach based on the local polynomial interpolation of the curved geometry. We define as mappings the transformations from the reference element to elements in the physical space. Let $\hat{\mathbf{x}}$ be a point in the reference element, namely the unit square $[0, 1] \times [0, 1] \subset \mathbb{R}^2$. Hence, each mapping is a function $\mathbf{F}_K$ such that $\mathbf{x} = \mathbf{F}_K(\hat{\mathbf{x}})$. We assume that the function $\mathbf{F}_K$ is invertible. A relevant quantity is the Jacobian of the transformation $\mathbf{J}_K(\hat{\mathbf{x}}) = \nabla \mathbf{F}_K(\hat{\mathbf{x}})$. Requiring that $\mathbf{F}_K$ is invertible is equivalent to assuming $\mathbf{J}_K(\hat{\mathbf{x}})$ has nonzero determinant [31]. For a generic function $\phi(\mathbf{x})$, which can be real-valued, vector-valued or tensor-valued, the mapping acts therefore as follows:

$$ \phi(\mathbf{x}) = \phi(\mathbf{F}_K(\hat{\mathbf{x}})) = \hat{\phi}(\hat{\mathbf{x}}). \qquad (5) $$

The DG method is characterized by integrals over elements and over faces that share two elements (see e.g. [27] for a general presentation of the method). Using a simple change of variables, it is possible to express integrals over a cell $K$ as integral over the reference element $\hat{K}$. More specifically, the following relation holds [4, 31]:

$$ \int_K \phi(\mathbf{x}) dx = \int_{\hat{K}} \hat{\phi}(\hat{\mathbf{x}}) |\det \mathbf{J}_K(\hat{\mathbf{x}})| d\hat{\mathbf{x}}, \qquad (6) $$

where $\det \mathbf{J}_K(\hat{\mathbf{x}})$ is the determinant of the Jacobian and $d\hat{\mathbf{x}}$ is the volume form in the reference element. Analogous considerations hold for face integrals. Note also that the transformation of differential forms, such as gradients of scalar functions, denoted by $\mathbf{v}$, and gradients of vector fields, denoted by $\mathbf{T}$, follows the general rule

$$ \mathbf{v}(\mathbf{x}) = \mathbf{A}(\hat{\mathbf{x}}) \hat{\mathbf{v}}(\hat{\mathbf{x}}) \qquad \mathbf{T}(\mathbf{x}) = \mathbf{A}(\hat{\mathbf{x}}) \hat{\mathbf{T}}(\hat{\mathbf{x}}) \mathbf{B}(\hat{\mathbf{x}}). \qquad (7) $$

The differential forms $\mathbf{A}$ and $\mathbf{B}$ are determined by the kind of object being transformed, as discussed in [4] (see also [3] for an analysis of finite element differential forms on curvilinear cubic meshes). For the sake of clarity, we report the transformation for the gradient of a shape function $\phi$, which reads as follows:

$$ \nabla \phi(\mathbf{x}) = \mathbf{J}_K^{-1}(\hat{\mathbf{x}}) \hat{\nabla} \hat{\phi}(\hat{\mathbf{x}}). \qquad (8) $$

The integral in (6) is then approximated by a quadrature rule, so as to obtain

$$ \int_K \phi(\mathbf{x}) dx \approx \sum_q \hat{\phi}(\hat{\mathbf{x}}_q) |\det \mathbf{J}_K(\hat{\mathbf{x}}_q)| w_q, \qquad (9) $$

where $\hat{\mathbf{x}}_q$ are the nodes and $w_q$ are the corresponding weights of the chosen quadrature formula. A visual impression of orography representation using high-order and low-order mapping is presented in Figure 1. Note that the transformation depends on the location of the physical element, therefore it is different for each element.

# 4 Numerical results

The technique outlined in Section 3 is employed in a number of two-dimensional benchmarks concerning flows over orography. Several analytical profiles for the bottom boundary have been considered in literature. We focus here on two of them:
<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison between linear mapping (left) and polynomial degree 3 mapping (right) (section of the Sierra profile in Section 4.5).</figcaption></figure>

• The versiera of Agnesi [7, 32]

$$h(x) = \frac{h_m}{1 + \left(\frac{x-x_c}{a_c}\right)^2}, \qquad (10)$$

with $h_m$ being the height of the orographic profile and $a_c$ representing the half-width.

• The five-peak mountain range profile [46]

$$h(x) = h_m \exp \left[ - \left( \frac{x}{a_c} \right)^2 \right] \cos^2 \left( \frac{\pi x}{\lambda_c} \right), \qquad (11)$$

where $\lambda_c$ is a characteristic length scale.

Wall boundary conditions are employed for the bottom boundary, whereas non-reflecting boundary conditions are needed for the top boundary and for lateral boundaries. These are achieved by applying Rayleigh damping in the regions close to the boundary, with the following coefficient [37, 42]:

$$\lambda = \begin{cases} 0, & \text{if } z < z_B \\ \frac{\pi}{\lambda} \sin^2 \left[ \frac{\pi}{2} \left( \frac{z-z_B}{z-z_T} \right) \right], & \text{if } z \ge z_B. \end{cases} \qquad (12)$$

Here, $z_B$ denotes the height at which the damping starts and $z_T$ is the top height of the considered domain. Analogous definitions apply for the two lateral boundaries. The classical Gal-Chen height-based coordinates [21] are adopted to obtain a terrain-following mesh in Cartesian coordinates:

$$z = \xi + \frac{z_T - \xi}{z_T} h(x), \qquad (13)$$

where $\xi$ is the vertical coordinate of the rectangular domain before the transformation. Transformation (13) leads to a terrain-following at bottom boundary, i.e. $\xi = 0$, and a rigid lid at top boundary, i.e. $\xi = z_T$.

A simple model for turbulent vertical diffusion for NWP applications, originally proposed in [36] and also discussed e.g. in [8], is also employed in some configurations to achieve stable solutions. The nonlinear diffusivity $\kappa$ has the form

$$\kappa \left( \frac{\partial u}{\partial z}, \frac{\partial \theta}{\partial z} \right) = l^2 \left| \frac{\partial u}{\partial z} \right| F(Ri), \qquad (14)$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Comparison between linear mapping (left) and polynomial degree 3 mapping (right) for a vertical section of the Sierra profile height in Section 4.5.</figcaption></figure>

• The versiera of Agnesi [7, 37]

$$h(x) = \frac{h_m}{1 + \left(\frac{x-x_c}{a_c}\right)^2}, \qquad (10)$$

with $h_m$ being the height of the orographic profile and $a_c$ representing its half-width.

• The five-peak mountain range profile [52]

$$h(x) = h_m \exp \left[ - \left( \frac{x}{a_c} \right)^2 \right] \cos^2 \left( \frac{\pi x}{\lambda_c} \right), \qquad (11)$$

where $\lambda_c$ is a characteristic length scale.

Wall boundary conditions are employed for the bottom boundary, whereas non-reflecting boundary conditions are needed for the top boundary and for lateral boundaries. These are achieved by applying Rayleigh damping in the regions close to the boundary, with the following coefficient [42, 47]:

$$\lambda = \begin{cases} 0, & \text{if } z < z_B \\ \frac{\pi}{\lambda} \sin^2 \left[ \frac{\pi}{2} \left( \frac{z-z_B}{z-z_T} \right) \right], & \text{if } z \ge z_B. \end{cases} \qquad (12)$$

Here, $z_B$ denotes the height at which the damping starts and $z_T$ is the top height of the considered domain. Analogous definitions apply for the two lateral boundaries. The classical Gal-Chen height-based coordinates [25] are adopted to obtain a terrain-following mesh in Cartesian coordinates:

$$z = \xi + \frac{z_T - \xi}{z_T} h(x), \qquad (13)$$

where $\xi$ is the vertical coordinate of the rectangular domain before the transformation. Transformation (13) leads to a terrain-following mesh at the bottom boundary, i.e. $\xi = 0$, and to a horizontal top boundary, i.e. $\xi = z_T$.

A simple model for turbulent vertical diffusion for NWP applications, originally proposed in [41] and also discussed e.g. in [8], is employed in some configurations to achieve stable solutions. The nonlinear diffusivity $\kappa$ has the form

$$\kappa \left( \frac{\partial u}{\partial z}, \frac{\partial \theta}{\partial z} \right) = l^2 \left| \frac{\partial u}{\partial z} \right| f(Ri), \qquad (14)$$
with *l* being a mixing length and *Ri* denoting the Richardson number given by

$$ Ri = \frac{g}{\theta_0} \frac{\frac{\partial \theta}{\partial z}}{|\frac{\partial u}{\partial z}|^2}. \qquad (15) $$

Here, $\theta_0$ is a reference temperature, while the function $F(Ri)$ is defined as

$$ F(Ri) = (1 + b|Ri|)^{\beta}, \qquad (16) $$

where

$$ \begin{cases} \beta = -2, b = 5 & \text{if } Ri > 0 \\ \beta = \frac{1}{2}, b = 20 & \text{if } Ri < 0. \end{cases} \qquad (17) $$

An important diagnostic quantity to verify if a correct orographic response is established is the vertical flux of horizontal momentum (henceforth “momentum flux”) [47]. It is defined as

$$ m(z) = \int_{-\infty}^{\infty} \bar{\rho}(z) u'(x, z) w'(x, z) dx. \qquad (18) $$

Here, $\bar{\rho}$ is the background density, whereas $u'$ and $w'$ denote the deviation from the background state of the horizontal and vertical velocity, respectively. Relation (18) is valid only for an isolated orography profile with a constant background field $\bar{u}$ [47]. Hence, we also compute the full vertical flux of horizontal momentum, given by

$$ m(z) = \int_{-\infty}^{\infty} (\bar{\rho}(z) + \rho'(x, z)) (\bar{u} + u'(x, z)) w'(x, z) dx, \qquad (19) $$

Table 1 reports the parameter values used in the test cases in this work. For time and space discretizations, we consider the IMEX-DG scheme proposed in [40, 41] and successfully employed for atmospheric dynamics in [42, 43]. We consider piecewise polynomials of degree $r = 4$ for the finite element space, unless differently stated. The solver is implemented in the framework of the deal.II library [2, 5], which natively supports high-order polynomial mappings in addition to linear mappings.

<table><caption>Table 1: Model parameters for the test cases, see main text for details. LHMW: linear hydrostatic mountain wave (Section 4.1). NLNHMW: nonlinear non-hydrostatic mountain wave (Section 4.2). Schär: Schär profile (Section 4.3). NST: non-smooth orography (Section 4.4). T-REX: Sierra profile, T-REX experiment (Section 4.5). The intervals where the damping layers are applied are in units of km.</caption><thead><tr><th>Test case</th><th>Δt<br>[s]</th><th>T<sub>f</sub><br>[h]</th><th>Domain<br>[km × km]</th><th>Damping<br>layer (x)</th><th>Damping<br>layer (z)</th><th>λΔt</th></tr></thead><tbody><tr><td>LHMW</td><td>2.5</td><td>15</td><td>240 × 30</td><td>(0,80), (160,240)</td><td>(15,30)</td><td>0.3</td></tr><tr><td>NLNHMW</td><td>1</td><td>5</td><td>40 × 20</td><td>(0,10), (30,40)</td><td>(9,20)</td><td>0.15</td></tr><tr><td>Schär</td><td>5</td><td>10</td><td>100 × 30</td><td>(-50,-20), (20,50)</td><td>(20,30)</td><td>1.2</td></tr><tr><td>NST</td><td>0.5</td><td>6</td><td>100 × 20</td><td>(0,20), (80,100)</td><td>(9,20)</td><td>0.15</td></tr><tr><td>T-REX</td><td>0.75</td><td>4</td><td>100 × 26</td><td>(0,50), (350,400)</td><td>(20,26)</td><td>0.15</td></tr></tbody></table>
with $l$ being a mixing length and $Ri$ denoting the Richardson number given by

$$ Ri = \frac{g}{\theta_0} \left| \frac{\partial \theta}{\partial z} \right|^2. \qquad (15) $$

Here, $\theta_0$ is a reference temperature, while the function $f(Ri)$ is defined as

$$ f(Ri) = (1 + b|Ri|)^{\beta}, \qquad (16) $$

where

$$ \begin{cases} \beta = -2, b = 5 & \text{if } Ri > 0 \\ \beta = \frac{1}{2}, b = 20 & \text{if } Ri < 0. \end{cases} \qquad (17) $$

We also recall the definition of the vertical flux of horizontal momentum (henceforth “momentum flux”) [53]. This is an important diagnostic quantity, often used to verify if a correct orographic response is established in numerical experiments. It is defined as

$$ m(z) = \int_{-\infty}^{\infty} \bar{\rho}(z) u'(x, z) w'(x, z) dx. \qquad (18) $$

Here, $\bar{\rho}$ is the background density, whereas $u'$ and $w'$ denote the deviation from the background state of the horizontal and vertical velocity, respectively. Relation (18) is valid only for an isolated orography profile with a constant background field $\bar{u}$ [53]. Hence, we also compute the full vertical flux of horizontal momentum, given by

$$ m(z) = \int_{-\infty}^{\infty} (\bar{\rho}(z) + \rho'(x, z)) (\bar{u} + u'(x, z)) w'(x, z) dx, \qquad (19) $$

Table 1 reports the parameter values used in the test cases in this work. For time and space discretizations, we consider the IMEX-DG scheme proposed in [45, 46] and successfully employed for atmospheric dynamics in [47, 48]. We consider piecewise polynomials of degree $r = 4$ for the finite element space, unless differently stated. The solver is implemented in the framework of the deal.II library [2, 5], which natively supports high-order polynomial mappings in addition to linear mappings.

## 4.1 Hydrostatic flow over a hill

We first consider the linear hydrostatic configuration employed, e.g., in [28, 47]. The initial state consists of a constant mean flow with $\bar{u} = 20 \text{ m s}^{-1}$ and of an isothermal background profile with temperature $\bar{T} = 250 \text{ K}$. Finally, the initial profile of the Exner pressure is given by

$$ \bar{\Pi} = \exp \left( - \frac{g}{c_p \bar{T}} z \right), \qquad (20) $$

with $c_p = \frac{\gamma}{\gamma-1} R$ denoting the specific heat at constant pressure. The bottom boundary is described by (10), with $h_m = 1 \text{ m}$, $x_c = 120 \text{ km}$, and $a_c = 10 \text{ km}$.

For the construction of the boundary elements, we consider polynomial degrees 2 (i.e., a parabolic profile) and 4, and we compare the results with those obtained employing a linear mapping. The computational mesh is composed by 100 × 60 elements, leading to a resolution of 600 m along the horizontal direction and of 125 m along the
## 4.1 Hydrostatic flow over a hill

We first consider the linear hydrostatic configuration employed e.g. in [23, 42]. The initial state consists of a constant mean flow with $\bar{u} = 20 \text{ m s}^{-1}$ and of an isothermal background profile with temperature $\bar{T} = 250 \text{ K}$. Finally, the initial profile of the Exner pressure is given by

$$ \bar{\Pi} = \exp\left(-\frac{g}{c_p \bar{T}} z\right), \qquad (20) $$

with $c_p = \frac{\gamma}{\gamma-1} R$ denoting the specific heat at constant pressure. The bottom boundary is described by (10), with $h_m = 1 \text{ m}$, $x_c = 120 \text{ km}$, and $a_c = 10 \text{ km}$. We consider for the construction of boundary elements a polynomial degree equal to 2, i.e. a parabolic profile, and equal to 4, and we compare the results with those obtained employing a linear mapping. The computational mesh is composed by 100 × 60 elements, leading to a resolution of 600 m along the horizontal direction and of 125 m along the vertical one. Note that here and in the following, the effective resolution is computed dividing the size of the element along each direction by the polynomial degree, $r = 4$ in this case. From linear theory [47], the analytical momentum flux is given by

$$ m^H = -\frac{\pi}{4} \bar{\rho}_s \bar{u}_s N h_m^2, \qquad (21) $$

where $\bar{\rho}_s$ and $\bar{u}_s$ denote the surface background density and velocity, respectively, and $N$ is the buoyancy frequency. The momentum flux at final time normalized by $m^H$ and computed using curved elements for the bottom boundary provides an improved description of the orographic response compared to the same quantity computed using linear mapping (Figure 2). Moreover, increasing the polynomial degree of the mapping does not lead to further improvement in the accuracy of the results, meaning that the parabolic profile is already an excellent approximation for orographic profile (10). Similar considerations apply to the vertical velocity variable, albeit with slight differences (Figure 3).

<figure><img src="image_3.png" /><figcaption>Figure 2: Linear hydrostatic flow over a hill, computed normalized momentum flux at $t = T_f = 15 \text{ h}$ using a mapping with polynomial degree 4 (black line), a mapping with polynomial degree 2 (blue dots), and a linear mapping (red line).</figcaption></figure>
<table><caption>Table 1: Model parameters for the test cases, see main text for details. LHMW: linear hydrostatic mountain wave (Section 4.1). NLNHMW: nonlinear non-hydrostatic mountain wave (Section 4.2). Schär: Schär profile (Section 4.3). NST: non-smooth orography (Section 4.4). T-REX: Sierra profile, T-REX experiment (Section 4.5). The intervals where the damping layers are applied are in units of km.</caption><thead><tr><th>Test case</th><th>Δt [s]</th><th>T<sub>f</sub> [h]</th><th>Domain [km × km]</th><th>Damping layer (x)</th><th>Damping layer (z)</th><th>λ&#773;Δt</th></tr></thead><tbody><tr><td>LHMW</td><td>2.5</td><td>15</td><td>240 × 30</td><td>(0,80), (160,240)</td><td>(15,30)</td><td>0.3</td></tr><tr><td>NLNHMW</td><td>1</td><td>5</td><td>40 × 20</td><td>(0,10), (30,40)</td><td>(9,20)</td><td>0.15</td></tr><tr><td>Schär</td><td>5</td><td>10</td><td>100 × 30</td><td>(-50,-20), (20,50)</td><td>(20,30)</td><td>1.2</td></tr><tr><td>NST</td><td>0.5</td><td>6</td><td>100 × 20</td><td>(0,20), (80,100)</td><td>(9,20)</td><td>0.15</td></tr><tr><td>T-REX</td><td>0.75</td><td>4</td><td>100 × 26</td><td>(0,50), (350,400)</td><td>(20,26)</td><td>0.15</td></tr></tbody></table>

vertical direction. Note that, here and in the following, the effective resolution is computed dividing the size of the element along each direction by the polynomial degree, $r = 4$ in this case. From linear theory [53], the analytical momentum flux is given by

$$m^H = -\frac{\pi}{4} \bar{\rho}_s \bar{u}_s N h_m^2, \qquad (21)$$

where $\bar{\rho}_s$ and $\bar{u}_s$ denote the surface background density and velocity, respectively, and $N$ is the buoyancy frequency.

The momentum flux at final time normalized by $m^H$ and computed using curved elements for the bottom boundary provides an improved description of the orographic response compared to the same quantity computed using linear mapping (Figure 2). Moreover, increasing the polynomial degree of the mapping does not lead to further improvement in the accuracy of the results, meaning that the parabolic profile is already an excellent approximation for orographic profile (10). Similar considerations apply to the vertical velocity variable, albeit with slight differences (Figure 3). The wall-clock time needed employing high-order mappings is essentially the same as that needed using the linear mapping (Table 2). We refer to Section 4.6 for further considerations on the computational cost of high-order mappings.

<table><caption>Table 2: Linear hydrostatic flow over a hill: wall-clock times (WT) for the linear mapping and the high-order mappings simulations using a computational mesh composed by $N_{el} = 100 \times 60 = 6000$ elements. The overhead using the high-order mapping is computed with respect to the WT of the simulation employing the linear mapping.</caption><thead><tr><th>Configuration</th><th>WT[s]</th><th>Overhead</th></tr></thead><tbody><tr><td>Linear mapping</td><td>10500</td><td>-</td></tr><tr><td>Degree 2 mapping</td><td>10400</td><td>-0.95%</td></tr><tr><td>Degree 4 mapping</td><td>10600</td><td>0.95%</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: Linear hydrostatic flow over a hill, vertical velocity at $t = T_f = 15$ h using the high-order mapping (continuous black lines) and the linear mapping (red dashed lines). Contours in $[-4.0, 4.0] \times 10^{-3} \text{ m s}^{-1}$ with a $5 \times 10^{-4} \text{ m s}^{-1}$ contour interval.</figcaption></figure>

## 4.2 Non-Hydrostatic flow over a hill

Next, we consider the nonlinear non-hydrostatic configuration employed e.g. in [42, 43]. Here, the background velocity is $\bar{u} = 13.28 \text{ m s}^{-1}$, whereas the initial state is described by the following potential temperature $\bar{\theta}$ and Exner pressure $\bar{\Pi}$:

$$ \bar{\theta} = T_{ref} \exp\left(\frac{N^2}{g}z\right) \qquad \bar{\Pi} = 1 + \frac{g^2}{c_p T_{ref} N^2} \left[\exp\left(-\frac{N^2}{g}z\right) - 1\right], \quad (22) $$

with $T_{ref} = 273 \text{ K}$ being the surface temperature and $N = 0.02 \text{ s}^{-1}$ being the constant Brunt-Väisälä frequency. The bottom boundary is described again by (10), with $h_m = 450 \text{ m}$, $x_c = 20 \text{ km}$, and $a_c = 1 \text{ km}$. We consider again elements of polynomial degree equal to 2 and to 4, and we compare the results with those achieved with a linear mapping. The computational mesh is composed by $50 \times 50$ elements, yielding a resolution of $200 \text{ m}$ along the horizontal direction and of $100 \text{ m}$ along the vertical direction.

The use of different polynomial mappings leads to different values of the normalized momentum flux at $t = T_f$ in the vicinity of the orography (Figure 4). The maximum relative difference is around 11.5%, meaning that the geometric error associated with the description of the bottom boundary with a linear mapping is not negligible. Analogous considerations to those reported for the linear hydrostatic test case in Section 4.1 are valid for the vertical velocity (Figure 3).

## 4.3 Flow over Schär profile

Next, we consider the more complex five-peak idealized mountain range (11) originally proposed in [46] (see also [37, 42]). We choose $h_m = 250 \text{ m}$, $a_c = 5 \text{ km}$, and $\lambda_c = 4 \text{ km}$. The initial state has the expression (22), with surface temperature $T_{ref} = 288 \text{ K}$, constant buoyancy frequency $N = 0.01 \text{ s}^{-1}$, and a background horizontal velocity $\bar{u} = 10 \text{ m s}^{-1}$. The domain is $\Omega = (-50, 50) \times (0, 30) \text{ km}$. The mesh is
<figure><img src="image_4.png" /><figcaption>Figure 2: Linear hydrostatic flow over a hill, computed normalized momentum flux at $t = T_f = 15$ h using a mapping with polynomial degree 4 (black line), a mapping with polynomial degree 2 (blue dots), and the linear mapping (red line).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Linear hydrostatic flow over a hill, vertical velocity at $t = T_f = 15$ h using the high-order mapping with polynomial degree 4 (continuous black lines) and the linear mapping (red dashed lines). Contours in the range [-4.0, 4.0] × 10⁻³ m s⁻¹ with a 5 × 10⁻⁴ m s⁻¹ contour interval.</figcaption></figure>

## 4.2 Non-Hydrostatic flow over a hill

Next, we consider the nonlinear non-hydrostatic configuration employed, e.g., in [47, 48]. Here, the background velocity is $\bar{u} = 13.28$ m s⁻¹, whereas the initial state is
<figure><img src="image_5.png" /><figcaption>Figure 4: Nonlinear non-hydrostatic flow over a hill, comparison of normalized momentum flux at $t = T_f = 5$ h using a mapping with polynomial degree 4 (black line), a mapping with polynomial degree 2 (blue dots), and a linear mapping (red line).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Nonlinear non-hydrostatic flow over a hill, vertical velocity at $t = T_f = 5$ h using the high-order mapping (continuous black lines) and the linear mapping (red dashed lines). Contours in [-4.2, 4.0] m s⁻¹ with a 0.586 m s⁻¹ contour interval.</figcaption></figure>

composed by 50 × 25 elements, leading to a resolution of 500 m along the horizontal direction and of 300 m along the vertical one. A reference solution is computed using a mesh composed by 200 × 50 elements with a linear map. We compare the momentum flux (18) at $t = T_f$ in the far-field and we normalize it by the value obtained with the reference solution at $z = 600$ m (Figure 6). It is important to remark that (18) is valid only for an isolated orography with a constant background field $\bar{u}$ [47]. Hence, following [25], we also compute the full vertical flux of horizontal momentum, given by

$$m(z) = \int_{x_{min}}^{x_{max}} (\bar{\rho}(z) + \rho'(x, z)) (\bar{u} + u'(x, z)) w'(x, z) dx, \quad (23)$$

with $x_{min} = 0$ km and $x_{max} = 30$ km.
described by the following potential temperature $\bar{\theta}$ and Exner pressure $\bar{\Pi}$:

$$ \bar{\theta} = T_{ref} \exp\left(\frac{N^2}{g} z\right) \qquad \bar{\Pi} = 1 + \frac{g^2}{c_p T_{ref} N^2} \left[\exp\left(-\frac{N^2}{g} z\right) - 1\right], \quad (22) $$

where $T_{ref} = 273$ K denotes the surface temperature and $N = 0.02 \text{ s}^{-1}$ denotes the constant Brunt-Väisälä frequency. The bottom boundary is described again by (10), with $h_m = 450 \text{ m}$, $x_c = 20 \text{ km}$, and $a_c = 1 \text{ km}$. We consider again elements of polynomial degree equal to 2 and to 4, and we compare the results with those achieved with a linear mapping. The computational mesh is composed by $50 \times 50$ elements, yielding a resolution of $200 \text{ m}$ along the horizontal direction and of $100 \text{ m}$ along the vertical direction.

The use of different polynomial mappings leads to different values of the normalized momentum flux at $t = T_f$ in the lower part of the computational domain (Figure 4). The maximum relative difference is around 8.2%, meaning that the geometric error associated with the description of the bottom boundary with a linear mapping is not negligible. Analogous considerations to those reported for the linear hydrostatic test case in Section 4.1 are valid for the vertical velocity (Figure 3). Finally, for what concerns the computational cost, analogous considerations to those reported in Section 4.1 still hold. We refer again to Section 4.6 for further considerations on the computational cost.

<table><caption>Table 3: Nonlinear non-hydrostatic flow over a hill: wall-clock times (WT) for the linear mapping and the high-order mappings simulations using a computational mesh composed by $N_{el} = 50 \times 50 = 2500$ elements. The overhead using the high-order mapping is computed with respect to the WT of the simulation employing the linear mapping.</caption><thead><tr><th>Configuration</th><th>WT[s]</th><th>Overhead</th></tr></thead><tbody><tr><td>Linear mapping</td><td>3320</td><td>-</td></tr><tr><td>Degree 2 mapping</td><td>3600</td><td>8.4%</td></tr><tr><td>Degree 4 mapping</td><td>3480</td><td>4.8%</td></tr></tbody></table>

## 4.3 Flow over Schär profile

Next, we consider the more complex five-peak idealized mountain range (11) originally proposed in [52] (see also [42, 47]). We choose $h_m = 250 \text{ m}$, $a_c = 5 \text{ km}$, and $\lambda_c = 4 \text{ km}$. The initial state has the expression (22), with surface temperature $T_{ref} = 288 \text{ K}$, constant buoyancy frequency $N = 0.01 \text{ s}^{-1}$, and a background horizontal velocity $\bar{u} = 10 \text{ m s}^{-1}$. The domain is $\Omega = (-50, 50) \times (0, 30) \text{ km}$. The mesh is composed by $50 \times 25$ elements, leading to a resolution of $500 \text{ m}$ along the horizontal direction and of $300 \text{ m}$ along the vertical direction. A reference solution is computed using a mesh composed by $200 \times 50$ elements and a linear mapping to describe the orography. We compute the momentum flux (18) at $t = T_f$ and we normalize it by the value obtained with the reference solution at $z = 600 \text{ m}$ (Figure 6). It is important to remark that (18) is valid only for an isolated orography with a constant background
Significant differences arise in the transfer of the momentum flux along the vertical direction (Figure 6). Moreover, for this test case, increasing the polynomial degree of the mapping yields significantly different results. This is due to the fact that (11) is described by a non-polynomial function and therefore the accuracy in the orography description improves as the polynomial degree of the mapping increases. We also notice that the use of a polynomial degree equal to 4 reduces the oscillations in the profile of the momentum flux. Finally, high-order mappings avoid small oscillations close to the orography for the horizontal velocity deviation (Figure 7) and this further emphasizes the importance of the use of high-order mappings.

<figure><img src="image_7.png" /><figcaption>Figure 6: Flow over Schär profile, normalized momentum flux computed with formula (18) (left) and formula (19) (right) at $t = T_f = 10$ h, using polynomial degree 4 mapping (blue line), polynomial degree 2 mapping (red line), linear mapping (green line). The black line denotes a reference solution with a 200 × 50 mesh using a linear map (black line). In all cases, the momentum flux is normalized using the value at $z = 600$ m obtained with the reference solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Flow over Schär profile, horizontal velocity deviation at $t = T_f = 10$ h using a mapping with polynomial degree 4 (solid black contours) and a linear mapping (dashed red contours). Contours in [-2, 2] m s⁻¹ with a 0.2 m s⁻¹ contour interval.</figcaption></figure>

## 4.4 Non-smooth orography

In this Section, we modify one of the previous classical benchmarks by adding a
<figure><img src="image_6.png" /><figcaption>Figure 4: Nonlinear non-hydrostatic flow over a hill, comparison of normalized momentum flux at $t = T_f = 5$ h using a mapping with polynomial degree 4 (black line), a mapping with polynomial degree 2 (blue dots), and a linear mapping (red line).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Nonlinear non-hydrostatic flow over a hill, vertical velocity at $t = T_f = 5$ h using the high-order mapping with polynomial degree 4 (continuous black lines) and the linear mapping (red dashed lines). Contours in the range [-4.2, 4.0] m s⁻¹ with a 0.586 m s⁻¹ contour interval.</figcaption></figure>

field $\bar{u}$ [53]. Hence, following [30], we also compute the full vertical flux of horizontal momentum, given by

$$m(z) = \int_{x_{min}}^{x_{max}} (\bar{\rho}(z) + \rho'(x, z)) (\bar{u} + u'(x, z)) w'(x, z) dx, \quad (23)$$

with $x_{min} = 0$ km and $x_{max} = 30$ km.
non-smooth perturbation. More specifically, we consider the versiera of Agnesi (10), with $h_m = 0.45$ km, $x_c = 50$ km, and $a_c = 4$ km, and we add the perturbation

$$h'(x) = 1 - 4 \left| x - \left\lfloor x + \frac{1}{2} \right\rfloor \right|, \qquad (24)$$

where the length is in units of kilometers. Hence, the resulting orography reads as follows (Figure 8):

$$h(x) = \begin{cases} \frac{h_m}{1 + \left(\frac{x-x_c}{a_c}\right)^2} + h_m \delta \left(1 - 4 \left| x - \left\lfloor x + \frac{1}{2} \right\rfloor \right| \right) & \text{if } |x - x_c| \le 2a_c \\ \frac{h_m}{1 + \left(\frac{x-x_c}{a_c}\right)^2} & \text{otherwise,} \end{cases} \qquad (25)$$

with $\delta = \frac{1}{40} = 0.025$. The perturbation has zero mean value, so that, if an increasingly strong filter is applied, the original smooth profile is recovered.

<figure><img src="image_9.png" /><figcaption>Figure 8: Non-smooth orography profile obtained with (25) using $\delta = 0.025$.</figcaption></figure>

The computational mesh is composed by 100 × 50 elements, yielding a resolution of 250 m along the horizontal direction and of 100 m along the vertical one. A reference solution is computed using a computational mesh composed by 300 × 50 elements with a linear mapping, so as to allow a correct description of the non-smooth perturbation. The value obtained at $z = 550$ m for the reference solution is adopted to normalize the momentum flux. Following the discussion in [47], the momentum flux (18) is associated with a pressure difference which results in a net drag force on the mountain. For an isolated orography with a constant background field, one can show that (18) is equivalent to

$$m(z) = \int_{-\infty}^{\infty} p'(x, z) \frac{dh(x)}{dx} dx, \qquad (26)$$

with $p'$ denoting the pressure perturbation with respect to the background state. Hence, the momentum flux is strongly influenced by the orography, and therefore it can be employed as an indicator to analyze the quality in the description of the orographic profile.
Significant differences arise in the transfer of the momentum flux along the vertical direction (Figure 6). Moreover, for this test case, increasing the polynomial degree of the mapping yields significantly different results. This is due to the fact that (11) is described by a non-polynomial function. Therefore, the accuracy in the orography description improves as the polynomial degree of the mapping increases. We also notice that the use of a polynomial degree equal to 4 reduces the oscillations in the profile of the momentum flux. Finally, high-order mappings avoid small oscillations close to the orography for the horizontal velocity deviation (Figure 7) and this further emphasizes the importance of the use of high-order mappings. Analogous considerations to those reported in Section 4.1 hold for the computational time (Table 4).

<table><caption>Table 4: Flow over Schär profile: wall-clock times (WT) for the linear mapping simulation and the high-order mappings simulations using a computational mesh composed by $N_{el} = 50 \times 25 = 1250$ elements. The overhead using the high-order mappings is computed with respect to the WT of the simulation employing the linear mapping.</caption><thead><tr><th>Configuration</th><th>WT[s]</th><th>Overhead</th></tr></thead><tbody><tr><td>Linear mapping</td><td>2030</td><td>-</td></tr><tr><td>Degree 2 mapping</td><td>2140</td><td>5.4%</td></tr><tr><td>Degree 4 mapping</td><td>1930</td><td>-4.9%</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 6: Flow over Schär profile, normalized momentum flux computed with formula (18) (left) and formula (19) (right) at $t = T_f = 10$ h, using polynomial degree 4 mapping (blue line), polynomial degree 2 mapping (red line), linear mapping (green line). The black line denotes a reference solution with a $200 \times 50$ mesh using a linear map (black line). In all cases, the momentum flux is normalized using the value at $z = 600$ m obtained with the reference solution.</figcaption></figure>

## 4.4 Non-smooth orography

In this Section, we modify one of the previous classical benchmarks by adding a non-smooth perturbation. More specifically, we consider the versiera of Agnesi (10), with $h_m = 0.45$ km, $x_c = 50$ km, and $a_c = 8$ km, and we add the perturbation

$$h'(x) = 1 - 4 \left| x - \left[ x + \frac{1}{2} \right] \right|, \qquad (24)$$
The use of high-order mappings leads to results in good agreement with the reference ones, while a visible discrepancy arises with the linear mapping (Figure 9). Quantitatively, we compute the $l^2$ relative error in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km, excluding the damping layer (Table 2). The results on the momentum flux computed with (18) show that the use of high-order mappings provides a better description of some small-scale features of the orography even at lower spatial resolution and leads to improved results in the development of lee waves and large-scale features. Analogous considerations hold for the vertical flux of horizontal momentum computed with (19).

<table><caption>Table 2: Flow over non-smooth orography with $\delta = 0.025$, $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19) using linear mapping and degree 3 (high-order) mapping. The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km.</caption><thead><tr><th rowspan="2">t</th><th colspan="2">Linear mapping</th><th colspan="2">High-order mapping</th></tr><tr><th>$l^2(m(z))$ (18)</th><th>$l^2(m(z))$ (19)</th><th>$l^2(m(z))$ (18)</th><th>$l^2(m(z))$ (19)</th></tr></thead><tbody><tr><td>$T_f/2 = 3$ h</td><td>$1.24 \times 10^{-1}$</td><td>$1.26 \times 10^{-1}$</td><td>$1.10 \times 10^{-2}$</td><td>$1.12 \times 10^{-2}$</td></tr><tr><td>$T_f = 6$ h</td><td>$6.08 \times 10^{-1}$</td><td>$6.15 \times 10^{-2}$</td><td>$2.00 \times 10^{-2}$</td><td>$2.00 \times 10^{-2}$</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 9: Flow over non-smooth orography with $\delta = 0.025$, normalized momentum flux computed with (18) at times $t = T_f/2 = 3$ h (left panel) and $t = T_f = 6$ h (right panel). Results with polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 300 × 50 elements mesh and linear mapping (black line). The momentum flux is normalized by the value at the surface obtained with the 300 × 50 elements mesh.</figcaption></figure>

Next, we increase the perturbation and we set $\delta = 0.15$ in (25) (Figure 10). We also employ the simple model for turbulent vertical diffusion and we take $l = 50$ m and $\theta_0 = 273$ K in (14).

A comparison of the contours of the horizontal velocity deviation and of the potential temperature shows that high-order mapping results more closely match the reference results compared to linear mapping results (Figure 11). This is confirmed by the result of the normalized momentum flux at $t = T_f/2$ (Figure 12 and Table 3). In terms of relative $l^2$ error, the high-order mapping gives a result one order of magnitude closer to the reference than the linear mapping. On the other hand, visible differences appear in the normalized momentum flux with the three different configurations at
<figure><img src="image_9.png" /><figcaption>Figure 7: Flow over Schär profile, horizontal velocity deviation at $t = T_f = 10$ h using a mapping with polynomial degree 4 (solid black contours) and a linear mapping (dashed red contours). Contours in the range [-2, 2] m s⁻¹ with a 0.2 m s⁻¹ contour interval.</figcaption></figure>

where $|x|$ denotes the largest integer smaller or equal to $x$ (floor function) and the length is in units of kilometers. Hence, the resulting orography reads as follows (Figure 8):

$$h(x) = \begin{cases} \frac{h_m}{1 + \left(\frac{x-x_c}{a_c}\right)^2} + h_m \delta\left(1 - 4|x - \lfloor x + \frac{1}{2} \rfloor|\right) & \text{if } |x - x_c| \le 2a_c \\ \frac{h_m}{1 + \left(\frac{x-x_c}{a_c}\right)^2} & \text{otherwise,} \end{cases} \quad (25)$$

with $\delta = 0.025$. The perturbation has zero mean value, so that, if an increasingly strong filter is applied, the original smooth profile is recovered. Notice that

$$\frac{Na_c}{\bar{u}} \approx 12 \gg 1, \quad (26)$$

which implies that we are considering a hydrostatic regime.

The computational mesh is composed by 100 × 50 elements, yielding a resolution of 250 m along the horizontal direction and of 100 m along the vertical direction. A reference solution is computed using a computational mesh composed by 300 × 50 elements with a linear mapping, so as to allow a correct description of the non-smooth perturbation. The value obtained at $z = 550$ m for the reference solution is adopted to normalize the momentum flux. Following the discussion in [53], the momentum flux (18) is associated with a pressure difference which results in a net drag force on the mountain. For an isolated orography with a constant background field, one can show that (18) is equivalent to

$$m(z) = \int_{-\infty}^{\infty} p'(x, z) \frac{dh(x)}{dx} dx, \quad (27)$$
<figure><img src="image_11.png" /><figcaption>Figure 10: Non-smooth orography profile obtained with (25) using $\delta = 0.15$.</figcaption></figure>

$t = T_f$ (12). This is likely due to the fact that there is low quantitative predictability of key characteristics for non-smooth orography (see also the discussion in [17]).

While a linear mapping might have seemed the most appropriate choice for the non-smooth, non-differentiable orography profile (25), the use of a linear mapping on coarse meshes yields a poor description of small-scale orographic patterns. It is well known [19, 30] that small-scale orographic patterns can have a significant impact on the development of lee waves and, more in general, on mountain wave-driven atmospheric processes and that the improvement in forecast skill for increasing resolutions is mainly due to the improved representation of the orography.

<table><caption>Table 3: Non-smooth orography with $\delta = 0.15$, $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km.</caption><thead><tr><th rowspan="2">t</th><th colspan="2">Linear mapping</th><th colspan="2">High-order mapping</th></tr><tr><th>l<sup>2</sup>(m(z)) (18)</th><th>l<sup>2</sup>(m(z)) (19)</th><th>l<sup>2</sup>(m(z)) (18)</th><th>l<sup>2</sup>(m(z)) (19)</th></tr></thead><tbody><tr><td>T<sub>f</sub> / 2 = 3 h</td><td>1.63</td><td>1.60</td><td>1.29 &times; 10<sup>-1</sup></td><td>1.32 &times; 10<sup>-1</sup></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>3.06 &times; 10<sup>-1</sup></td><td>2.79 &times; 10<sup>-1</sup></td><td>2.62 &times; 10<sup>-1</sup></td><td>2.63 &times; 10<sup>-1</sup></td></tr></tbody></table>

### 4.4.1 Filtered non-smooth orography

We consider now the results obtained smoothing the orography using filters. This approach seeks to avoid the resolution of small-scale features at a resolution close to that of the mesh, so as to avoid spurious oscillations and problems with physical parametrizations, see e.g. the discussion in [54]. Subgrid-scale orographic (gravity wave) drag parameterizations are employed in NWP and climate models to compensate the insufficient resolution of orographic features [38, 44]. However, the interplay between resolved and parameterized orographic effects is critical and improved results in simulations of mountain atmospheric processes can be obtained increasing
<figure><img src="image_10.png" /><figcaption>Figure 8: Non-smooth orography profile obtained with formula (25) using $\delta = 0.025$.</figcaption></figure>

with $p'$ denoting the pressure perturbation with respect to the background state. Hence, the momentum flux is strongly influenced by the orography, and therefore it can be employed as an indicator to analyze the quality in the description of the orographic profile.

The use of high-order mappings leads to results in good agreement with the reference solution, while a visible discrepancy arises with the linear mapping (Figure 9). Quantitatively, we compute the $l^2$ relative error in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km, excluding the damping layer (Table 5). The results on the momentum flux computed with (18) show that the use of high-order mappings provides a better description of some small-scale features of the orography even at lower spatial resolution and leads to improved results in the development of lee waves and large-scale features. Analogous considerations hold for the vertical flux of horizontal momentum computed with (19). The overhead of the simulation employing the high-order mapping with respect to that using the linear mapping amounts to around 4.8% in terms of wall-clock time (Table 5).

<table><caption>Table 5: Flow over non-smooth orography with $\delta = 0.025$ in (25), $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19) using linear mapping and degree 3 (high-order) mapping. The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km. The overhead using the high-order mapping is computed with respect to the WT of the simulation employing the linear mapping.</caption><thead><tr><th rowspan="2">t</th><th colspan="2">Linear mapping</th><th colspan="3">High-order mapping</th></tr><tr><th>error m(z) (18)</th><th>error m(z) (19)</th><th>error m(z) (18)</th><th>error m(z) (19)</th><th>Overhead</th></tr></thead><tbody><tr><td>T<sub>f</sub> = 3 h</td><td>6.37 × 10<sup>-2</sup></td><td>6.43 × 10<sup>-2</sup></td><td>1.88 × 10<sup>-2</sup></td><td>1.90 × 10<sup>-2</sup></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>3.48 × 10<sup>-2</sup></td><td>3.64 × 10<sup>-2</sup></td><td>6.86 × 10<sup>-3</sup></td><td>7.46 × 10<sup>-3</sup></td><td>4.8%</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 11: Flow over non-smooth orography with $\delta = 0.15$ in (25) at $t = T_f = 6$ h. Top: Reference solution. Middle: Polynomial degree 3 mapping. Bottom: Linear mapping. Horizontal velocity perturbation (colors), contours in the range [-15, 15] m s$^{-1}$ with a 2.5 m s$^{-1}$ interval. Potential temperature (dashed lines), contours in the range [273, 403] K with a 10 K interval.</figcaption></figure>

the resolution, see [19, 30] and the recent contribution of the authors on the use of non-conforming meshes [43].

In order to smooth the orography, we consider a discrete set of $N$ points $h_i$ computed from (25) and we apply different filtering techniques. More specifically, we take $N = 601$ equispaced points and we consider two different approaches, namely a mov-
<figure><img src="image_11.png" /><figcaption>Figure 9: Flow over non-smooth orography with $\delta = 0.025$ in (25), normalized momentum flux computed with (18) at times $t = \frac{T_f}{2} = 3$ h (left panel) and $t = T_f = 6$ h (right panel). Results with polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 300 × 50 elements mesh and linear mapping (black line). The momentum flux is normalized by the value obtained with the 300×50 elements mesh at $z = 550$ m.</figcaption></figure>

Next, we increase the perturbation and we set $\delta = 0.15$ in (25) (Figure 10). We also employ the simple model for turbulent vertical diffusion (14) and we take $l = 50$ m and $\theta_0 = 273$ K in (14).

<figure><img src="image_12.png" /><figcaption>Figure 10: Non-smooth orography profile obtained with formula (25) using $\delta = 0.15$ in (25).</figcaption></figure>

A comparison of the contours of the horizontal velocity deviation and of the potential temperature shows that high-order mapping results more closely match the reference results compared to linear mapping results (Figure 11). This is confirmed by the values of the normalized momentum flux at $t = \frac{T_f}{2}$ and at $t = T_f$ (Figure 12 and Table 6). In terms of relative $l^2$ error with respect to the reference solution, the high-order mapping outperforms the linear mapping by an order of magnitude.
<figure><img src="image_13.png" /><figcaption>Figure 12: Flow over non-smooth orography with $\delta = 0.15$, normalized momentum flux computed with (18) at $t = \frac{T_f}{2}$ (left panel), $t = T_f$ (right panel). Results with polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 300×50 elements mesh and linear mapping (black line). The momentum flux is normalized by the value at $z = 550$ m obtained with the 300×50 elements mesh.</figcaption></figure>

ing average filter and the so-called Raymond filter described in [54]. For the moving average, we define the new discrete orography as

$$ \hat{h}_i = \frac{1}{M+1} \sum_{i-\frac{M}{2}}^{i+\frac{M}{2}} h_i, \qquad (27) $$

with $M + 1$ being the number of values employed for the average. We consider here $M = 4, 6$. The Raymond filter is instead a spectral filter and it is defined by specifying the Fourier coefficients of its response function as

$$ \tilde{F}(\hat{\kappa}) = \frac{1}{1 + \varepsilon \tan^m\left(\frac{1}{2}\hat{\kappa}\Delta\right)}, \qquad (28) $$

with $\hat{\kappa}$ denoting the wave number, $\varepsilon$ being a suitable constant and $m$ being a suitable exponent. The variable $\Delta$ is related to the resolution in the description of the orography. Following [54], we take $m = 6$ and $\varepsilon = 1$, whereas we consider $\Delta = \frac{3}{20}\Delta x, \frac{1}{5}\Delta x, \frac{1}{4}\Delta x$, with $\Delta x = \frac{L}{N-1} = \frac{100}{600} = \frac{1}{6}$. Recall that the units are in kilometers and $L = 100$ km is the length of our computational domain. Since the orography changes after the application of the filter, a reference solution using the linear mapping on the 300 × 50 elements mesh is computed for each filtered configuration. A cubic spline interpolation [14] is employed to evaluate the orography at the points of this finer mesh.

First, we analyze the results obtained with the moving-average filter (27) (Figure 13). For $M = 4$ and in terms of $l^2$ relative error, a comparison of the normalized momentum flux (18) shows that the high-order mapping performs an order of magnitude better than the linear mapping at $t = \frac{T_f}{2}$ and almost three times better at $t = \frac{T_f}{2}$ (Figure 14 and Table 4). Analogous results are obtained for $M = 6$, for which the better accuracy established by the high-order is even more evident. This is due to the fact that few small-scale topographic features are preserved by the filter, and such features, at a coarse resolution, can be properly captured only by means of a high-order mapping.

Next, we consider the spectral filter (28). We immediately notice that the filtered orography strongly depends on the parameter $\Delta$ (Figure 15, 16). For $\Delta = 0.2\Delta x$, the
Slightly larger values of the error of the momentum flux using the high-order mapping are obtained for larger values of the final time, due to the severity of the artificial orography profile (not shown). Moreover, the overhead of the simulation employing the high-order mapping with respect to that using the linear mapping amounts to ≈ 71% in terms of wall-clock time (Table 6). The possible reasons for this overhead will be discussed in Section 4.6.

While a linear mapping might seem the most appropriate choice for the non-smooth, non-differentiable orography profile (25), the use of a linear mapping on coarse meshes yields a poor description of small-scale orographic patterns. It is well known [22, 35] that small-scale orographic patterns can have a significant impact on the development of lee waves and, more in general, on mountain wave-driven atmospheric processes. Improved representation of the orography significantly contributes to improvements in forecast skill at increasing resolutions.

<table><caption>Table 6: Flow over non-smooth orography with $\delta = 0.15$ in (25), $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km. The overhead is computed with respect to the WT of the simulation employing the linear mapping.</caption><thead><tr><th rowspan="2">t</th><th colspan="2">Linear mapping</th><th colspan="3">High-order mapping</th></tr><tr><th>error m(z) (18)</th><th>error m(z) (19)</th><th>error m(z) (18)</th><th>error m(z) (19)</th><th>Overhead</th></tr></thead><tbody><tr><td>T<sub>f</sub> = 3 h</td><td>1.30</td><td>1.31</td><td>1.72 × 10<sup>-1</sup></td><td>1.69 × 10<sup>-1</sup></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>3.79 × 10<sup>-1</sup></td><td>3.76 × 10<sup>-1</sup></td><td>3.26 × 10<sup>-2</sup></td><td>3.19 × 10<sup>-2</sup></td><td>71%</td></tr></tbody></table>

### 4.4.1 Filtered non-smooth orography

Next, we consider the results obtained smoothing the orography using filters. This approach seeks to avoid the resolution of small-scale features at a resolution close to that of the mesh, so as to avoid spurious oscillations and problems with physical parametrizations, see, e.g., the discussion in [60]. Subgrid-scale orographic (gravity wave) drag parametrizations are employed in NWP and climate models to compensate the insufficient resolution of orographic features [43, 49]. However, the interplay between resolved and parameterized orographic effects is critical, and improved results in simulations of mountain atmospheric processes can be obtained by increasing the resolution, see [22, 35], and the recent contribution of the authors on the use of non-conforming meshes [48].

In order to smooth the orography, we consider a discrete set of $N$ values $h_i$ computed from (25) and we apply two different filtering techniques: the moving average filter and the Raymond filter described in [60]. We take $N = 601$ equispaced points. For the moving average, we define the new discrete orography as

$$ \hat{h}_i = \frac{1}{M+1} \sum_{i-\frac{M}{2}}^{i+\frac{M}{2}} h_i, \qquad (28) $$
<table><caption>Table 4: Flow over non-smooth orography with $\delta = 0.15$ filtered using the moving-average filter (27), $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km.</caption><thead><tr><th rowspan="2">M</th><th rowspan="2">t</th><th colspan="2">Linear mapping</th><th colspan="2">High-order mapping</th></tr><tr><th>l<sup>2</sup>(m(z)) (18)</th><th>l<sup>2</sup>(m(z)) (19)</th><th>l<sup>2</sup>(m(z)) (18)</th><th>l<sup>2</sup>(m(z)) (19)</th></tr></thead><tbody><tr><td rowspan="2">4</td><td>T<sub>f</sub> = 3 h</td><td>1.41 &times; 10<sup>-1</sup></td><td>1.43 &times; 10<sup>-1</sup></td><td>1.05 &times; 10<sup>-2</sup></td><td>1.06 &times; 10<sup>-2</sup></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>7.42 &times; 10<sup>-2</sup></td><td>7.51 &times; 10<sup>-2</sup></td><td>2.65 &times; 10<sup>-2</sup></td><td>2.63 &times; 10<sup>-2</sup></td></tr><tr><td rowspan="2">6</td><td>T<sub>f</sub> = 3 h</td><td>5.69 &times; 10<sup>-2</sup></td><td>5.62 &times; 10<sup>-2</sup></td><td>2.19 &times; 10<sup>-3</sup></td><td>2.29 &times; 10<sup>-3</sup></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>1.01 &times; 10<sup>-1</sup></td><td>1.01 &times; 10<sup>-1</sup></td><td>4.98 &times; 10<sup>-3</sup></td><td>4.79 &times; 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Figure 13: Non-smooth orography with $\delta = 0.15$, smoothed profiles with moving-average filter (27). Top: $M = 4$. Bottom: $M = 6$.</figcaption></figure>

high-order mapping provides a higher accuracy than the linear mapping. On the other hand, for $\Delta = 0.15\Delta x$ and $\Delta = 0.25\Delta x$, the linear mapping performs globally better than the high-order mapping. This is likely due to the fact that the filtered orography for these values of $\Delta$ is characterized by spurious high frequency features, which are not resolved by the linear mapping and induce instead a spurious response of the high order mapping. The obtained results display a very strong sensitivity with respect to
<figure><img src="image_13.png" /><figcaption>Figure 11: Flow over non-smooth orography with $\delta = 0.15$ in (25) at $t = T_f = 6$ h. Top: Reference solution. Middle: Polynomial degree 3 mapping. Bottom: Linear mapping. Horizontal velocity perturbation (colors), contours in the range $[-15, 15]$ m s$^{-1}$ with a 2.5 m s$^{-1}$ interval. Potential temperature (dashed lines), contours in the range [273, 403] K with a 10 K interval.</figcaption></figure>

with $M + 1$ being the number of values employed for the average. We consider here $M = 4, 6$. The Raymond filter is instead a spectral filter and it is defined by specifying the Fourier coefficients of its response function as

$$ \tilde{F}(\hat{\kappa}) = \frac{1}{1 + \varepsilon \tan^m\left(\frac{1}{2}\hat{\kappa}\Delta\right)}, \qquad (29) $$
<figure><img src="image_15.png" /><figcaption>Figure 14: Flow over non-smooth orography with $\delta = 0.15$ filtered using the moving-average filter (27) and $M = 4$ (top row), $M = 4$ (bottom row), normalized momentum flux computed with (18) at $t = \frac{T_f}{2}$ (left column) and $t = T_f$ (right column). Results with polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 300 × 50 elements mesh and a linear mapping (black line). The momentum flux is normalized by the value at the $z = 550$ m obtained with the 300 × 50 elements mesh.</figcaption></figure>

the parameter $\Delta$ and highlight the difficulty of applying the filtering approach (28) for mesoscale applications.

<table><caption>Table 5: Flow over non-smooth orography with $\delta = 0.15$ filtered using the spectral filter (28), $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km.</caption><thead><tr><th rowspan="2">Δ</th><th rowspan="2">t</th><th colspan="2">Linear mapping</th><th colspan="2">High-order mapping</th></tr><tr><th>l²(m(z)) (18)</th><th>l²(m(z)) (19)</th><th>l²(m(z)) (18)</th><th>l²(m(z)) (19)</th></tr></thead><tbody><tr><td rowspan="2">0.15Δx</td><td>T<sub>f</sub> / 2 = 3 h</td><td>3.05 × 10<sup>-2</sup></td><td>3.17 × 10<sup>-2</sup></td><td>4.22 × 10<sup>-2</sup></td><td>4.32 × 10<sup>-2</sup></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>2.46 × 10<sup>-2</sup></td><td>2.38 × 10<sup>-2</sup></td><td>2.61 × 10<sup>-2</sup></td><td>2.60 × 10<sup>-2</sup></td></tr><tr><td rowspan="2">0.2Δx</td><td>T<sub>f</sub> / 2 = 3 h</td><td>3.38 × 10<sup>-2</sup></td><td>3.52 × 10<sup>-2</sup></td><td>1.18 × 10<sup>-2</sup></td><td>1.23 × 10<sup>-2</sup></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>2.43 × 10<sup>-2</sup></td><td>2.33 × 10<sup>-3</sup></td><td>8.61 × 10<sup>-3</sup></td><td>8.49 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">0.25Δx</td><td>T<sub>f</sub> / 2 = 3 h</td><td>6.18 × 10<sup>-2</sup></td><td>6.33 × 10<sup>-2</sup></td><td>4.48 × 10<sup>-2</sup></td><td>4.54 × 10<sup>-2</sup></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>9.70 × 10<sup>-3</sup></td><td>9.18 × 10<sup>-3</sup></td><td>4.97 × 10<sup>-2</sup></td><td>4.99 × 10<sup>-2</sup></td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>Figure 12: Flow over non-smooth orography with $\delta = 0.15$ in (25), normalized momentum flux computed with (18) at $t = \frac{T_f}{2}$ (left panel), $t = T_f$ (right panel). Results with polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 300 × 50 elements mesh and linear mapping (black line). The momentum flux is normalized by the value at $z = 550$ m obtained with the 300 × 50 elements mesh.</figcaption></figure>

where $\hat{\kappa}$ denotes the wave number, $\varepsilon$ a suitable constant, and $m$ a suitable exponent. The variable $\Delta$ is related to the resolution in the description of the orography. Since the orography changes after the application of the filter, a reference solution using the linear mapping on the 300 × 50 elements mesh is computed for each filtered configuration. A cubic spline interpolation [17] is employed to evaluate the orography at the points of this finer mesh.

First, we analyze the results obtained with the moving-average filter (28) (Figure 13). For $M = 4$ and in terms of $l^2$ relative error, a comparison of the normalized momentum flux (18) shows that the high-order mapping outperforms the linear mapping by a factor 3 at $t = \frac{T_f}{2}$ and by a factor 4 at $t = T_f$ (Figure 14 and Table 7). Analogous results are obtained for $M = 6$, for which the better accuracy established by the high-order mapping is even more evident. This is due to the fact that only few small-scale topographic features are preserved by the filter, and such features, at a coarse resolution, can be properly captured only by means of a high-order mapping. The overhead of the simulations employing the high-order mapping with respect to those using the linear mapping amounts around 15 – 20% in terms of wall-clock time (Table 7).

Next, we consider the spectral filter (29). Following [60], we take $m = 6$ and $\varepsilon = 1$, whereas we consider

$$ \Delta = \frac{3}{20}\Delta x, \frac{1}{5}\Delta x, \frac{1}{4}\Delta x, $$

with

$$ \Delta x = \frac{L}{N-1} = \frac{100}{600} = \frac{1}{6}. $$

We immediately notice that the filtered orography strongly depends on the parameter $\Delta$ (Figure 15, 16). For $\Delta = 0.2\Delta x$, the high-order mapping provides a higher accuracy than the linear mapping. On the other hand, for $\Delta = 0.15\Delta x$ and $\Delta = 0.25\Delta x$, the linear mapping performs globally better than the high-order mapping. This is likely due to the fact that the filtered orography for these values of $\Delta$ is characterized
<figure><img src="image_16.png" /><figcaption>Figure 15: As in Figure 13, but using the spectral filter (28) as defined in [54]. Top: $\Delta = 0.15\Delta x$. Middle: $\Delta = 0.2\Delta x$. Bottom: $\Delta = 0.25\Delta x$.</figcaption></figure>

## 4.5 T-REX Mountain-Wave

In a final test, we consider simulations of a flow over a realistic complex orography [17, 43], whose profile is shown in Figure 17. The initial state, reported in Figure 18, is horizontally homogeneous and it is based on conditions observed during the Intensive Observation Period (IOP) 6 of the Terrain-Induced Rotor Experiment (T-REX) [17].
<table><caption>Table 7: Flow over non-smooth orography with $\delta = 0.15$ in (25) filtered using the moving-average filter (28), $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km. Here, and in the following tables and figures, $M + 1$ is the number of values employed for the moving-average filter. The overhead using the high-order mapping is computed with respect to the WT of the simulation employing the linear mapping.</caption><thead><tr><th rowspan="2">M</th><th rowspan="2">t</th><th colspan="4">Linear mapping</th><th colspan="4">High-order mapping</th></tr><tr><th>error<br>(18)</th><th>m(z)</th><th>error<br>(19)</th><th>m(z)</th><th>error<br>(18)</th><th>m(z)</th><th>error<br>(19)</th><th>m(z)</th><th>Overhead</th></tr></thead><tbody><tr><td rowspan="2">4</td><td>T<sub>f</sub>/2 = 3 h</td><td>7.97 &times; 10<sup>-2</sup></td><td></td><td>8.04 &times; 10<sup>-2</sup></td><td></td><td>2.52 &times; 10<sup>-2</sup></td><td></td><td>2.55 &times; 10<sup>-2</sup></td><td></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>4.16 &times; 10<sup>-2</sup></td><td></td><td>4.05 &times; 10<sup>-2</sup></td><td></td><td>7.39 &times; 10<sup>-3</sup></td><td></td><td>9.97 &times; 10<sup>-3</sup></td><td></td><td>12%</td></tr><tr><td rowspan="2">6</td><td>T<sub>f</sub>/2 = 3 h</td><td>5.39 &times; 10<sup>-2</sup></td><td></td><td>5.42 &times; 10<sup>-2</sup></td><td></td><td>6.00 &times; 10<sup>-3</sup></td><td></td><td>5.97 &times; 10<sup>-3</sup></td><td></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>2.77 &times; 10<sup>-2</sup></td><td></td><td>2.69 &times; 10<sup>-2</sup></td><td></td><td>2.69 &times; 10<sup>-3</sup></td><td></td><td>4.23 &times; 10<sup>-3</sup></td><td></td><td>20%</td></tr></tbody></table>

by spurious high frequency features, which are not resolved by the linear mapping and induce instead a spurious response of the high order mapping. The obtained results display a very strong sensitivity with respect to the parameter $\Delta$ and highlight the difficulty of applying the filtering approach (29) for mesoscale applications. The overhead of the simulations employing the high-order mapping with respect to those using the linear mapping amounts to around 10% in terms of wall-clock time (Table 8).

<table><caption>Table 8: Flow over non-smooth orography with $\delta = 0.15$ in (25) filtered using the spectral filter (29), $l^2$ relative errors on the normalized momentum flux computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 1$ km and $z_2 = 9$ km. The overhead is computed with respect to the WT time of the simulation employing the linear mapping.</caption><thead><tr><th rowspan="2">Δ</th><th rowspan="2">t</th><th colspan="4">Linear mapping</th><th colspan="4">High-order mapping</th></tr><tr><th>error<br>(18)</th><th>m(z)</th><th>error<br>(19)</th><th>m(z)</th><th>error<br>(18)</th><th>m(z)</th><th>error<br>(19)</th><th>m(z)</th><th>Overhead</th></tr></thead><tbody><tr><td rowspan="2">0.15Δx</td><td>T<sub>f</sub>/2 = 3 h</td><td>5.80 &times; 10<sup>-3</sup></td><td></td><td>5.83 &times; 10<sup>-3</sup></td><td></td><td>2.14 &times; 10<sup>-2</sup></td><td></td><td>2.15 &times; 10<sup>-2</sup></td><td></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>6.41 &times; 10<sup>-3</sup></td><td></td><td>6.61 &times; 10<sup>-3</sup></td><td></td><td>2.22 &times; 10<sup>-2</sup></td><td></td><td>2.20 &times; 10<sup>-2</sup></td><td></td><td>7.1%</td></tr><tr><td rowspan="2">0.2Δx</td><td>T<sub>f</sub>/2 = 3 h</td><td>6.41 &times; 10<sup>-3</sup></td><td></td><td>6.58 &times; 10<sup>-3</sup></td><td></td><td>4.52 &times; 10<sup>-3</sup></td><td></td><td>4.49 &times; 10<sup>-3</sup></td><td></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>1.02 &times; 10<sup>-2</sup></td><td></td><td>1.01 &times; 10<sup>-2</sup></td><td></td><td>1.16 &times; 10<sup>-2</sup></td><td></td><td>1.18 &times; 10<sup>-2</sup></td><td></td><td>9.6%</td></tr><tr><td rowspan="2">0.25Δx</td><td>T<sub>f</sub>/2 = 3 h</td><td>1.13 &times; 10<sup>-2</sup></td><td></td><td>1.13 &times; 10<sup>-2</sup></td><td></td><td>3.61 &times; 10<sup>-2</sup></td><td></td><td>3.63 &times; 10<sup>-2</sup></td><td></td><td></td></tr><tr><td>T<sub>f</sub> = 6 h</td><td>6.65 &times; 10<sup>-3</sup></td><td></td><td>7.02 &times; 10<sup>-3</sup></td><td></td><td>2.78 &times; 10<sup>-2</sup></td><td></td><td>2.79 &times; 10<sup>-2</sup></td><td></td><td>13.3%</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>Figure 16: As in Figure 14, but using the spectral filter (28) with $\Delta = 0.15\Delta x$ (top row), $\Delta = 0.2\Delta x$ (middle row), $\Delta = 0.25\Delta x$ (bottom row).</figcaption></figure>

The pressure is computed from the hydrostatic balance, namely

$$p(z) = p_0 \exp \left( -\frac{g}{R} \int_{0}^{z} \frac{1}{T(s)} ds \right), \qquad (29)$$

with $p_0 = 10^5$ Pa. Linear interpolation is used to evaluate both temperature and horizontal velocity outside the given data points. We consider a DG spatial discretization using degree $r = 2$ polynomials. The computational mesh is composed by 150 × 52 elements, yielding a resolution of around 133.33 m along the horizontal direction and of 250 m along the vertical one. A reference solution is computed using a computational mesh composed by 400 × 52 elements. Notice that, as reported in [17], the orographic profile has already been filtered to remove high frequency variations. Therefore, we do not employ the filtering approaches described in Section 4.4. We take $l = 100$ m and $\theta_0 = 273$ K in (14). The vertical turbulent diffusion model is necessary to obtain a stable numerical solution.

A comparison of the contours of the horizontal velocity deviation and of the potential temperature shows that high-order mapping results are more similar to the reference ones with respect to those obtained with the linear mapping (Figure 19).
<figure><img src="image_15.png" /><figcaption>Figure 13: Non-smooth orography with $\delta = 0.15$ in formula (25), smoothed profiles with moving-average filter (28). Top: $M = 4$. Bottom: $M = 6$.</figcaption></figure>

## 4.5 T-REX Mountain-Wave

In a final test, we consider simulations of a flow over a real orography profile (Figure 17, see also [20, 48]). The initial state, reported in Figure 18, is horizontally homogeneous and it is based on conditions observed during the Intensive Observation Period 6 of the Terrain-Induced Rotor Experiment (T-REX) [20]. The pressure is computed from the temperature using hydrostatic balance, namely

$$p(z) = p_0 \exp \left( -\frac{g}{R} \int_0^z \frac{1}{T(s)} ds \right), \quad (30)$$

with $p_0 = 10^5$ Pa. Linear interpolation is used to evaluate both temperature and horizontal velocity outside the given data points. We consider a DG spatial discretization using degree $r = 2$ polynomials. The computational mesh is composed by 150 × 52 elements, yielding a resolution of around 133.33 m along the horizontal direction and of 250 m along the vertical direction. A reference solution is computed using a computational mesh composed by 400 × 52 elements. Note that, as reported in [20], the
<figure><img src="image_18.png" /><figcaption>Figure 17: T-REX mountain-wave test, Sierra profile.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: T-REX mountain-wave test case, initial conditions. Left: horizontal velocity. Right: potential temperature.</figcaption></figure>

Finally, from a more quantitative point of view, we compute the momentum flux. Following [17], we use the average value of the velocity along the horizontal and vertical direction to compute $u'$ and $w'$ in (18). As discussed in [17], there is low predictability of key characteristics such as the strength of downslope winds or the location and intensity of stratospheric wave breaking. Moreover, the change in the resolution of the orography has been shown to modify the representation of mountain wave-driven middle atmosphere processes [30]. In spite of this, the values of the momentum flux obtained with the high-order mapping are much closer to the reference ones than those obtained with the linear mapping, see Figure 20 and Table 6. In terms of $l^2$ relative error, the high-order mapping is almost three times smaller than the error using the linear mapping. These results corroborate those obtained in the non-smooth orography test case, and confirm the advantage of high-order mapping over linear mapping in resolving small-scale orographic features.

# 5 Conclusions

We have performed a quantitative study on the impact of curved elements for flows over orography. To the best of our knowledge, this study has been never performed for
<figure><img src="image_16.png" /><figcaption>Figure 14: Flow over non-smooth orography with $\delta = 0.15$ in (25) filtered using the moving-average filter (28) and $M = 4$ (top row), $M = 6$ (bottom row), normalized momentum flux computed with (18) at $t = \frac{T_f}{2}$ (left column) and $t = T_f$ (right column). Results with polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 300 × 50 elements mesh and a linear mapping (black line). The momentum flux is normalized by the value at the $z = 550$ m obtained with the 300 × 50 elements mesh.</figcaption></figure>

orographic profile has already been filtered to remove high frequency variations. Therefore, we do not employ the filtering approaches described in Section 4.4. The vertical turbulent diffusion model (14) is again necessary to obtain a stable numerical solution and it is used with $l = 100$ m and $\theta_0 = 273$ K.

A comparison of the contours of the horizontal velocity deviation and of the potential temperature shows that high-order mapping results are more similar to the reference results with respect to those obtained with the linear mapping (Figure 19).

Finally, we provide a more quantitative comparison by computing the vertical momentum flux. Following [20], we use the average value of the velocity along the horizontal and vertical direction to compute $u'$ and $w'$ in (18). As discussed in [20], there is low predictability of key characteristics such as the strength of downslope winds or the location and intensity of stratospheric wave breaking. Moreover, the change in the resolution of the orography has been shown to modify the representation of mountain wave-driven middle atmosphere processes [35]. Nevertheless, the values of the momentum flux obtained with the high-order mapping are much closer to the reference results than those obtained with the linear mapping, see Figure 20 and Table 9. In terms of $l^2$ relative error, the high-order mapping is almost three times smaller than the error using the linear mapping. These results corroborate those obtained in the
<figure><img src="image_20.png" /><figcaption>Figure 19: T-REX mountain wave test case at $t = T_f = 4$ h. Top: Reference solution. Middle: Polynomial degree 3 mapping. Bottom: Linear mapping. Horizontal velocity perturbation (colors), contours in the range [−50, 50] m s⁻¹ with a 5 m s⁻¹ interval. Potential temperature (dashed lines), contours in the range [270, 800] K with a 10 K interval.</figcaption></figure>

atmospheric applications. We have employed the IMEX-DG solver originally proposed in [41] and validated for atmospheric applications in [42, 43], which has been implemented in the framework of the open-source library *deal.II* [2]. The software natively supports high-order polynomials for the mapping between the reference element and elements in the physical space.
<figure><img src="image_17.png" /><figcaption>Figure 15: As in Figure 13, but using the spectral filter (29) as defined in [60]. Top: $\Delta = 0.15\Delta x$. Middle: $\Delta = 0.2\Delta x$. Bottom: $\Delta = 0.25\Delta x$.</figcaption></figure>
<table><caption>Table 6: T-REX mountain wave test case, $l^2$ relative errors on the normalized momentum flux at $t = T_f = 4$ h computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 5$ km and $z_2 = 20$ km.</caption><thead><tr><th colspan="2">Linear mapping</th><th colspan="2">High-order mapping</th></tr><tr><th>l<sup>2</sup>(m(z)) (18)</th><th>l<sup>2</sup>(m(z)) (19)</th><th>l<sup>2</sup>(m(z)) (18)</th><th>l<sup>2</sup>(m(z)) (19)</th></tr></thead><tbody><tr><td>1.66</td><td>1.78</td><td>5.01 &times; 10<sup>-1</sup></td><td>6.10 &times; 10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>Figure 20: T-REX mountain wave test case, comparison of normalized momentum flux at $t = T_f = 4$ h. Left: momentum flux computed with (18) using the average values of the velocity to compute $u'$ and $w'$. Right: momentum flux computed with (19). Results using polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 400 × 52 elements mesh and a linear mapping (black line). The momentum flux is normalized by the value at $z = 5$ km obtained with the 400 × 52 elements mesh.</figcaption></figure>

Numerical experiments on a number of classical benchmarks show that, at a given resolution, results obtained with a high-order mapping significantly outperform those obtained using a linear mapping in terms of error with respect to reference solutions. In fact, the use of high-order mapping leads to results which are analogous of those obtained at very high resolution and using a linear mapping, which well resolve the orographic profile. Hence, employing high-order mapping is to some extent equivalent to consider a sub-tessellation to evaluate the integrals. Taking into account the improved representation of the geometry, such results should be considered as the reference ones for future works. Applications of the use of high-order mappings to non-smooth orography profiles have also been presented, showing that the use of high-order mappings provides in general better results with respect to those obtained with a linear mapping. This is valid also for orography profiles obtained with standard filtering procedures. At a given spatial resolution, high-order mappings capture small-scale topographic features inaccessible to the standard linear mapping. Performing very high resolution simulations is computationally expensive, but it is beneficial for mountain atmospheric processes and forecast skill. The use of high-order mappings can be considered as a valuable alternative tool to resolve orographic features at lower spatial resolutions and at a feasible computational cost.
<figure><img src="image_18.png" /><figcaption>Figure 16: As in Figure 14, but using the spectral filter (29) with $\Delta = 0.15\Delta x$ (top row), $\Delta = 0.2\Delta x$ (middle row), $\Delta = 0.25\Delta x$ (bottom row).</figcaption></figure>

non-smooth orography test case, and confirm the advantage of high-order mapping over linear mapping in resolving small-scale orographic features. The overhead of the simulations employing the high-order mapping with respect to those using the linear mapping amounts to around the 37% in terms of wall-clock time (Table 9).

## 4.6 Considerations on the computational cost

In this Section, we report some considerations on the computational cost. The use of high-order mappings has an impact on the quadrature rule (9). More specifically, the determinant of the Jacobian is not a constant, but a polynomial. In this work, following, e.g., [26], we consider exact integration. In order to achieve this goal, we employ the so-called over-integration or consistent integration. This means that we use more than $r + 1$ Gauss quadrature points along each coordinate direction and,
# Acknowledgements

The simulations have been partly run at CINECA thanks to the computational resources made available through the ISCRA-C projects FEMTUF - HP10CTQ8X7 and FEM-GPU - HP10CQYKJ1. We acknowledge the CINECA award, for the availability of high-performance computing resources and support. This work has been partly supported by the ESCAPE-2 project, European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 800897). The authors also gratefully acknowledge Dr. Doyle and Dr. Gaberšek for having provided the orography profile and the initial data of the T-REX mountain wave.

# References

[1] P.F. Antonietti et al. *lymph: discontinuous poLYtopal methods for Multi-PHysics differential problems*. 2024. arXiv: 2401.13376 [math.NA].

[2] D. Arndt et al. “The deal.II library, Version 9.5”. *Journal of Numerical Mathematics* 31.3 (2023), pp. 231–246.

[3] D.N. Arnold, D. Boffi, and F. Bonizzoni. “Finite element differential forms on curvilinear cubic meshes and their approximation properties”. *Numerische Mathematik* 129.1 (2015), pp. 1–20.

[4] D.N. Arnold, R. Falk, and R. Winther. “Finite element exterior calculus: from Hodge theory to numerical stability”. *Bulletin of the American mathematical society* 47.2 (2010), pp. 281–354.

[5] W. Bangerth, R. Hartmann, and G. Kanschat. “deal II: a general-purpose object-oriented finite element library”. *ACM Transactions on Mathematical Software* 33 (2007), pp. 24–51.

[6] F. Bassi and S. Rebay. “A High-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations”. *Journal of Computational Physics* 131 (1997), pp. 267–279.

[7] L. Bonaventura. “A Semi-Implicit, Semi-Lagrangian Scheme Using the Height Coordinate for a Nonhydrostatic and Fully Elastic Model of Atmospheric Flows.” *Journal of Computational Physics* 158 (2000), pp. 186–213.

[8] L. Bonaventura and R. Ferretti. “Semi-Lagrangian methods for parabolic problems in divergence form.” *SIAM Journal of Scientific Computing* 36 (2014), A2458–A2477.

[9] M. Ciallella et al. “Shifted boundary polynomial corrections for compressible flows: high order on curved domains using linear meshes”. *Applied Mathematics and Computation* 441 (2023), p. 127698.

[10] M. Ciallella et al. *Very high order treatment of embedded curved boundaries in compressible flows: ADER discontinuous Galerkin with a space-time Reconstruction for Off-site data*. 2023. arXiv: 2312.07170 [math.NA].
<figure><img src="image_19.png" /><figcaption>Figure 17: T-REX mountain-wave test, Sierra height profile.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18: T-REX mountain-wave test case, initial conditions. Left: horizontal velocity. Right: potential temperature.</figcaption></figure>

<table><caption>Table 9: T-REX mountain wave test case, $l^2$ relative errors on the normalized momentum flux at $t = T_f = 4$ h computed with both (18) and (19). The relative error is computed with respect to the reference solution in the region $[z_1, z_2]$, with $z_1 = 5$ km and $z_2 = 20$ km. The overhead using the high-order mapping is computed with respect to the WT time of the simulation employing the linear mapping.</caption><thead><tr><th colspan="2">Linear mapping</th><th colspan="3">High-order mapping</th></tr><tr><th>error m(z) (18)</th><th>error m(z) (19)</th><th>error m(z) (18)</th><th>error m(z) (19)</th><th>Overhead</th></tr></thead><tbody><tr><td>1.64</td><td>1.84</td><td>4.42 × 10<sup>-1</sup></td><td>5.85 × 10<sup>-1</sup></td><td>37%</td></tr></tbody></table>

in particular, we use $2r + 1$ quadrature points in the case of linear mapping and $2r + 1 + [\hat{r} - 2]$ quadrature points in the case of high-order mapping, where $\hat{r}$ is the polynomial degree of the high-order mapping. We also recall that $r$ denotes the polynomial degree employed for the spatial discretization. One can easily notice that, for $\hat{r} = 2$, the same number of quadrature points employed for linear mapping are sufficient to obtain exact integration. The overhead in the use of high-order mapping is therefore due to the use of higher order quadrature formulas for $\hat{r} > 2$.

However, one can notice from the results in Section 4.1-4.5 that the overhead typ-
[11] P.G. Ciarlet. *The finite element method for elliptic problems*. SIAM, 2002.

[12] R. Costa et al. "Very high-order accurate finite volume scheme on curved boundaries for the two-dimensional steady-state convection-diffusion equation with Dirichlet condition". *Applied Mathematical Modelling* 54 (2018), pp. 752-767.

[13] R. Costa et al. "Very high-order accurate finite volume scheme for the convection-diffusion equation with general boundary conditions on arbitrary curved boundaries". *International Journal for Numerical Methods in Engineering* 117.2 (2019), pp. 188-220.

[14] C. De Boor. *A practical guide to splines*. Vol. 27. springer-verlag New York, 1978.

[15] S. Dey, R.M. O'Bara, and M.S. Shephard. "Towards curvilinear meshing in 3D: the case of quadratic simplices". *Computer-Aided Design* 33.3 (2001), pp. 199-209.

[16] V. Dolejší. "On the discontinuous Galerkin method for the numerical solution of the Navier-Stokes equations". *International Journal for Numerical Methods in Fluids* 45 (2004), pp. 1083-1106.

[17] J. D. Doyle et al. "An intercomparison of T-REX mountain-wave simulations and implications for mesoscale predictability". *Monthly Weather Review* 139.9 (2011), pp. 2811-2831.

[18] M. Fortunato and P.-O. Persson. "High-order unstructured curved mesh generation using the Winslow equations". *Journal of Computational Physics* 307 (2016), pp. 1-14.

[19] D.C. Fritts et al. "Impacts of limited model resolution on the representation of mountain wave and secondary gravity wave dynamics in local and global models. 1: Mountain waves in the stratosphere and mesosphere". *Journal of Geophysical Research: Atmospheres* 127.9 (2022), e2021JD035990.

[20] E. Gaburro et al. "High order entropy preserving ADER-DG schemes". *Applied Mathematics and Computation* 440 (2023), p. 127644.

[21] T. Gal-Chen and R.C.J. Somerville. "On the use of a coordinate transformation for the solution of the Navier-Stokes equations". *Journal of Computational Physics* 17 (1975), pp. 209-228.

[22] F.X. Giraldo. *An Introduction to Element-Based Galerkin Methods on Tensor-Product Bases*. Springer Nature, 2020.

[23] F.X. Giraldo and M. Restelli. "A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases". *Journal of Computational Physics* 227 (2008), pp. 3849-3877.
<figure><img src="image_21.png" /><figcaption>Figure 19: T-REX mountain wave test case at $t = T_f = 4$ h. Top: Reference solution. Middle: Polynomial degree 3 mapping. Bottom: Linear mapping. Horizontal velocity perturbation (colors), contours in the range [−50, 50] m s⁻¹ with a 5 m s⁻¹ interval. Potential temperature (dashed lines), contours in the range [270, 800] K with a 10 K interval.</figcaption></figure>

ically amounts only to around ≈ 5 – 10% and, for the test cases in Section 4.1 and Section 4.3, the simulations using high-order mappings can be even faster. This is likely due to the fact that the use of higher order quadrature formulas is limited to the curved boundary, namely on a very small number of faces. Moreover, an efficient implementation of quadrature rules in the deal.II library helps to avoid overhead
[24] W.J. Gordon and C.A. Hall. "Transfinite element methods: blending-function interpolation over arbitrary curved element domains". *Numerische Mathematik* 21.2 (1973), pp. 109–129.

[25] J.E. Guerra. "A Development Framework for High Fidelity Non-hydrostatic Simulations of Cross Mountain". PhD thesis. University of California, Davis, 2018.

[26] M.E. Gurtin. *An introduction to continuum mechanics*. Academic press, 1982.

[27] Jan S Hesthaven and Tim Warburton. "Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations". *Journal of Computational Physics* 181.1 (2002), pp. 186–221.

[28] F. Hindenlang, T. Bolemann, and C.-D. Munz. "Mesh Curving Techniques for High Order Discontinuous Galerkin Simulations". Vol. 128. Jan. 2015, pp. 133–152. ISBN: 978-3-319-12885-6.

[29] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. "Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement". *Computer methods in applied mechanics and engineering* 194.39-41 (2005), pp. 4135–4195.

[30] T. Kanehama et al. "Which orographic scales matter most for medium-range forecast skill in the Northern Hemisphere winter?" *Journal of Advances in Modeling Earth Systems* 11.12 (2019), pp. 3893–3910.

[31] J.B. Klemp and D.R. Durran. "An Upper Boundary Condition Permitting Internal Gravity Wave Radiation in Numerical Mesoscale Models". *Monthly Weather Review* 111 (1983), pp. 430–444.

[32] J.B. Klemp and D.K. Lilly. "Numerical simulation of hydrostatic mountain waves". *Journal of the Atmospheric Sciences* 35 (1978), pp. 78–107.

[33] D.A. Kopriva. *Implementing spectral methods for partial differential equations: Algorithms for scientists and engineers*. Springer Science & Business Media, 2009.

[34] L. Krivodonova and M. Berger. "High-order accurate implementation of solid wall boundary conditions in curved geometries". *Journal of Computational Physics* 211.2 (2006), pp. 492–512.

[35] I. Lomtev, R.M. Kirby, and G.E. Karniadakis. "A discontinuous Galerkin ALE method for compressible viscous flows in moving domains". *Journal of Computational Physics* 155.1 (1999), pp. 128–159.

[36] J.F. Louis. "A parametric model of vertical eddy fluxes in the atmosphere." *Boundary Layer Meteorology* 17 (1979), pp. 197–202.

[37] T. Melvin et al. "A mixed finite-element, finite-volume, semi-implicit discretization for atmospheric dynamics: Cartesian geometry". *Quarterly Journal of the Royal Meteorological Society* 145 (2019), pp. 2835–2853.
<figure><img src="image_22.png" /><figcaption>Figure 20: T-REX mountain wave test case, comparison of normalized momentum flux at $t = T_f = 4$ h computed with (19). Results using polynomial degree 3 mapping (blue dots), linear mapping (red line), and reference results with a 400 × 52 elements mesh and a linear mapping (black line). The momentum flux is normalized by the value at $z = 5$ km obtained with the 400 × 52 elements mesh.</figcaption></figure>

when high-order mappings are employed. More specifically, a variable $JxW$, which stands for “Jacobian determinant times weight”, is stored. The two factors, i.e. the determinant of the Jacobian and the weights of the quadrature formula, are always coupled and, therefore, storing only their product allows an efficient evaluation of integrals and saves computational operations. We refer to [5] for further details on the specific implementation in the library. Hence, for $\hat{r} = 2$, the computational cost is the same of the linear mapping (see Table 2), while for $\hat{r} > 2$ the overhead in terms of wall-clock time is reduced and, because of caching effects, the computational cost can be even lower.

The only test cases in which a sizeable overhead appears are the unfiltered non-smooth topography with perturbation factor $\delta = 0.15$ in (25) and, with a reduced overhead, the T-REX mountain wave (Tables 6 and 9). This result is likely due to the fact that the high-order mapping captures and resolves small-scale flow patterns that lead to a more complex flow field and strongly affect the development of lee waves (see Figure 12). This is confirmed by the fact that the overhead factor is significantly reduced employing a filtered orography. These considerations further support the use of high-order mappings, which significantly outperform linear mapping in terms of accuracy with an increase in computational cost that only depends on the more complex flow features induced by the better resolved orography. However, the wall-clock time needed for the high-order mapping is more than twice smaller than the wall-clock needed for the reference simulation employing the linear mapping.
[38] M.J. Miller, T.N. Palmer, and R. Swinbank. "Parametrization and influence of subgridscale orography in general circulation and numerical weather prediction models". *Meteorology and Atmospheric Physics* 40 (1989), pp. 84–109.

[39] D. Moxey et al. "High-order curvilinear meshing using a thermo-elastic analogy". *Computer-Aided Design* 72 (2016), pp. 130–139.

[40] G. Orlando. "Modelling and simulations of two-phase flows including geometric variables". http://hdl.handle.net/10589/198599. PhD thesis. Politecnico di Milano, 2023.

[41] G. Orlando, P.F. Barbante, and L. Bonaventura. "An efficient IMEX-DG solver for the compressible Navier-Stokes equations for non-ideal gases". *Journal of Computational Physics* 471 (2022), p. 111653.

[42] G. Orlando, T. Benacchio, and L. Bonaventura. "An IMEX-DG solver for atmospheric dynamics simulations with adaptive mesh refinement". *Journal of Computational and Applied Mathematics* 427 (2023), p. 115124.

[43] G. Orlando, T. Benacchio, and L. Bonaventura. *Robust and accurate simulations of flows over orography using non-conforming meshes*. 2024. arXiv: 2402.07759 [physics.ao-ph].

[44] T.N. Palmer, G.J. Shutts, and R. Swinbank. "Alleviation of a systematic westerly bias in general circulation and numerical weather prediction models through an orographic gravity wave drag parametrization". *Quarterly Journal of the Royal Meteorological Society* 112 (1986), pp. 1001–1039.

[45] O. Sahni et al. "Curved boundary layer meshing for adaptive viscous flow simulations". *Finite Elements in Analysis and Design* 46.1-2 (2010), pp. 132–139.

[46] C.M. Schär et al. "A new terrain-following vertical coordinate formulation for atmospheric prediction models". *Monthly Weather Review* 130 (2002), pp. 2459–2480.

[47] R.B. Smith. "The Influence of Mountains on the Atmosphere". *Advances in Geophysics* 21 (1979), pp. 87–230.

[48] G. Strang and A.E. Berger. "The change in solution due to change in domain". 1973, pp. 199–205.

[49] V. Thomée. *Polygonal domain approximation in Dirichlet's problem*. Mathematics Research Center, University of Wisconsin, 1971.

[50] T. Toulorge et al. "Robust untangling of curvilinear meshes". *Journal of Computational Physics* 254 (2013), pp. 8–26.

[51] M. Turner, J. Peiró, and D. Moxey. "Curvilinear mesh generation using a variational framework". *Computer-Aided Design* 103 (2018), pp. 73–91.
# 5 Conclusions

The accurate resolution of orographic profiles is of paramount importance for atmospheric applications [48, 50]. Here we have performed a quantitative study on the impact of curved elements for flows over orography. To the best of our knowledge, this study has never been performed for atmospheric applications. We have employed the IMEX-DG solver originally proposed in [46] and validated for atmospheric applications in [47, 48]. The software, implemented in the framework of the open-source library deal.II [2], natively supports high-order polynomials for the mapping between the reference element and elements in the physical space.

Numerical experiments on a number of classical benchmarks of flow over idealised and real orography have shown that, at a given resolution, results obtained with a high-order mapping significantly outperform those obtained using a linear mapping in terms of error with respect to reference solutions. A sizeable increase in computational cost is only observed when more complex flow features are induced by the better resolved orography. In fact, the use of high-order mapping leads to results which are analogous to those obtained using a linear mapping at orography-resolving high resolution. Hence, employing high-order mapping is to some extent equivalent to considering a sub-tessellation to evaluate the integrals. Applications of the use of high-order mappings to non-smooth orography profiles have also been presented, showing that the use of high-order mappings generally provides better results with respect to those obtained with a linear mapping. This is valid also for orography profiles obtained with standard filtering procedures. At a given spatial resolution, high-order mappings capture small-scale topographic features inaccessible to the standard linear mapping. Performing very high resolution simulations is computationally expensive, but it is beneficial for mountain atmospheric processes and forecast skill. The use of high-order mappings can be considered as a valuable alternative tool to resolve orographic features at lower spatial resolutions and at a feasible computational cost. Furthermore, based on the numerical evidence presented in this paper, we believe that only results taking into account improved representations of the geometry should be considered as the reference for future works.

# Acknowledgements

G.O. is part of the INdAM-GNCS National Research Group. The simulations have been partly run at CINECA thanks to the computational resources made available through the ISCRA-C projects FEMTUF - HP10CTQ8X7 and FEM-GPU - HP10CQYKJ1. We acknowledge the CINECA award, for the availability of high-performance computing resources and support. This work has been partly supported by the ESCAPE-2 project, European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 800897). The authors also gratefully acknowledge Dr. Doyle and Dr. Gaberšek for having provided the orography profile and the initial data of the T-REX mountain wave.
[52] Z.J. Wang and H. Gao. "A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids". *Journal of Computational Physics* 228.21 (2009), pp. 8161-8186.

[53] Z.J. Wang and Y. Sun. "Curvature-based wall boundary condition for the Euler equations on unstructured grids". *AIAA journal* 41.1 (2003), pp. 27-33.

[54] S. Webster et al. "Improvements to the representation of orography in the Met Office Unified Model". *Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography* 129.591 (2003), pp. 1989-2010.

[55] J. Yin et al. "A curved boundary treatment for discontinuous Galerkin method applied to Euler equations on triangular and tetrahedral grids". *Computer Physics Communications* 258 (2021), p. 107549.

[56] X. Zhang. "A curved boundary treatment for discontinuous Galerkin schemes solving time dependent problems". *Journal of Computational Physics* 308 (2016), pp. 153-170.

[57] O.C. Zienkiewicz and P.B. Morice. *The finite element method in engineering science*. Vol. 1977. McGraw-hill London, 1971.
References

[1] P.F. Antonietti et al. *lymph: discontinuous poLYtopal methods for Multi-PHysics differential problems*. 2024. arXiv: 2401.13376 [math.NA].

[2] D. Arndt et al. “The deal.II library, Version 9.5”. *Journal of Numerical Mathematics* 31.3 (2023), pp. 231–246.

[3] D.N. Arnold, D. Boffi, and F. Bonizzoni. “Finite element differential forms on curvilinear cubic meshes and their approximation properties”. *Numerische Mathematik* 129.1 (2015), pp. 1–20.

[4] D.N. Arnold, R. Falk, and R. Winther. “Finite element exterior calculus: from Hodge theory to numerical stability”. *Bulletin of the American mathematical society* 47.2 (2010), pp. 281–354.

[5] W. Bangerth, R. Hartmann, and G. Kanschat. “deal II: a general-purpose object-oriented finite element library”. *ACM Transactions on Mathematical Software* 33 (2007), pp. 24–51.

[6] F. Bassi and S. Rebay. “A High-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations”. *Journal of Computational Physics* 131 (1997), pp. 267–279.

[7] L. Bonaventura. “A Semi-Implicit, Semi-Lagrangian Scheme Using the Height Coordinate for a Nonhydrostatic and Fully Elastic Model of Atmospheric Flows.” *Journal of Computational Physics* 158 (2000), pp. 186–213.

[8] L. Bonaventura and R. Ferretti. “Semi-Lagrangian methods for parabolic problems in divergence form.” *SIAM Journal of Scientific Computing* 36 (2014), A2458–A2477.

[9] Walter Boscheri et al. “A geometrically and thermodynamically compatible finite volume scheme for continuum mechanics on unstructured polygonal meshes”. *Journal of Computational Physics* 507 (2024), p. 112957.

[10] Andrea Cangiani, Zhaonan Dong, and Emmanuil Georgoulis. “hp-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements”. *Mathematics of Computation* 91.333 (2022), pp. 1–35.

[11] M. Ciallella et al. “Shifted boundary polynomial corrections for compressible flows: high order on curved domains using linear meshes”. *Applied Mathematics and Computation* 441 (2023), p. 127698.

[12] M. Ciallella et al. *Very high order treatment of embedded curved boundaries in compressible flows: ADER discontinuous Galerkin with a space-time Reconstruction for Off-site data*. 2023. arXiv: 2312.07170 [math.NA].

[13] P.G. Ciarlet. *The finite element method for elliptic problems*. SIAM, 2002.
[14] M. Corti et al. "Numerical modeling of the brain poromechanics by high-order discontinuous Galerkin methods". *Mathematical Models and Methods in Applied Sciences* 33.08 (2023), pp. 1577-1609.

[15] R. Costa et al. "Very high-order accurate finite volume scheme on curved boundaries for the two-dimensional steady-state convection-diffusion equation with Dirichlet condition". *Applied Mathematical Modelling* 54 (2018), pp. 752-767.

[16] R. Costa et al. "Very high-order accurate finite volume scheme for the convection-diffusion equation with general boundary conditions on arbitrary curved boundaries". *International Journal for Numerical Methods in Engineering* 117.2 (2019), pp. 188-220.

[17] C. De Boor. *A practical guide to splines*. Vol. 27. springer-verlag New York, 1978.

[18] S. Dey, R.M. O'Bara, and M.S. Shephard. "Towards curvilinear meshing in 3D: the case of quadratic simplices". *Computer-Aided Design* 33.3 (2001), pp. 199-209.

[19] V. Dolejší. "On the discontinuous Galerkin method for the numerical solution of the Navier-Stokes equations". *International Journal for Numerical Methods in Fluids* 45 (2004), pp. 1083-1106.

[20] J. D. Doyle et al. "An intercomparison of T-REX mountain-wave simulations and implications for mesoscale predictability". *Monthly Weather Review* 139.9 (2011), pp. 2811-2831.

[21] M. Fortunato and P.-O. Persson. "High-order unstructured curved mesh generation using the Winslow equations". *Journal of Computational Physics* 307 (2016), pp. 1-14.

[22] D.C. Fritts et al. "Impacts of limited model resolution on the representation of mountain wave and secondary gravity wave dynamics in local and global models. 1: Mountain waves in the stratosphere and mesosphere". *Journal of Geophysical Research: Atmospheres* 127.9 (2022), e2021JD035990.

[23] I. Fumagalli et al. "Polytopal discontinuous Galerkin discretization of brain multiphysics flow dynamics". *Journal of Computational Physics* 513 (2024), p. 113115.

[24] E. Gaburro et al. "High order entropy preserving ADER-DG schemes". *Applied Mathematics and Computation* 440 (2023), p. 127644.

[25] T. Gal-Chen and R.C.J. Somerville. "On the use of a coordinate transformation for the solution of the Navier-Stokes equations". *Journal of Computational Physics* 17 (1975), pp. 209-228.

[26] F.X. Giraldo. "High-order triangle-based discontinuous Galerkin methods for hyperbolic equations on a rotating sphere". *Journal of Computational Physics* 214.2 (2006), pp. 447-465.
[27] F.X. Giraldo. *An Introduction to Element-Based Galerkin Methods on Tensor-Product Bases*. Springer Nature, 2020.

[28] F.X. Giraldo and M. Restelli. "A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases". *Journal of Computational Physics* 227 (2008), pp. 3849-3877.

[29] W.J. Gordon and C.A. Hall. "Transfinite element methods: blending-function interpolation over arbitrary curved element domains". *Numerische Mathematik* 21.2 (1973), pp. 109-129.

[30] J.E. Guerra. "A Development Framework for High Fidelity Non-hydrostatic Simulations of Cross Mountain". PhD thesis. University of California, Davis, 2018.

[31] M.E. Gurtin. *An introduction to continuum mechanics*. Academic press, 1982.

[32] Jan S Hesthaven and Tim Warburton. "Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations". *Journal of Computational Physics* 181.1 (2002), pp. 186-221.

[33] F. Hindenlang, T. Bolemann, and C.-D. Munz. "Mesh Curving Techniques for High Order Discontinuous Galerkin Simulations". Vol. 128. Jan. 2015, pp. 133-152. ISBN: 978-3-319-12885-6.

[34] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. "Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement". *Computer methods in applied mechanics and engineering* 194.39-41 (2005), pp. 4135-4195.

[35] T. Kanehama et al. "Which orographic scales matter most for medium-range forecast skill in the Northern Hemisphere winter?" *Journal of Advances in Modeling Earth Systems* 11.12 (2019), pp. 3893-3910.

[36] J.B. Klemp and D.R. Durran. "An Upper Boundary Condition Permitting Internal Gravity Wave Radiation in Numerical Mesoscale Models". *Monthly Weather Review* 111 (1983), pp. 430-444.

[37] J.B. Klemp and D.K. Lilly. "Numerical simulation of hydrostatic mountain waves". *Journal of the Atmospheric Sciences* 35 (1978), pp. 78-107.

[38] D.A. Kopriva. *Implementing spectral methods for partial differential equations: Algorithms for scientists and engineers*. Springer Science & Business Media, 2009.

[39] L. Krivodonova and M. Berger. "High-order accurate implementation of solid wall boundary conditions in curved geometries". *Journal of Computational Physics* 211.2 (2006), pp. 492-512.
[40] I. Lomtev, R.M. Kirby, and G.E. Karniadakis. "A discontinuous Galerkin ALE method for compressible viscous flows in moving domains". *Journal of Computational Physics* 155.1 (1999), pp. 128–159.

[41] J.F. Louis. "A parametric model of vertical eddy fluxes in the atmosphere." *Boundary Layer Meteorology* 17 (1979), pp. 197–202.

[42] T. Melvin et al. "A mixed finite-element, finite-volume, semi-implicit discretization for atmospheric dynamics: Cartesian geometry". *Quarterly Journal of the Royal Meteorological Society* 145 (2019), pp. 2835–2853.

[43] M.J. Miller, T.N. Palmer, and R. Swinbank. "Parametrization and influence of subgridscale orography in general circulation and numerical weather prediction models". *Meteorology and Atmospheric Physics* 40 (1989), pp. 84–109.

[44] D. Moxey et al. "High-order curvilinear meshing using a thermo-elastic analogy". *Computer-Aided Design* 72 (2016), pp. 130–139.

[45] G. Orlando. "Modelling and simulations of two-phase flows including geometric variables". http://hdl.handle.net/10589/198599. PhD thesis. Politecnico di Milano, 2023.

[46] G. Orlando, P.F. Barbante, and L. Bonaventura. "An efficient IMEX-DG solver for the compressible Navier-Stokes equations for non-ideal gases". *Journal of Computational Physics* 471 (2022), p. 111653.

[47] G. Orlando, T. Benacchio, and L. Bonaventura. "An IMEX-DG solver for atmospheric dynamics simulations with adaptive mesh refinement". *Journal of Computational and Applied Mathematics* 427 (2023), p. 115124.

[48] G. Orlando, T. Benacchio, and L. Bonaventura. "Robust and accurate simulations of flows over orography using non-conforming meshes" (2024). Accepted for publication in Quarterly Journal of the Royal Meteorological Society. DOI: 10.1002/qj.4839.

[49] T.N. Palmer, G.J. Shutts, and R. Swinbank. "Alleviation of a systematic westerly bias in general circulation and numerical weather prediction models through an orographic gravity wave drag parametrization". *Quarterly Journal of the Royal Meteorological Society* 112 (1986), pp. 1001–1039.

[50] J.M. Prusa and W.J. Gutowski. "MPDATA and grid adaptivity in geo-physical fluid flow models". *International Journal for Numerical Methods in Fluids* 50 (2006), pp. 1207–1228.

[51] O. Sahni et al. "Curved boundary layer meshing for adaptive viscous flow simulations". *Finite Elements in Analysis and Design* 46.1-2 (2010), pp. 132–139.

[52] C.M. Schär et al. "A new terrain-following vertical coordinate formulation for atmospheric prediction models". *Monthly Weather Review* 130 (2002), pp. 2459–2480.
[53] R.B. Smith. "The Influence of Mountains on the Atmosphere". *Advances in Geophysics* 21 (1979), pp. 87-230.

[54] G. Strang and A.E. Berger. "The change in solution due to change in domain". 1973, pp. 199-205.

[55] V. Thomée. *Polygonal domain approximation in Dirichlet's problem*. Mathematics Research Center, University of Wisconsin, 1971.

[56] T. Toulorge et al. "Robust untangling of curvilinear meshes". *Journal of Computational Physics* 254 (2013), pp. 8-26.

[57] M. Turner, J. Peiró, and D. Moxey. "Curvilinear mesh generation using a variational framework". *Computer-Aided Design* 103 (2018), pp. 73-91.

[58] Z.J. Wang and H. Gao. "A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids". *Journal of Computational Physics* 228.21 (2009), pp. 8161-8186.

[59] Z.J. Wang and Y. Sun. "Curvature-based wall boundary condition for the Euler equations on unstructured grids". *AIAA journal* 41.1 (2003), pp. 27-33.

[60] S. Webster et al. "Improvements to the representation of orography in the Met Office Unified Model". *Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography* 129.591 (2003), pp. 1989-2010.

[61] J. Yin et al. "A curved boundary treatment for discontinuous Galerkin method applied to Euler equations on triangular and tetrahedral grids". *Computer Physics Communications* 258 (2021), p. 107549.

[62] X. Zhang. "A curved boundary treatment for discontinuous Galerkin schemes solving time dependent problems". *Journal of Computational Physics* 308 (2016), pp. 153-170.

[63] O.C. Zienkiewicz and P.B. Morice. *The finite element method in engineering science*. Vol. 1977. McGraw-hill London, 1971.