1 Introduction↩︎

The absorption of vertically propagating waves near the upper boundary is a key component of numerical models of atmospheric fluid flow. As solving the governing equations on a computer necessitates a truncated vertical domain, some form of artificial boundary is required at the model top. Models with a fixed height coordinate often apply a rigid lid condition that enforces zero vertical velocity at the stationary model top, whilst methods that use floating Lagrangian pressure levels [1] often apply a free-slip boundary where there is no mass flux through the uppermost pressure level. However, both boundary conditions permit reflections at the model top, causing waves that should have exited the truncated domain to instead propagate downwards and pollute the simulated dynamics.

To mitigate model top wave reflections, most atmospheric models introduce a sponge layer, which damps the waves that enter the uppermost subdomain [2]. Sponge layers often use Rayleigh damping, which relaxes the targeted fields back to a reference state. Alternatively, a Laplacian diffusion type damping might be used, or a combination of Laplacian and Rayleigh damping. However, typical sponge layers are often insufficient in idealised tests, or require many levels to provide effective damping; examples include the 10 km sponge layer required for the mesoscale, orographically driven test cases of [3], or the 30 sponge levels recommended for the mountain test case of [4]. Other approaches for dealing with vertically propagating waves include the modified Rayleigh damping method of [5], or the use of radiation or absorbing boundary conditions, e.g. [6], [7]. However, these types of boundary conditions are often complicated or infeasible to accurately derive for the governing equations, which contain a spectrum of wave types and frequencies. This work proposes and examines a new possibility for absorbing vertically propagating waves for atmospheric flow problems: the perfectly matched layer (PML).

Perfectly matched layers have been widely used for electromagnetic applications since their introduction in [8], [9]. From there, the PML has been applied to many fields of computational modelling, such as aerospace, seismic, photonic and elastodynamical applications; a summary of important contributions for hyperbolic PDEs can be found in [10]. The PML is a layer that damps all waves that enter, irrespective of their wavelength and angle of incidence. The term perfectly matched describes that, in the absence of numerical error, there is no reflection of waves at the boundary between the PML and the rest of the domain. By contrast, Rayleigh damping sponge layers are often very sensitive to the damping strength, with a too large damping coefficient inducing reflection at the damping layer boundary [3], [7]. Whilst the PML requires the introduction of auxiliary variables and equations which couple to the original equation set, it typically only requires a small region to provide effective damping.

It is significant to note that whilst there are instances of PMLs for geophysical fluid applications, these have been in the context of damping waves at the horizontal boundaries of limited area models of the shallow water equations [11]–[13]. A related method for the shallow water equations is the absorbing boundary condition approach of [14]. To our knowledge, the use of a PML to absorb waves in the vertical dimension of atmospheric models has not been investigated. The distinction between horizontal and vertical dynamics is a critical feature of atmospheric dynamics. Where the shallow water equations are the common simplified model for the horizontal dynamics, the equivalent in the vertical dimension is the Boussinesq equations, which examine the impact of stratification on buoyancy oscillations and vertically propagating waves. This motivates our use of the Boussinesq equations in a two-dimensional vertical slice to test the PML for model top damping.

We will work with a compressible version of the Boussinesq equations to test the PML for damping two types of atmospheric waves: i) acoustic waves, and ii) orographic gravity waves, which are a subset of internal gravity waves that are excited by flow over a varying surface height. Although many other types of waves may be present in a full atmospheric model, acoustic and orographic gravity waves are two key types that can reflect off the model top. Acoustic waves operate on a very fast timescale, whilst orographic gravity waves propagate much more slowly. In this work, we provide a single PML formulation to damp both acoustic and orographic gravity waves.

We derive PML equations for both linear and nonlinear versions of the compressible Boussinesq equations. As the PML is essentially a linear mechanism, for the linear equations we simply apply the existing PML theory for linear hyperbolic systems [10]. Applications of the PML to nonlinear hyperbolic systems have received much less attention in the literature, aside from the notable exception of [15]. Thus, the derivation of the PML for the nonlinear Boussinesq equations is a novel feature of this study. To succeed in the nonlinear system, we formulate the PML to absorb perturbations from a hydrostatically balanced state.

Another novelty of this study is the formulation of compatible finite element approximations of the PML systems. The compatible finite element codebase Gusto is used for numerical tests of the PML for the linear and nonlinear Boussinesq equations. Gusto is a toolbox for solving geophysical fluid equations, designed for rapid prototyping and built from the finite element codebase Firedrake [16]. We will use simple timestepping methods in Gusto to most clearly highlight the impact of adding PML terms to the original equations. It is important to emphasise that whilst we use the compatible finite element method in this study, the PML equations could be approximated by any numerical method of choice, such as finite difference, finite volume, or continuous and discontinuous Galerkin finite element methods.

After this introduction, Section 2 provides a derivation of the PML for first-order hyperbolic systems, before describing our specific PML for the linear and nonlinear Boussinesq equations. Section 3 describes the compatible finite element discretisation used for the numerical tests of section 4: A linear Boussinesq test for acoustic wave damping, and a nonlinear test with a mountain to examine acoustic and orographic gravity wave damping. Section 5 ends with conclusions and future research directions.

2 The perfectly matched layer↩︎

This section presents a short derivation of the PML for general linear hyperbolic PDEs using the complex coordinate stretching technique of [9], as outlined in [10]. This approach is then applied to linear Boussinesq equations to form a PML for this system. Subsequently, we extend the PML to the nonlinear Boussinesq equations. Since the PML theory is fundamentally linear, the primary novelty of this section lies in the construction of an accurate PML model in the presence of advective nonlinearities. Last, the parameters of the PML are briefly discussed.

2.1 Derivation for first-order hyperbolic PDEs↩︎

We first consider the application of the PML to a linear, first-order, hyperbolic system of equations. The form of PML we consider here is the complex frequency shifted (CFS) variant introduced by [17] and analysed by [10].

Consider the initial value, linear partial differential equation (PDE) for a prognostic variable set \(\mathbf{U}\), \[\frac{\partial{\mathbf{U}}}{\partial{t}} + A \frac{\partial{\mathbf{U}}}{\partial{x}} + B \frac{\partial{\mathbf{U}}}{\partial{z}} + C \mathbf{U} = \mathbf{0}, ~\mathbf{U}(x,z,t=0) = \mathbf{U}_0,\]

defined on a vertical slice domain with Cartesian coordinates \((x,z)\), with matrices \(A,B,\) and \(C\) containing constant coefficients. We seek to damp waves propagating vertically, i.e. in the \(z\)-coordinate.

Consider a Laplace transformation of the prognostic variables from physical space to the complex frequency domain, defined by \[\mathbf{\widetilde{U}}(x,z,s) = \int_{t=0}^{\infty} e^{-st} \mathbf{U}(x,z,t) ~\text{d}t, ~s \in \mathbb{C}, ~\text{Re} \{s\} > 0.\]

The Laplace-transformed PDE system is \[s \mathbf{\widetilde{U}} + A \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{x}} + B \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{z}} + C \mathbf{\widetilde{U}} = \mathbf{0}, \label{eq:laplace95transformed95pde}\tag{1}\]

where we have tacitly assumed homogenous initial conditions.

Next, we introduce a scaled vertical coordinate of \(\zeta=\int_0^z S_z(z') ~\text{d}z'\), with the PML complex metric \[S_z(z) = \gamma_z(z) \left( 1 + \frac{\sigma(z)}{s + \alpha} \right). \label{eq:Sz}\tag{2}\]

Here, \(\gamma_z(z)>0\) is a real coordinate stretching [10], [18], [19], \(\sigma(z)\ge 0\) is the PML damping function that is only nonzero in the PML damping region near the model top, and \(\alpha \geq 0\) is a real-valued constant defining a complex frequency shift which improves stability of the PML [10], [17], [18], [20]. Using the coordinate transform (2 ) leads to \[s \mathbf{\widetilde{U}} + A \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{x}} + \frac{1}{S_z}B \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{\zeta}} + C \mathbf{\widetilde{U}} = \mathbf{0}.\]

Now, we introduce a vector of auxiliary PML variables, \(\mathbf{Q}\). The auxiliary PML variables are defined in the complex frequency domain as \[\mathbf{\widetilde{Q}} = \frac{1}{(s + \alpha)S_z} B \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{\zeta}}.\]

We then identify a coupled system of equations for the original and PML variables in the complex frequency domain, \[\begin{align} s \mathbf{\widetilde{U}} + A \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{x}} + \frac{1}{\gamma_z}B \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{\zeta}} - \sigma \mathbf{\widetilde{Q}} + C \mathbf{\widetilde{U}} &= \mathbf{0}, \\ s \mathbf{\widetilde{Q}} - \frac{1}{\gamma_z} B \frac{\partial{\mathbf{\widetilde{U}}}}{\partial{\zeta}} + (\sigma + \alpha) \mathbf{\widetilde{Q}} &= \mathbf{0}, \end{align}\]

using the identity \[\frac{1}{S_z} = \frac{1}{\gamma_z} - \frac{1}{S_z} \frac{\sigma}{s + \alpha},\]

which follows from the definition of \(S_z\) (2 ). Performing the inverse Laplace transform yields the PML equations in physical space: \[\begin{align} \frac{\partial{\mathbf{U}}}{\partial{t}} + A \frac{\partial{\mathbf{U}}}{\partial{x}} + \frac{1}{\gamma_z} B \frac{\partial{\mathbf{U}}}{\partial{\zeta}} + C \mathbf{U} - \sigma \mathbf{Q} &= \mathbf{0}, ~\mathbf{U}(x,z,t=0) = \mathbf{U}_0, \tag{3} \\ \frac{\partial{\mathbf{Q}}}{\partial{t}} - \frac{1}{\gamma_z} B \frac{\partial{\mathbf{U}}}{\partial{\zeta}} + (\sigma + \alpha) \mathbf{Q} &= \mathbf{0}, ~\mathbf{Q}(x,z,t=0) = \mathbf{0}. \tag{4} \end{align} \tag{5}\]

Observe that the initial condition for all PML variables is zero.

We introduce two important variants of the gradient operator for use in the PML equations, which are a scaled gradient of \({\boldsymbol{\nabla_S}}\), and a scaled gradient that only acts in the vertical dimension, \({\boldsymbol{\nabla_V}}= (\hat{\boldsymbol{\mathrm{k}}}\cdot {\boldsymbol{\nabla_S}})\hat{\boldsymbol{\mathrm{k}}}\), where \(\hat{\boldsymbol{\mathrm{k}}}\) denotes the unit vector in the \(z\)-coordinate. These operators are \[{\boldsymbol{\nabla_S}}= \left[ \frac{\partial{}}{\partial{x}},\frac{1}{\gamma_z}\frac{\partial{}}{\partial{z}} \right]^T, ~{\boldsymbol{\nabla_V}}= \left[ 0,\frac{1}{\gamma_z}\frac{\partial{}}{\partial{z}} \right]^T. \label{eq:PML95grad}\tag{6}\]

2.2 A PML for the Boussinesq system↩︎

We now use the theory for linear hyperbolic systems to derive a PML for the linear Boussinesq equations, then discuss the modifications required to apply this to the nonlinear Boussinesq equations. The prognostic variables of the Boussinesq system are a velocity vector of \(\mathbf{u} = [u,w]\), for horizontal velocity \(u\) and vertical velocity \(w=\mathbf{u} \cdot \mathbf{\hat{k}}\), a kinematic pressure \(p\) (units of m\(^2\) s\(^{-2}\)), and a buoyancy of \(b\) (units of m s\(^{-2}\)). These pressure and buoyancy fields are deviations from a background, hydrostatically balanced state that is used to derive the Boussinesq equations. There are two parameters in this system: the sound speed \(c_{\text{s}}\), and the Brunt-Väisälä frequency \(N\), which measures the timescale of the buoyancy oscillations.

2.2.1 Linear Boussinesq equations↩︎

We first apply the PML to a linearised version of the compressible Boussinesq equations, which support acoustic and gravity waves. These equations are described in Section 8.2 of [21] and were used by [22] to investigate the choice of vertical staggering for the compatible finite element dynamical core, GungHo [23], [24]. The linear Boussinesq equations are \[\begin{align} \frac{\partial{\mathbf{u}}}{\partial{t}} + {\boldsymbol{\nabla}{p}} - b \mathbf{\hat{k}} &= \mathbf{0}, \tag{7} \\ \frac{\partial{p}}{\partial{t}} + c_{\text{s}}^2 (\boldsymbol{\nabla} \cdot \mathbf{u}) &= 0 \tag{8}, \\ \frac{\partial{b}}{\partial{t}} + N^2 (\mathbf{u} \cdot \mathbf{\hat{k}}) &= 0. \tag{9} \end{align} \tag{10}\]

We now derive a PML for the linear Boussinesq equations (10 ). As \(z\) derivatives only appear in the equations for \(\mathbf{u}\) and \(p\), we only require the introduction of two PML variables, \(\mathbf{q_u}\) and \(q_p\). This results in the PML equations of \[\begin{align} \frac{\partial{\mathbf{u}}}{\partial{t}} + \boldsymbol{\nabla_S} p - b \mathbf{\hat{k}} - \sigma \mathbf{q_u} &= \mathbf{0}, \\ \frac{\partial{p}}{\partial{t}} + c_{\text{s}}^2 (\boldsymbol{\nabla_S} \cdot \mathbf{u}) - \sigma q_p &= 0, \\ \frac{\partial{b}}{\partial{t}} + N^2 (\mathbf{u} \cdot \mathbf{\hat{k}}) &= 0, \\ \frac{\partial{\mathbf{q_u}}}{\partial{t}} - \boldsymbol{\nabla_V} p + (\sigma + \alpha) \mathbf{q_u}&= \mathbf{0}, \tag{11} \\ \frac{\partial{q_p}}{\partial{t}} - c_{\text{s}}^2 (\boldsymbol{\nabla_V} \cdot \mathbf{u}) + (\sigma + \alpha) q_p&= 0 \tag{12}, \end{align} \tag{13}\]

with the introduction of the scaled PML gradient operators of \({\boldsymbol{\nabla_S}}\) and \({\boldsymbol{\nabla_V}}\), defined in 6 . The PML terms of \(- \boldsymbol{\nabla_V} p\) and \(-c_{\text{s}}^2 (\boldsymbol{\nabla_V} \cdot \mathbf{u})\) damp the vertical wave propagation constructed by the pressure gradient and divergence, which are the key drivers of acoustic wave motion.

2.2.2 Nonlinear Boussinesq equations↩︎

We now move to the nonlinear Boussinesq system. To apply the theory for linear hyperbolic PDEs, the PML is applied to perturbations from a hydrostatically balanced base state of the Boussinesq fields (noting that the Boussinesq fields are themselves perturbations from a background hydrostatic state), derived from linearisation of the nonlinear equations. We use the compressible form of the Boussinesq equations as described by [25], which permits acoustic waves that are not present with the standard Boussinesq approximation.

The nonlinear Boussinesq equations are \[\begin{align} \frac{\partial{\mathbf{u}}}{\partial{t}} + (\mathbf{u} \cdot {\boldsymbol{\nabla}{)}} \mathbf{u} + {\boldsymbol{\nabla}{p}} - b \mathbf{\hat{k}} &= \mathbf{0}, \\ \frac{\partial{p}}{\partial{t}} + (\mathbf{u} \cdot {\boldsymbol{\nabla}{)}} p + c_{\text{s}}^2 (\boldsymbol{\nabla} \cdot \mathbf{u}) &= 0, \\ \frac{\partial{b}}{\partial{t}} + (\mathbf{u} \cdot {\boldsymbol{\nabla}{)}} b &= 0. \end{align} \label{eq:bous}\tag{14}\]

The main difference from the linear Boussinesq equations (14 ) is the presence of nonlinear advective transport terms for each prognostic variable. In the buoyancy equation, the advective term replaces the constant buoyancy frequency term of \(N^2 (\mathbf{u} \cdot \hat{\boldsymbol{\mathrm{k}}})\) (9 ).

To linearise the nonlinear Boussinesq equations, we decompose the variables into mean and perturbation components, which are denoted by overbars and primes, respectively: \[\mathbf{u}(x,z) = \mathbf{u}^{'}, ~p(x,z) = \overline{p}(z) + p^{'}, ~b(x,z) = \overline{b}(z) + b^{'}.\]

Note that the velocity is linearised about a state of rest.

The linearisation of the nonlinear Boussinesq equations is then \[\begin{align} \frac{\partial{\mathbf{u}^{'}}}{\partial{t}} + {\boldsymbol{\nabla}{p}}^{'} - b^{'} \mathbf{\hat{k}} &= \mathbf{0}, \\ \frac{\partial{p^{'}}}{\partial{t}} + (\mathbf{u}^{'} \cdot {\boldsymbol{\nabla}{)}} \overline{p} + c_{\text{s}}^2 (\boldsymbol{\nabla} \cdot \mathbf{u}^{'}) &= 0, \\ \frac{\partial{b^{'}}}{\partial{t}} + (\mathbf{u}^{'} \cdot {\boldsymbol{\nabla}{)}} \overline{b} &= 0. \end{align} \label{eq:bous95linearised}\tag{15}\]

We apply PML damping to the pressure gradient and divergence terms, as with the linear Boussinesq system, which means that a PML buoyancy variable is again not required. Although the advective terms contain \(z\) derivatives corresponding to vertical transport, we neglect these in this study for a couple of reasons. First, the application of the PML to advective terms is often challenging and can be prone to instability, e.g. [26]–[29]. Second, the vertical velocity is typically much smaller than the horizontal velocity, especially near the model top, where no vertical velocity may be prescribed at the upper boundary. Hence, we neglect the vertical transport terms for this proof-of-concept use of the PML, but will revisit this assumption in future work.

The nonlinear Boussinesq PML equations are then \[\begin{align} \frac{\partial{\mathbf{u}}}{\partial{t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} + \boldsymbol{\nabla_S} p^{'} - b^{'} \mathbf{\hat{k}} - \sigma \mathbf{q_u}&= \mathbf{0}, \tag{16} \\ \frac{\partial{p}}{\partial{t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) p + c_{\text{s}}^2 (\boldsymbol{\nabla_S} \cdot \mathbf{u}) - \sigma q_p&= 0, \\ \frac{\partial{b}}{\partial{t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) b &= 0, \\ \frac{\partial{\mathbf{q_u}}}{\partial{t}} - \boldsymbol{\nabla_V} p^{'} + (\sigma + \alpha) \mathbf{q_u}&= \mathbf{0}, \\ \frac{\partial{q_p}}{\partial{t}} - c_{\text{s}}^2 (\boldsymbol{\nabla_V} \cdot \mathbf{u}) + (\sigma + \alpha) q_p&= 0. \end{align} \tag{17}\]

Note the explicit use of hydrostatic balance of the reference state, \(\boldsymbol{\nabla} \overline{p} = \overline{b}\mathbf{\hat{k}}\), to simplify the momentum equation (16 ). Setting \(\sigma(z) = 0\) restores the original nonlinear Boussinesq equations (14 ).

2.3 PML parameters↩︎

We now describe the key parameters in the PML, which are summarised in Table 1. We consider a PML with a sine-squared damping profile to match the most common choice for Rayleigh damping in dynamical cores. Other choices include the cubic profile that is common in PML literature [10], or quadratic or hyperbolic tangent functions that have alternatively been used in dynamical cores [2].

Our PML damping function is defined as \[\label{eq:damp95coeff} \begin{align} &\sigma(z) = \begin{dcases} 0, & \text{if} \quad z < z_{\text{damp}},\\ \sigma_0 \sin^2 \left( \frac{\pi}{2}\left( \frac{z - z_{\text{damp}}}{\delta} \right) \right), & \text{if} \quad z \geq z_{\text{damp}} , \end{dcases} \end{align}\tag{18}\]

with \(z_{\text{damp}} = H_z - \delta\) the height at which the PML starts, with \(H_z\) the model top height and \(\delta\) the thickness of the damping layer, and \(\sigma_0\) an inverse damping timescale defined as \[\sigma_0 = \frac{2c_{\max}}{\delta} \ln \left( \frac{1}{\beta} \right), \label{eq:sigma950}\tag{19}\]

with \(c_{\text{max}}\) the fastest resolvable wave speed in the domain and \(\beta\) a tolerance for the domain truncation error. We set \(c_{\text{max}}\) equal to the acoustic wave speed, \(c_{\text{s}}\). Our default PML parameters (Table 1) and a model top at \(H_z = 20 ~\text{km}\) lead to a damping strength of \(\sigma_0 = 2.42 ~\text{s}^{-1}\). A good choice of \(\sigma_0\) balances being large enough to successfully damp incoming waves, but not too large to degrade the numerical accuracy.

For the stretched PML vertical coordinate (2 ), we use \[\gamma_z(z) = 1 + \gamma_0 \sigma(z), \label{eq:gamma95z}\tag{20}\]

with \(\gamma_0 \geq 0\in \mathbb{R}\) a real grid stretching factor (with units of s) and \(\gamma_0 = 0 ~\text{s}\) reducing to an unstretched grid. A non-zero \(\gamma_0\) typically improves the damping ability of the PML, but can lead to accuracy reduction if too large [30]. We note that both the unstretched (\(\gamma_0=0 ~\text{s}\)) and stretched (\(\gamma_0>0 ~\text{s}\)) PMLs use the same vertical coordinate \(z\) in the computational mesh. The grid stretching modifies the transformation in the complex frequency domain, but does not require a modification of the physical domain. Hence, varying the grid stretching factor is a simple modification that only scales derivatives, as reflected in the modified PML operators of (6 ).

The other key PML parameter is the complex frequency shift (CFS) parameter, \(\alpha\), which only appears in the equations for the auxiliary PML variables (4 ). Increasing \(\alpha\) improves the stability of the PML but reduces its absorption properties [31]. Following [32], we set \(\alpha=0.05 \sigma_0 ~\text{s}^{-1} = 0.121 ~\text{s}^{-1}\).

Our numerical tests, using a finite element method, will use elements of height \(\Delta z = 500 ~\text{m}\). This means that for a model top at \(H_z = 20 ~\text{km}\) and PML thickness of \(\delta = 0.1 H_z = 2 ~\text{km}\), there are only four elements in the vertical extent of the PML. As the auxiliary PML variables only couple to the prognostic equations in the PML region, the extra computation with a PML is restricted to only a small subset of the domain. However, for ease of implementation in these initial studies, we define PML variables that span the entire domain. We will term the region of the domain excluding the PML, \(z\in [0,18] ~\text{km}\), the undamped domain.

3 Numerical methods↩︎

The PML derived for the Boussinesq equations can be approximated by any numerical method of choice, such as finite difference, finite volume or continuous and discontinuous Galerkin finite element methods. For testing purposes, we will discretise the PML using the compatible finite element method. This section overviews the numerical model and weak forms of the PML Boussinesq equations.

3.1 Spatial discretisation↩︎

We will work with a compatible finite element model containing quadrilateral cells and next-to-lowest-order elements, which allows for second-order accuracy. A Charney-Phillips vertical staggering is applied, where thermodynamic variables are offset from the dry density in the vertical coordinate. We select the Charney-Phillips grid over the Lorenz grid, where the dry density and thermodynamic variables are instead colocated, to avoid Lorenz grid computational modes [22], [33]. An advantage of using the compatible finite element method is that it does not require explicit diffusion for stability [34], [35], meaning that the damping layer can be tested in the absence of additional diffusion mechanisms.

The quadrilateral elements are constructed as the tensor product, \(\otimes\), of a one-dimensional finite element in the horizontal with a one-dimensional finite element in the vertical [36]. For degree \(l\), we denote a discontinuous element as \(\text{d}Q_l\) and a continuous element as \(Q_l\). We then use the following types of function space (visualised in Figure 1):

• \(\mathbb{V}_{\mathbf{u}}:= \text{RT}_1\), the next-to-lowest-order Raviart-Thomas function space on quadrilateral elements, which enforces continuity of normal components of the vector field at the facets, whilst the tangential component is permitted to be discontinuous. The horizontal component of the vector field is represented by continuous order 2 polynomials in the horizontal and discontinuous order 1 polynomials in the vertical (\(Q_2 \otimes \text{d}Q_1\)). Conversely, the vertical component is represented by discontinuous order 1 polynomials in the horizontal and continuous order 2 polynomials in the vertical (\(\text{d}Q_1 \otimes Q_2\)).

• \(\mathbb{V}_{p} := \text{d}Q_1 \otimes \text{d}Q_1\), the fully discontinuous next-to-lowest order function space.

• \(\mathbb{V}_{b} := \text{d}Q_1 \otimes Q_2\), the next-to-lowest-order Charney-Phillips function space for thermodynamic variables, i.e. buoyancy, which is discontinuous in the horizontal and continuous in the vertical.

For a Charney-Phillips staggering, we take \((\mathbf{u},p,b) \in (\mathbb{V}_{\mathbf{u}}, \mathbb{V}_{p}, \mathbb{V}_{b})\). We define the PML variables in the same function space as their prognostic counterpart, so \(\mathbf{q_u} \in \mathbb{V}_{\mathbf{u}}, q_p \in \mathbb{V}_{p}\).

3.2 Weak forms↩︎

Here we present the weak forms used for our compatible finite element discretisations. We will consider first the weak form of the PML for the linear Boussinesq equations, then proceed to the nonlinear Boussinesq equations.

3.2.1 Linear Boussinesq equations↩︎

To obtain a weak form of the linear Boussinesq PML equations (17 ), we take the inner product of each equation with a corresponding test function, using \(\tau\) to denote prognostic test functions and \(\chi\) for PML test functions. Integrating over the domain, \(\Omega\), arrives at \[\begin{align} \int_{\Omega} \boldsymbol{\tau_u} \cdot \frac{\partial{\mathbf{u}}}{\partial{t}} ~ \text{d} \boldsymbol{x} - \int_{\Omega} p {\boldsymbol{\nabla}{\cdot }}(\boldsymbol{\tau_u} \odot \boldsymbol{\gamma_S}) ~ \text{d} \boldsymbol{x} - \int_{\Omega} \boldsymbol{\tau_u} \cdot b \hat{\mathbf{k}} ~ \text{d} \boldsymbol{x}& \nonumber \\ - \int_{\Omega} \boldsymbol{\tau_u} \cdot \sigma \mathbf{q_u} ~ \text{d} \boldsymbol{x}&= 0, ~\forall \boldsymbol{\tau_u} \in \mathbb{V}_{\mathbf{u}}, \tag{21} \\ \int_{\Omega} \tau_p \frac{\partial{p}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Omega} c_{\text{s}}^2 \tau_p (\boldsymbol{\nabla_S} \cdot \mathbf{u}) ~ \text{d} \boldsymbol{x} - \int_{\Omega} \sigma \tau_p q_p ~ \text{d} \boldsymbol{x} &= 0, ~\forall \tau_p \in \mathbb{V}_{p}, \tag{22} \\ \int_{\Omega} \tau_b \frac{\partial{b}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Omega} \tau_b N^2 (\mathbf{u} \cdot \hat{\mathbf{k}}) ~ \text{d} \boldsymbol{x} &= 0, ~\forall \tau_b \in \mathbb{V}_{b}, \tag{23} \\ \int_{\Omega} \boldsymbol{\chi_u} \cdot \frac{\partial{\mathbf{q_u}}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Omega} p {\boldsymbol{\nabla}{\cdot }}(\boldsymbol{\chi_u} \odot \boldsymbol{\gamma_V)} ~ \text{d} \boldsymbol{x} + \int_{\Omega} \boldsymbol{\chi_u} \cdot (\sigma + \alpha) \mathbf{q_u} ~ \text{d} \boldsymbol{x} &= 0, ~\forall \boldsymbol{\chi_u} \in \mathbb{V}_{\mathbf{u}}, \\ \int_{\Omega} \chi_p \frac{\partial{q_p}}{\partial{t}} ~ \text{d} \boldsymbol{x} - \int_{\Omega} c_{\text{s}}^2 \chi_p (\boldsymbol{\nabla_V} \cdot \mathbf{u}) ~ \text{d} \boldsymbol{x} + \int_{\Omega} (\sigma + \alpha) \chi_p q_p ~ \text{d} \boldsymbol{x} &= 0, ~\forall \chi_p \in \mathbb{V}_{p}, \end{align} \tag{24}\]

where \(\boldsymbol{x}\) denotes the spatial coordinate. The application of integration by parts on the pressure gradient means that the PML grid stretching term \(\gamma_z(z)\) is operated upon by the divergence operator. We express this by introducing PML scaling terms of \(\boldsymbol{\gamma_S} = [1, \gamma_z]^T\) and \(\boldsymbol{\gamma_V} = [0, \gamma_z]^T\), and the Hadamard (element-wise) product of two vectors, denoted as \(\odot\) (note that \({\boldsymbol{\nabla_S}}= \boldsymbol{\nabla} \odot \boldsymbol{\gamma_S}, {\boldsymbol{\nabla_V}}= \boldsymbol{\nabla} \odot \boldsymbol{\gamma_V}\)). There are no surface terms left after integration by parts due to the continuity of normal components of vector fields in the Raviart-Thomas function space (Figure 1).

The linear Boussinesq equations conserve the energy \[\mathcal{H} = \frac{1}{2} \int_{\Omega} \left( \mathbf{u} \cdot \mathbf{u} + c_{\text{s}}^{-2} p^2 + N^{-2} b^2 \right) ~ \text{d} \boldsymbol{x}, \label{eq:energy}\tag{25}\]

which can be decomposed into kinetic energy (\(E_{\text{K}}\)), internal energy (\(E_{\text{I}}\)) and potential energy (\(E_{\text{P}}\)) components: \[\begin{align} E_{\text{K}} &= \frac{1}{2} \int_{\Omega} \mathbf{u} \cdot \mathbf{u} ~ \text{d} \boldsymbol{x}, \tag{26} \\ E_{\text{I}} &= \frac{1}{2} \int_{\Omega} c_{\text{s}}^{-2} p^2 ~ \text{d} \boldsymbol{x}, \tag{27} \\ E_{\text{P}} &= \frac{1}{2} \int_{\Omega} N^{-2} b^2 ~ \text{d} \boldsymbol{x}. \tag{28} \end{align} \tag{29}\]

To examine energy conservation in the Boussinesq equations, we take the weak form PML system and scale the pressure equation (22 ) by \(c_{\text{s}}^{-2}\) and the buoyancy equation by \(N^{-2}\) (23 ), so that \[\begin{align} \int_{\Omega} \boldsymbol{\tau_u} \cdot \frac{\partial{\mathbf{u}}}{\partial{t}} ~ \text{d} \boldsymbol{x} - \int_{\Omega} p {\boldsymbol{\nabla}{\cdot }}(\boldsymbol{\tau_u} \odot \boldsymbol{\gamma_S}) ~ \text{d} \boldsymbol{x} - \int_{\Omega} \boldsymbol{\tau_u} \cdot b \hat{\mathbf{k}} ~ \text{d} \boldsymbol{x} &= \int_{\Omega} \boldsymbol{\tau_u} \cdot \sigma\mathbf{q_u} ~ \text{d} \boldsymbol{x}, &~\forall \boldsymbol{\tau_u} \in \mathbb{V}_{\mathbf{u}}, \\ \int_{\Omega} c_{\text{s}}^{-2} \tau_p \frac{\partial{p}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Omega} \tau_p (\boldsymbol{\nabla_S} \cdot \mathbf{u}) ~ \text{d} \boldsymbol{x} &= \int_{\Omega} c_{\text{s}}^{-2}\sigma \tau_p q_p ~ \text{d} \boldsymbol{x}, &~\forall \tau_p \in \mathbb{V}_{p}, \\ \int_{\Omega} N^{-2} \tau_b \frac{\partial{b}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Omega} \tau_b (\mathbf{u} \cdot \hat{\mathbf{k}}) ~ \text{d} \boldsymbol{x} &= 0, &~\forall \tau_b \in \mathbb{V}_{b}. \end{align}\]

Choosing test functions of \(\boldsymbol{\tau_u} = \mathbf{u}, ~\tau_p=p,\) and \(\tau_b=b\) shows that the energy (25 ) is conserved in time in the case of no PML, i.e. \(d\mathcal{H}/dt =0\) when \(\sigma=0\). When including the PML, we have a source term which modifies the energy as \[\frac{\mathrm{d} \mathcal{H}}{\mathrm{d} t} = P_{\text{PML}} = \int_{\Omega} \sigma (\mathbf{u} \cdot \mathbf{q_u} + c_{\text{s}}^{-2} p q_p) ~ \text{d} \boldsymbol{x}, \label{eq:P95PML}\tag{30}\]

with \(P_{\text{PML}}\) denoting the PML power. Energy-conserving formulations for mixed finite element discretisations have been presented previously for other geophysical systems such as the compressible Euler equations [35], [37], [38].

3.2.2 Nonlinear Boussinesq equations↩︎

For the weak form of the nonlinear Boussinesq PML equations (17 ), discontinuous Galerkin (DG) upwinding is applied to the advective terms [39], [40], \[\begin{align} \int_{\Omega} \boldsymbol{\tau_u} \cdot \frac{\partial{\mathbf{u}}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Gamma} (\mathbf{u} \cdot \mathbf{\hat{n}^+}) \mathbf{u}^{\dagger} \cdot [[\boldsymbol{\tau_u}]] ~ \text{d}S & \nonumber \\ - \int_{\Omega} \mathbf{u} \cdot [\nabla \cdot (\boldsymbol{\tau_u} \otimes \mathbf{u})] ~ \text{d} \boldsymbol{x} - \int_{\Omega} p^{'} \boldsymbol{\nabla} \cdot (\boldsymbol{\tau_u} \odot \boldsymbol{\gamma_S}) ~ \text{d} \boldsymbol{x} \nonumber \\ - \int_{\Omega} \boldsymbol{\tau_u} \cdot b^{'}\hat{\boldsymbol{\mathrm{k}}} ~ \text{d} \boldsymbol{x} - \int_{\Omega} \boldsymbol{\tau_u} \cdot \sigma\mathbf{q_u} ~ \text{d} \boldsymbol{x} &= 0, ~\forall \boldsymbol{\tau_u} \in \mathbb{V}_{\mathbf{u}}, \tag{31}\\ \int_{\Omega} \tau_p \frac{\partial{p}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Gamma} (\mathbf{u} \cdot \mathbf{\hat{n}}^+) p^{\dagger} [[\tau_p]] ~ \text{d}S - \int_{\Omega} p \boldsymbol{\nabla} \cdot (\tau_p \mathbf{u}) ~ \text{d} \boldsymbol{x} & \nonumber \\ + \int_{\Omega} c_{\text{s}}^2 \tau_p (\boldsymbol{\nabla_S} \cdot \mathbf{u}) ~ \text{d} \boldsymbol{x} - \int_{\Omega} \sigma \tau_p q_p ~ \text{d} \boldsymbol{x} &= 0, ~\forall \tau_p \in \mathbb{V}_{p}, \\ \int_{\Omega} \tau_b \frac{\partial{b}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Gamma} (\mathbf{u} \cdot \mathbf{\hat{n}}^+) b^{\dagger} [[\tau_b]] ~ \text{d}S - \int_{\Omega} b \boldsymbol{\nabla} \cdot (\tau_b \mathbf{u}) ~ \text{d} \boldsymbol{x} &= 0, ~\forall \tau_b \in \mathbb{V}_{b}, \\ \int_{\Omega} \boldsymbol{\chi_u} \cdot \frac{\partial{\mathbf{q_u}}}{\partial{t}} ~ \text{d} \boldsymbol{x} + \int_{\Omega} p^{'} {\boldsymbol{\nabla}{\cdot }}(\boldsymbol{\chi_u} \odot \boldsymbol{\gamma_V)} ~ \text{d} \boldsymbol{x} + \int_{\Omega} (\sigma + \alpha) \boldsymbol{\chi_u} \cdot \mathbf{q_u} ~ \text{d} \boldsymbol{x} &= 0, ~\forall \boldsymbol{\chi_u} \in \mathbb{V}_{\mathbf{u}}, \\ \int_{\Omega} \chi_p \frac{\partial{q_p}}{\partial{t}} ~ \text{d} \boldsymbol{x} - \int_{\Omega} c_{\text{s}}^2 \chi_p ({\boldsymbol{\nabla_V}}\cdot \mathbf{u}) ~ \text{d} \boldsymbol{x} + \int_{\Omega} (\sigma + \alpha) \chi_p q_p ~ \text{d} \boldsymbol{x} &= 0, ~\forall \chi_p \in \mathbb{V}_{p}, \end{align} \tag{32}\]

with \(\Gamma\) denoting the set of interior facets, and \(^\dagger\) denoting the upwind value on a facet, e.g. \[\label{def:upwind} \mathbf{q}^\dagger := \left\lbrace \begin{matrix} \mathbf{q}^+ & \mathrm{if} \;\mathbf{u} \cdot \hat{\mathbf{n}}^+ \geq 0, \\ \mathbf{q}^- & \mathrm{if} \;\mathbf{u} \cdot \hat{\mathbf{n}}^+ < 0, \end{matrix}\right.\tag{33}\]

where the ‘+’ and ‘-’ superscripts refer to arbitrarily labelled sides of the facet, \(\hat{\mathbf{n}}\) are outwardly facing normals, and \(\mathbf{q}\) could be a scalar or vector. The double square brackets denote the jump of a field across a facet, \[[[\mathbf{q}]] = \mathbf{q}^+ - \mathbf{q}^-.\]

3.3 Rayleigh damping sponge layer↩︎

To evaluate the PML for damping vertically propagating waves, we compare its performance to that of a Rayleigh damping sponge layer. Following [41], we introduce a sponge term of \(\mu \hat{\boldsymbol{\mathrm{k}}}(\mathbf{u} \cdot \hat{\boldsymbol{\mathrm{k}}})\). The corresponding weak form of \[\int_{\Omega} \mu(\boldsymbol{\tau_u} \cdot \hat{\boldsymbol{\mathrm{k}}})(\mathbf{u} \cdot \hat{\boldsymbol{\mathrm{k}}}) ~ \text{d} \boldsymbol{x}, \label{eq:sponge95term}\tag{34}\]

is added to the left-hand sides of the momentum equations, so (21 ) for the linear Boussinesq system and (31 ) for the nonlinear Boussinesq system. This Rayleigh damping term acts to damp the vertical velocity back to zero. Rayleigh damping can also be applied to horizontal velocities [2], [7], [42], [43], with damping of the vertical velocity being more common in nonhydrostatic models, such as MPAS [44], ICON [45], and GungHo [24].

The sponge layer damping function \(\mu(z)\) is defined analogously to \(\sigma(z)\) for the PML (18 ), with an inverse damping timescale of \(\mu_0\) used in place of \(\sigma_0\). For the acoustic tests, we examine sponge strengths for values of \(\mu_0\) separated by \(\Delta \mu_0 = 0.5 ~\text{s}^{-1}\). We note that more vertical levels can improve the efficacy of a sponge layer, but we use the same thickness as the PML for fair comparison, which is four elements in the vertical. The use of only a few models levels in the sponge layer is standard for dynamical cores, where it is important to minimise the additional computation, e.g. typically three layers are used in the finite volume cubed-sphere (FV3) model [42] and the Community Atmosphere Model (CAM) [46] version of the spectral element model [47], [48].

4 Numerical experiments↩︎

We now describe tests of the PML in the linear and nonlinear Boussinesq equations, to show the effective PML absorption properties for outgoing waves and verify numerical accuracy and stability. In particular, we compare the performance of the PML to the Rayleigh damping sponge layer that is the default for dynamical cores.

In these tests, we use a sound speed of \(c_{\text{s}} = 350 ~\text{m} ~\text{s}^{-1}\). The Brunt-Väisälä frequency is set to a constant \(N=0.01 ~\text{s}^{-1}\), which mimics an isothermal atmosphere. The domain is \([0,L_x] \times [0,H_z]\), with length \(L_x\) = 100 km and height \(H_z\) = 20 km. Next-to-lowest order elements of width \(\Delta x = 500 ~\text{m}\) and height \(\Delta z = 500 ~\text{m}\) are used.

4.1 Linear Boussinesq equations↩︎

We begin with the linear system, where we test the damping of vertically propagating acoustic waves.

4.1.1 Test setup↩︎

We consider an initial state that contains a perturbation to the pressure field. Often, the buoyancy field is instead perturbed in Boussinesq systems, most famously in the gravity wave test of [49]. However, the buoyancy perturbation excites predominantly horizontally propagating gravity waves that are much less prone to model top reflection in this configuration.

The pressure perturbation takes the form of a compactly supported Gaussian function at the centre of the domain, \[p(x,z,t=0)= \begin{cases} p_0 \exp\left( -\left( \frac{(x-L_x/2)^2+(z-H_z/2)^2}{d^2} \right) \right), &\text{if} ~x \in [0.3L_x, 0.7L_x] ~\text{and} ~z \in [0.3H_z, 0.7H_z], \\ 0, &\text{else}. \end{cases}\]

The maximum pressure perturbation magnitude is \(p_0 = 10 ~\text{Pa}\) and the Gaussian half-width is \(d = 1 ~\text{km}\). The initial velocity and buoyancy fields are zero. A flat surface at \(z_{\text{s}}=0\) is used. The acoustic waves take 28.6 s to travel \(10 ~\text{km}\), i.e. from the centre of the domain to the model top or surface. An explicit RK4 timestepper is used with \(\Delta t = 0.25 ~\text{s}\), which leads to an acoustic wave Courant number of \(\textrm{Cr}_{\text{ac}}= c_{\text{s}} \Delta t/\Delta x=0.175\). The simulation length is \(T_{\text{end}} = 100 ~\text{s}\).

We compute a reference solution for quantitative comparisons of no model top damping, a sponge layer, or a PML. The reference solution applies no model top damping, but uses a higher domain with \(H_z = 50 ~\text{km}\) to avoid wave reflections. Figure 2 shows the end state of the reference solution. There are two wave fronts, as the initial pressure perturbation excites upwardly and downwardly travelling components. The downwardly travelling component reflects off the bottom surface and then travels upwards, thus ending at a lower height than the upwardly travelling crest. In the next section, errors relative to this reference solution are computed as the total \(L_2\) error of the vertical velocity within the undamped domain of \(z \in [0,18] ~\text{km}\).

4.1.2 Results↩︎

The vertical velocity and pressure fields from simulations with no damping, a sponge layer, and the PML with \(\gamma_0=0.25 ~\text{s}\) are compared with the reference solution in Figure 3. The best sponge strength leading to the smallest vertical velocity error is identified as \(\mu_0 = 1.5 ~\text{s}^{-1}\), which is weaker than the PML coefficient of \(\sigma_0 = 2.42 ~\text{s}^{-1}\) (19 ). In Figure 3, the undamped solution has clear model top reflections from both the upwardly and downwardly propagating parts of the acoustic wave, the latter after it has reflected off the bottom surface. With a sponge layer, the reflection is weaker and the amplitudes of \(w\) are reduced, but there is still a noticeable spurious impact on the dynamics. With the PML, both the upwardly and downwardly propagating acoustic crests are successfully damped once they enter the PML, leaving no visual evidence of reflection. Accordingly, the solution is visually indistinguishable from the reference solution in the undamped domain.

The improvement with the PML is quantified through a time series of the vertical velocity error in Figure 4. Both the unstretched (\(\gamma_0=0 ~\text{s}\)) and stretched (\(\gamma_0 \neq 0~\text{s}\)) PML solutions are two orders of magnitude more accurate than solutions with no damping or a sponge layer. The use of a grid stretching factor of \(\gamma_0=0.25 ~\text{s}\) further reduces the error compared to the unstretched PML when \(t > 75 ~\text{s}\); this stretching factor was identified as the choice with lowest vertical velocity error from \(\gamma_0 \in \{0, 0.05, 0.1, 0.25, 0.5 \} ~\text{s}\). Hence, the use of real grid stretching can improve the damping ability of the PML, although making this too large begins to degrade the accuracy.

We now examine the energetics of the linear Boussinesq system (25 ) in Figure 5,,. The initial condition only contains internal energy (due to the linear equation set solving for perturbations around a stably stratified state) which is primarily exchanged with kinetic energy during the simulation. In the undamped case, there are three bumps in the kinetic and internal energies after the initial disturbance, which are due to i) the upwardly travelling crest of the acoustic wave hitting the model top and the downwardly travelling crest hitting the bottom surface (\(t \approx 28.6\) s), ii) the upwardly and downwardly travelling crests meeting in the middle of the domain (\(t \approx 57.1\) s), and iii), the initially downwardly travelling wave hitting the model top (\(t \approx 85.7\) s).

The undamped solution maintains a fairly constant total energy throughout the simulation, which is a consequence of a conserved energy in the linear Boussinesq equations (25 ). With a sponge layer, there is a reduction in total energy once the upwardly travelling crest enters the damping region. The second and third bumps in \(E_K\) and \(E_I\) are still present with a sponge, although with a smaller amplitude than in the undamped case. With the PML, there is again a reduction in total energy as the upwardly travelling crest enters the damping region. A faster energy drop off with the PML compared to the sponge indicates a more effective removal of acoustic wave energy. The first bump in the energies is present, as there is still a reflection of the downwardly travelling crest at the bottom surface, but there are no second and third bumps, indicating the absence of wave reflection at the model top.

Lastly, Figure 6 shows the rate of energy removal in the PML for different grid stretching factors. The PML power (30 ) is zero until the upwardly travelling crest enters the PML region and is damped. The PML power peaks when damping this first acoustic wave. There is also a second region of strong PML power when the downwardly travelling crest reaches the PML region after its reflection off the bottom surface. The most energy removed is with \(\gamma_0=0.1 ~\text{s}\), although \(\gamma_0 = 0.25 ~\text{s}\) leads to the smallest vertical velocity error in this test.

4.2 Nonlinear Boussinesq equations↩︎

We now move to the nonlinear Boussinesq equations, where we test the damping of both acoustic and orographic gravity waves.

4.2.1 Test setup↩︎

We apply an analytical initial condition of \[\begin{align} \mathbf{u}(x,z,t=0)=\overline{\mathbf{u}} &= [10,0]^T, \\ b(x,z,t=0)=\overline{b}(z) &= N^2 z, \\ p(x,z,t=0)=\overline{p}(z) &= \frac{N^2 z^2}{2} - p_0, \end{align}\]

with \(p_0 = 200 ~\text{hPa}\). This initial condition satisfies the Boussinesq hydrostatic balance of \({\boldsymbol{\nabla}{p}}=b\mathbf{\hat{k}}\). Alternatively, one could enforce hydrostatic balance numerically with the method of [50], which iteratively computes a \(\overline{p}\) that is in discrete hydrostatic balance with \(\overline{b}\) in each column to some error tolerance. We choose an analytical expression as this enables the computation of a reference solution with a higher model top, whereas the discrete hydrostatic balance method leads to a different \(\overline{p}(z)\) with higher columns.

We use a varied surface height to enable orographic gravity waves in this test. We define a surface height containing a Gaussian mountain, \[z_{\text{s}} = h_0 \exp \left( -\left( \frac{x-x_{\text{c}}}{a} \right)^2 \right),\]

with a mountain height of \(h_0 = 1 ~\text{km}\), and Gaussian half-width of \(a = 10 ~\text{km}\). The mountain is centred at the middle of the domain, \(x_{\text{c}}=L_x/2\). The mountain is sufficiently high to enter the nonlinear ‘flow-around’ regime, which leads to more complicated dynamics than the linear ‘flow-over’ regime of small mountains [51], [52].

We use a hybrid terrain-following vertical coordinate of \[z = \overline{z} + \frac{H_z-\overline{z}}{H_z} z_{\text{s}}, \label{eq:hybrid95z95coord}\tag{35}\]

with \(\overline{z}\) a reference vertical coordinate with constant height levels evenly spaced by \(\Delta \overline{z} = 500 ~\text{m}\). This choice ensures that \(z=z_{\text{s}}\) at the surface and the model top at \(z=H_z\) is flat.

The same initial condition is used for the acoustic and orographic gravity wave configurations of this test. The acoustic wave is generated through an incompatibility of the horizontal flow in the initial condition with the impermeability boundary condition of \(\mathbf{u} \cdot \mathbf{\hat{n}}=0\) at the bottom surface. Unlike the linear Boussinesq test, the acoustic wave starts at the surface and only contains an upwardly travelling crest. The orographic gravity wave is generated from the horizontal flow over the mountain and propagates on a much slower timescale than the acoustic wave. To specifically test the damping of each type of wave, we vary the time discretisation and simulation length:

• For the acoustic test: Explicit RK4 timestepping, \(\Delta t = 0.25 ~\text{s}, ~T_{\text{end}} = 150 ~\text{s}\). The acoustic wave Courant number is \(\textrm{Cr}_{\text{ac}}=0.175\).

• For the orographic gravity wave test: Implicit midpoint timestepping (also called the Crank-Nicolson scheme or trapezoidal rule), \(\Delta t = 10 ~\text{s}, ~T_{\text{end}} = 10,000 ~\text{s}\). The acoustic Courant number is \(\textrm{Cr}_{\text{ac}}=7\) and advective Courant number is \(\textrm{Cr}_{\text{adv}}=u_0 \Delta t/\Delta x = 0.2\), with \(u_0 = 10 ~\text{m} ~\text{s}^{-1}\).

4.2.2 Acoustic test results↩︎

For the acoustic test, we again compute a reference solution with a higher model top of 50 km and no damping layer. The reference solution uses the terrain following vertical coordinate (35 ) for \(z \in [0, 20] ~\text{km}\) with equally spaced levels of element height \(\Delta z=500 ~\text{m}\) for \(z\in [20,50] ~\text{km}\). With the simulation time of \(T_{\text{end}} = 150 ~\text{s}\), the acoustic waves travel 52.5 km and do not interfere with the solution in the undamped domain of \(z \in [0, 18] ~\text{km}\).

Figure 7 shows vertical velocity fields at \(t=100 ~\text{s}\) and the final time of \(t=150 ~\text{s}\) in the acoustic test. Without damping, there is clear evidence of acoustic wave reflection: at \(t=100 ~\text{s}\) the wave has reflected off the model top and interferes with the flow near the mountain, whilst at \(t=150 ~\text{s}\) the wave has reflected off the bottom surface and is travelling upwards again. The use of a sponge layer of \(\mu_0=3.5 ~\text{s}^{-1}\), which is the choice of damping timescale with the lowest vertical velocity error in this test, damps the acoustic wave to some degree, but there remain spurious features in the solution from wave reflection. With the PML, the acoustic waves are successfully damped, and the flow appears very close to the reference solution.

A time series of the \(L_2\) vertical velocity error is shown in Figure 8. There is an order of magnitude accuracy improvement with a PML compared to the sponge and undamped solutions. Like with the linear Boussinesq test, a nonzero vertical grid stretching factor reduces the error, although there is minimal difference between \(\gamma_0=0 ~\text{s}\) and \(\gamma_0=0.25 ~\text{s}\) until \(t > 110 ~\text{s}\). Grid stretching factors of \(\gamma_0 \in \{0.25, 0.5, 1\} ~\text{s}\) lead to similar solutions.

4.2.3 Orographic gravity wave test results↩︎

We now move to the longer simulation that tests the PML’s ability to damp orographic gravity waves. An approximate reference solution is computed using a 50 km model top and a 30 km sponge layer, using a strength of \(\mu_0 = 3.5 ~\text{s}^{-1}\). The increased sponge thickness leads to smaller gradients in the damping profile (18 ) compared to a thickness of \(\delta = 2 ~\text{km}\). During the simulation length, the orographic gravity wave will not travel far enough to reflect off the high model top, and the large sponge should sufficiently damp the acoustic waves.

Figure 9 compares vertical velocity fields at the end time of 10,000 s. The undamped solution contains a ‘checkerboard’ standing wave pattern due to reflection at the model top, which degrades the representation of the orographic gravity wave throughout the domain. When using a sponge layer, there is a reduction of this checkerboarding, but it is still present. A weaker sponge coefficient is needed compared to the acoustic test, as \(\mu_0 = 3.5 ~\text{s}^{-1}\) leads to strong reflections at the onset of the damping layer at \(z=18 ~\text{km}\). This highlights a drawback of sponge layers, in that they are difficult to tune to damp both acoustic and gravity waves; this means that sponge layers are typically only designed to damp gravity waves. We identify \(\mu_0 = 0.1 ~\text{s}^{-1}\) as a suitable choice for damping the orographic gravity wave in this test.

Simulations with a PML in the orographic gravity wave test are visually free of the checkerboard reflection pattern, with a clearer representation of the orographic gravity wave signal for \(z \leq 18 ~\text{km}\). However, the use of different grid stretching factors \(\gamma_0 \in \{0, 0.25, 0.5, 1 \}\) s modifies the representation of the orographic gravity wave; each choice of \(\gamma_0\) replicates its general structure, but there are notable discrepancies relative to the reference solution. Hence, the currently formulated PML successfully avoids the checkerboarding reflection but still needs improvement for a more accurate representation of orographic gravity waves. This will likely require additional terms in the PML equations pertaining to advection or buoyancy.

This test still excites an acoustic wave, albeit one that is not properly resolved, due to the implicit timestepping method and a long timestep size of acoustic Courant number \(\textrm{Cr}_{\text{ac}}=7\). To investigate the impact of the acoustic waves, Figure 10a plots time series of the integrated vertical kinetic energy, \(E_{\text{K,V}} = \int_{\Omega} 0.5w^2 ~ \text{d} \boldsymbol{x}\). In the undamped and sponge solutions, there are high frequency oscillations in the vertical kinetic energy corresponding to the acoustic wave bouncing between the model top and bottom surface. Such oscillations in the vertical kinetic energy have also been observed for energy-conserving discretisations of the compressible Euler equations on the Lorenz grid [37], [38]. The use of the PML leads to a smoother time evolution of \(E_{\text{K,V}}\), indicating an absence of this trapped acoustic wave. Hence, the PML can damp the vertically propagating waves even when they are not resolved by the timestepping method. A plot of spectral energy (Figure 10b), computed with the temporal Fourier transform of \(E_{\text{K,V}}\), shows that there are high frequency components in the undamped and sponge solutions that are removed with the PML.

We lastly examine the impact of the PML on the energy in the nonlinear Boussinesq equations. Whilst we have not identified a conserved energy for the nonlinear system, we nevertheless use the same PML power term as for the linear system (30 ) as an approximate measure of energy removal from the PML. We note that for all cases, irrespective of damping choice, there is energy dissipation in the nonlinear equations from upwinding of the advective terms (31 ).

Figure 11 examines the net energy removal from the PML by plotting the time-integrated PML power, approximated as \[\Delta E_{\text{PML}}(t_n) = \int_{t=0}^{t_n} ~P_{\text{PML}} ~\textrm{d} t \approx \sum_{i=1}^{n} \Delta t ~P_{\text{PML}}(t_i),\]

where \(t_i\) are the times used in the discrete solution. Figure 11 shows that a larger grid stretching factor \(\gamma_0\) leads to more energy removal during the simulation; the difference is clearest at the later times, when the orographic gravity wave has reached the PML. Increasing the complex frequency shift factor, \(\alpha\), reduces the total energy removed.

5 Conclusions and discussion↩︎

This paper introduced a novel application of the perfectly matched layer for damping vertically propagating waves in atmospheric models. In the continuous equations, the PML acts to damp all waves that enter the damping region, irrespective of angle of incidence and wavelength, and avoids reflection at the onset of the damping layer. This study derived PML equations for the compressible Boussinesq equations, which are a valuable model of vertical dynamics in the atmosphere and support acoustic and gravity waves.

Numerical tests of the PML equations were performed with a compatible finite element model. We set the PML to cover four vertical levels of next-to-lowest-order elements, so although auxiliary PML variables need to be introduced, the increase in computational cost is somewhat mitigated by the thin damping layer. The proposed PML damps motions created by the pressure gradient and divergence terms in the Boussinesq system. Accordingly, this requires velocity and pressure PML variables but not a buoyancy PML variable.

Tests with the linear Boussinesq equations showed that the PML is more effective at damping acoustic waves than a standard vertical velocity sponge layer. We then provided a novel use of linear PML theory to the nonlinear Boussinesq equations, which damps perturbations from a hydrostatically balanced state. Numerical tests in the nonlinear system showed that the PML can avoid model top reflections from both acoustic and orographic gravity waves. This is in contrast to the sponge, which still permitted model top reflections in the acoustic test with a strong coefficient, whilst needing a much weaker coefficient for the orographic gravity wave test.

The use of a real grid stretching factor in the PML improved the damping of acoustic waves for both the linear and nonlinear Boussinesq systems. However, different stretch factors led to modified representations of the orographic gravity wave in the nonlinear test, although all choices successfully avoided model top reflection. Hence, an important next step is improving our PML for orographic gravity wave applications. This could be achieved by including additional terms in the PML, such as a measure of the buoyancy oscillation frequency or a linearisation of the advective terms, possibly about a nonzero mean flow. Additionally, improvements may also be obtained through different selections of PML parameters.

There are a number of future research topics following this study. First is the development of stability proofs for our Boussinesq PML formulation, using error estimates within the Laplace domain. Second, we aim to derive a PML for the compressible Euler equations in the vertical slice domain. The Euler equations are more complex than the Boussinesq equations, due to a nonlinear pressure gradient term and a nonconstant acoustic wave speed. For the compressible Euler PML, we would like to incorporate a linearisation of internal gravity wave propagation. Third, more sophisticated timestepping methods could be used with the PML equations. Specifically, the PML could be incorporated into the semi-implicit quasi-Newton method that is often used for the compressible Euler equations with compatible finite element models, e.g. [53]. Last, we eventually aim to implement the PML in a fully three-dimensional model on the sphere, thus testing the potential for a PML to replace Rayleigh damping or Laplacian sponge layers in dynamical cores.

Code Availability↩︎

The Gusto code base is available at the GitHub repository https://github.com/firedrakeproject/gusto. Python scripts used to run the simulations and create the figures are available at the GitHub repository https://github.com/ta440/vertical_slice_PML. Functions from the tomplot Python code library were used for the figures, and these can be found at the GitHub repository https://github.com/tommbendall/tomplot.

Acknowledgments↩︎

TCA thanks Christiane Jablonowski for the discussions about Rayleigh damping sponge layers, which led him down the path of exploring other options like the PML.

References↩︎