Performance of Adaptive Unstructured Mesh Modelling in Idealized Advection Cases over Steep Terrains

Abstract

Advection errors are common in basic terrain-following (TF) coordinates. Numerous methods, including the hybrid TF coordinate and smoothing vertical layers, have been proposed to reduce the advection errors. Advection errors are affected by the directions of velocity fields and the complexity of the terrain. In this study, an unstructured adaptive mesh together with the discontinuous Galerkin finite element method is employed to reduce advection errors over steep terrains. To test the capability of adaptive meshes, five two-dimensional (2D) idealized tests are conducted. Then, the results of adaptive meshes are compared with those of cut-cell and TF meshes. The results show that using adaptive meshes reduces the advection errors by one to two orders of magnitude compared to the cut-cell and TF meshes regardless of variations in velocity directions or terrain complexity. Furthermore, adaptive meshes can reduce the advection errors when the tracer moves tangentially along the terrain surface and allows the terrain to be represented without incurring in severe dispersion. Finally, the computational cost is analyzed. To achieve a given tagging criterion level, the adaptive mesh requires fewer nodes, smaller minimum mesh sizes, less runtime and lower proportion between the node numbers used for resolving the tracer and each wavelength than cut-cell and TF meshes, thus reducing the computational costs.

1. Introduction

The basic terrain-following (BTF) coordinate has been widely employed in numerical weather prediction (NWP) models due to its simplicity and ability to accurately represent the terrain [1,2,3,4,5,6,7]. However, advection errors can occur in a BTF coordinate system when the velocity field is not aligned with terrain-following (TF) vertical layers (VLs). With increasingly improved horizontal model resolutions enabled by advances in computing technology, the slopes of the terrains represented in NWP models are becoming progressively steeper. This resolution improvement is making it increasingly difficult to accurately represent advection processes over steep terrains. This difficulty may even cause numerical instabilities [8]. To reduce this source of advection errors, model developers have proposed the use of either special vertical coordinates or computational grid/mesh approaches to address this problem.

Regarding the use of special vertical coordinates, two main categories of methods are available: hybrid terrain-following (HTF) coordinates [9,10] and smoothed VLs based on BTF coordinates. In an HTF coordinate system, the BTF coordinate is used in the vicinity of the terrain, and a height (or pressure) coordinate is used in the upper levels of the atmosphere, thereby smoothing the VLs and coupling them to the upper-level atmospheric flow. Alternatively, smoothed VLs can be achieved by constructing mathematical functions with the aim of reducing both the misalignment of velocity fields with VLs and the grid distortions. Many typical methods are utilized for this purpose, including the smooth level vertical (SLEVE) coordinate [4], the smoothed terrain-following (STF) coordinate [11], the cosine-smoothed terrain-following (COS) coordinate [12], and the orthogonal terrain-following (OTF) coordinate [2].

Regarding the use of special computational grid/mesh approaches, alternative non-TF coordinates with local boundary mesh schemes have been developed to improve the interaction of mountains with the flow. These can broadly be divided into conformal meshes, which follow the boundary contours and embedded boundary methods, in which the presence of the boundary is modelled in a non-conformal mesh by adding a distribution of source terms [13,14,15]. The cut-cell method [5,16,17,18], based on a height coordinate, is another approach for representing terrains by approximately inserting the terrain into the domain through a bilinear spline or a quadratic function. The computational cells are cut by the terrain. The cut-cell method has been shown to effectively reduce upper-level advection errors and spurious winds in the TF coordinate [6,17,19]. Both gentle and steep 2D terrains have been utilized to verify the accuracy and stability of the cut-cell method [16,18,19,20,21]. However, the small cut-cells around the terrain impose constraints on the time step, introducing issues related to stability and computational efficiency, which can be alleviated by a thin-wall approximation [18] and a cell-merging technique [19]. Moreover, several 3D applications of cut-cells in atmospheric modelling studies have been undertaken [22,23,24,25,26].

Substantial improvements in the alignment of the velocity field with the VLs and in reducing advection errors have been achieved by modifying TF coordinates and generating new non-TF coordinates. Nevertheless, in conformal mesh schemes, the bottom-most grid layer is always forced to follow the underlying terrain either approximately or accurately, which results in misalignments between the mesh face normals and the flow direction at the peaks of terrain profiles, thus inevitably inducing numerical errors. Moreover, Shaw and Weller [5] confirmed that the advection accuracy depends on the alignment of the velocity field with the VLs, indicating that a TF coordinate should be able to effectively simulate TF advection processes, while the cut-cell method is preferred for horizontal advection problems. Finally, Zou et al. [27] demonstrated that advection errors are affected by the complexity of the terrain; that is, for steeper terrains, TF coordinates will induce larger advection errors. Consequently, it may be difficult to find a coordinate that is suitable for all types of meteorological flows [6]. It is also noted that with regard to the cut-cell method, only a few studies have focused on the advection along the surfaces of complex terrains, with the exception of [6,28], who utilized a constant slope terrain to test advection errors and the effect of grid stagger.

From a computational mesh perspective, another alternative for simulating atmospheric flows is proposed using unstructured adaptive meshes based on the discontinuous Galerkin (DG) finite element method (FEM) [29,30,31,32,33,34,35,36]. The common adaptive meshes are the nested adaptive, stretched adaptive and dynamically adaptive mesh [32]. The nested adaptive mesh has found widespread use in many models, for example, in Weather Research and Forecasting Model. However, the interaction between the high- and low-resolution meshes is still a concern owing to the mismatch across their boundary. Iselin [37] utilized the stretched adaptive mesh to deal with 1D and 2D advection problems. Bacon et al. [29] used horizontal dynamically adaptive mesh to test the advection process in an Operational Multiscale Environment Model with Grid Adaptivity. Weller et al. [38] used the optimal transport and numerical solution of a Monge–Ampère type equation to adapt meshes. The dynamically adaptive mesh technique was introduced by [36] to accurately represent air pollutant transport processes. In this study, the dynamically adaptive mesh and a corresponding optimization-based adaptive technique are used because the computational meshes can be adapted (in both time and space) with specified variables to resolve small-scale flow features as the flow evolves. Such an adaptive mesh provides several advantages [33,39,40,41]. First, the effective adaptation of the mesh allows a good representation of complex topographic and flow features, thereby eliminating spurious oscillations and reducing spatial and temporal interpolation errors. Second, the adaptive mesh technique allows the meshes to automatically become denser where the gradients of variable solutions are higher and to remain relatively sparse in other parts of the domain. Thus, an acceptable global accuracy is achieved while saving the computational resources. Third, the meshes can be temporally adjusted and are therefore able to both resolve multiscale processes and capture flow features over time.

In this study, adaptive unstructured meshes based on the DG-FEM are introduced to investigate their effects on advection simulations. The adaptive modelling framework known as Fluidity and its adaptive mesh technique are described in Section 2. Section 3 presents the governing equations and parameters employed in all subsequent advection experiments. Section 4 presents a series of 2D idealized advection test cases to investigate the effects of the adaptive meshes. Classic horizontal and TF advection tests are used to test the capabilities of the adaptive meshes in reducing advection errors. Then, variations in the velocity direction and increasing terrain complexity are considered to investigate whether the simulation accuracy of the adaptive meshes is sensitive to variations in the velocity directions and terrains. Moreover, the tracer for the advection test is positioned next to the surface of the terrains to test the effect of the lower boundary in the adaptive meshes. In addition, a short comparison of the computational cost is given in Section 4.2. Finally, Section 5 presents a discussion and concludes this study.

2. Adaptive Model: Fluidity

This study employs an unstructured adaptive mesh scheme known as Fluidity developed by the Applied Modelling and Computation Group at Imperial College London. Fluidity supports many advanced numerical techniques, including numerical stabilization schemes, anisotropic fully unstructured meshes and conservative adaptive mesh optimization schemes, in addition to new mixed finite elements that exactly represent key balances within the governing equations. In this study, the DG-FEM is used for spatial discretization which is well suited for h-r adaptivity and parallelization problems (see details in [36]).

This study introduces a type of dynamically adaptive unstructured mesh [41] to reduce advection errors and both horizontal and vertical anisotropic unstructured meshes are used to capture the details of flows in all three directions [33,34,41]. Such dynamically adaptive meshes have numerous advantages over other widely used meshes (e.g., the cut-cell and TF meshes depicted in Figure 1a,b) and existing adaptive mesh refinement (AMR) methods (e.g., locally nested dynamic mesh methods); for example, a dynamically adaptive mesh can accurately represent the boundary conditions over complex bathymetries (i.e., steep terrains). Moreover, anisotropic adaptive mesh technology enables such a mesh to be aligned with the bathymetry and additionally adapts the mesh to advection flows.

The refinement tagging criterion is defined by

$$\widehat{\epsilon }={\upsilon }^T\left|H\right|\upsilon$$

(1)

where $H=\left(\begin{array}{cc}{\partial }^{2}q/\partial x^2& {\partial }^{2}q/\partial x\partial y\\ {\partial }^{2}q/\partial x\partial y& {\partial }^{2}q/\partial y^2\end{array}\right)$ is the Hessian matrix for a specified field (in this study, the field of the tracer concentration q) and the vector $\upsilon$ is the dimension of a given mesh element. The relationship between the tagging criterion and the size of the element can be clearly seen in 1D. The tagging criterion in 1D is ${\widehat{\epsilon }}_{1D}=h^2\left|{\partial }^{2}q/\partial x^2\right|$ where h is the length of the 1D element. Therefore, the tagging criterion is directly proportional to the size of the element.

With this tagging criterion, a metric tensor

$$\widehat{M}=\frac{1}{\widehat{\epsilon }}H$$

(2)

is used to calculate the required edge lengths and the orientation of the mesh elements. Then $\widehat{M}$ is further modified by considering the maximum and minimum element sizes:

$$M=V^T\Lambda V$$

(3)

where V is the right eigenvectors matrix of the metric $\widehat{M}$. The eigenvalues ${\lambda }_j,\forall j\in \left\{1,2\right\}$ of $M$ are modified as follows:

$${\lambda }_j=\max \left\{{\lambda }_j^{\prime },\frac{1}{a^2}\underset{i=1,2}{\max }{\lambda }_j^{\prime }\right\},\forall j\in \left\{1,2\right\}$$

(4)

with

$${\lambda }_j^{\prime }=\min \left\{\frac{1}{h_{\min }^2},\max \left\{\left|{\widehat{\lambda }}_j\right|,\frac{1}{h_{\max }^2}\right\}\right\},\forall j\in \left\{1,2\right\}$$

(5)

where a is the maximum aspect ratio, h_min and h_max are the minimum and maximum element sizes and ${\widehat{\lambda }}_j,\forall j\in \left\{1,2\right\}$ are the eigenvalues of metric $\widehat{M}$ [41]. In this study, a range of thresholds for the tagging criterion in the advection tests is given to demonstrate the capability of the adaptive mesh.

Local mesh operations are based on a series of mesh connectivity and node position searches (e.g., edge collapsing/splitting, face-to-edge and edge-to-face swapping, edge-to-edge swapping and local node movement or mesh smoothing). The quality of adaptive meshes can be defined by imposing suitable tolerances on the interpolation errors as well as minimum and maximum element sizes and aspect ratios and by establishing a cut-off point for the magnitudes of state variables below which a fine mesh resolution is deemed unnecessary. The Galerkin project of discrete fields between two meshes is implemented via a supermeshing algorithm [42], which ensures the accuracy and the conservation of the solutions.

3. Governing Equations and Parameters

3.1. Governing Equation and Discretization

In this study, the following linear advection equation is used in all the tests:

$$\frac{\partial q}{\partial t}+\mathbf{u}\cdot \nabla q=0$$

(6)

where q represents a scalar tracer concentration and $\mathbf{u}={\left(u,w\right)}^T$ is the velocity vector.

The DG-FEM and a Cartesian height coordinate are utilized to obtain the discretized form of Equation (6). In the DG-FEM, the variable q can be expressed as follows:

$$q\left(X\right)={\displaystyle \sum _{i=1}^N{\varphi }_i\left(X\right)q_i}$$

(7)

where φ_i and q_i are the values of the test function φ and the variable q, respectively, at node i; N is the number of nodes; and X is the location of a node. By multiplying Equation (6) by the test function φ_i and integrating over the domain Ω, the integral form of Equation (6) is obtained:

$${\displaystyle {\int }_{\Omega }\left(\varphi \frac{\partial q}{\partial t}+\mathbf{u}\cdot \nabla \left(\varphi q\right)\right)d\Omega }=0$$

(8)

The semi-discrete matrix form of Equation (8) can be written as follows:

$$\mathbf{M}\frac{\partial \mathbf{q}}{\partial t}+\mathbf{A}\left(\mathbf{u}\right)\mathbf{q}=\mathbf{r}$$

(9)

where the vector $\mathbf{q}={\left(q_1,q_2,\dots ,q_N\right)}^T$ contains the solutions for the variable q at various nodes, M is the mass matrix, A(u) is the advection operator with first-order upwind scheme, and r is the right-hand side vector containing the boundary, source and absorption terms, with

$${\mathbf{M}}_{ij}={\displaystyle {\int }_{\Omega }{\varphi }_i{\varphi }_{j}d\Omega },{\mathbf{A}}_{ij}=-{\displaystyle {\int }_{\Omega }\mathbf{u}\cdot \nabla \left({\varphi }_i{\varphi }_j\right)d\Omega },i,j\in \left(1,2,\dots ,N\right)$$

(10)

In this study, the test functions φ_i (i = 1, 2, …, N) are first-order linear polynomials (see details in [36]).

The time derivative term at time level n + 1 is treated using the θ-method [43,44], yielding

$$\mathbf{M}\frac{{\mathbf{q}}^{n+1}-{\mathbf{q}}^n}{\mathsf{\Delta }t}+\mathbf{A}\left({\mathbf{u}}^n\right){\mathbf{q}}^{n+\theta }={\mathbf{r}}^{n+\theta }$$

(11)

where Δt is the time-step, $\theta \in [0,1]$ and the vector ${\mathbf{q}}^{n+\theta }$ is given by

$${\mathbf{q}}^{n+\theta }=\theta {\mathbf{q}}^{n+1}+\left(1-\theta \right){\mathbf{q}}^n$$

(12)

Equation (11) can be rearranged to solve for the unknown vector ${\mathbf{q}}^{n+\frac{1}{2}}$:

$$\left[\mathbf{M}+\frac{1}{2}\mathsf{\Delta }t\mathbf{A}\left({\mathbf{u}}^n\right)\right]{\mathbf{q}}^{n+1}=\left[\mathbf{M}-\frac{1}{2}\mathsf{\Delta }t\mathbf{A}\left({\mathbf{u}}^n\right)\right]{\mathbf{q}}^n+{\mathbf{r}}^{n+\frac{1}{2}}$$

(13)

In this study, the semi-implicit time-stepping scheme and the upwind spatial discretization scheme are utilized in the DG-FEM method where θ = 0.5 in Equation (13). The free lateral boundary condition is applied to the horizontal boundaries, the non-slip boundary condition is applied to the bottom boundary, and the free-slip condition is applied to the top boundary [45,46]. The Generalized Minimum Residual Method is used to solve the large sparse linear systems of discrete advection equation. For the mesh adaptivity, the mesh is adapted at every two time-steps for computational efficiency, and the adaptation process is based on the mesh adaptivity technique of [41] using the criterion $\widehat{\epsilon }$ greater than O(10⁻³).

3.2. Parameters and Meshes

The computational domain is 300 km wide in the horizontal direction and 25 km deep in the vertical direction. The horizontal and vertical resolutions in the computational space are dX = 1 km and dZ = 0.5 km for regular grids and vary from 0.1 km to 10 km for adaptive mesh modelling. Owing to the semi-implicit time-stepping scheme, the Courant number is limited to 7. Therefore, the modelling period is 10,000 s, with an integration time step of 25 s, which is stable for the simulations.

At the start of the computation, the computational domain interior is initialized by imposing the tracer concentration distribution

$$q\left(x,z\right)=\{\begin{array}{ll}{q}_0{\cos }^2\left(\frac{\pi }{2}r\right),& r\le 1\\ 0,& r>1\end{array}$$

(14)

where q₀ = 1.0, $r=\sqrt{{\left(\frac{x-x_0}{A_x}\right)}^2+{\left(\frac{z-z_0}{A_z}\right)}^2}$, A_x = 25 km is the horizontal radius, A_z = 3 km is the vertical radius, and the tracer is centered at x₀ = 100 km and z₀ = 7 km. The tracer is uniform in the y-direction. For comparative purposes, three types of meshes, namely, cut-cell meshes, TF meshes, and the proposed adaptive unstructured meshes based on a height coordinate, are used in all numerical tests. Examples of adaptive unstructured meshes are shown in Figure 1c–f (the difference in the mesh resolutions is the thresholds of the criterion $\widehat{\epsilon }$ where $\widehat{\epsilon }$ is O(10⁻²) in Figure 1c and O(10⁻³) in Figure 1d). The distribution of the absolute values of the absolute error (E_a) and l₂ error norms on all the computational meshes are defined as follows:

$$E_a=\left|q_i\left(x,z\right)-q_a\right|$$

(15)

$$l_2=\sqrt{\frac{{\displaystyle \sum _i{\left[q_i\left(x,z\right)-q_a\right]}^2{V}_i}}{{\displaystyle \sum _i{q}_a^2{V}_i}}}$$

(16)

where i is the node index, q_i(x,z) and q_a are the corresponding numerical and analytical tracer concentration distribution [4,5], and V_i is the cell volume.

4. Results

To investigate the advection performance achieved using adaptive meshes in 2D numerical simulations, five advection test cases are implemented with different types of velocity fields and underlying terrains. The results and corresponding analysis are shown in Section 4.1, and a short comparison of the computational cost is discussed in Section 4.2.

4.1. 2D Advection Tests

In this section, the following five advection tests are investigated:

• Case I: a horizontal advection test similar to [4] which is not meant to represent a real atmospheric boundary layer scenario but just aims to evaluate the scalar advection performance in a simple analytical velocity field;
• Case II: a TF advection test similar to Case I except with the TF velocity field as proposed by [5];
• Case III: advection tests similar to Case I except with spatially non-uniform velocity fields, including an upward linear (III-1), a downward linear (III-2) and a wavelike (III-3) velocity field. In these ascending and descending velocity test scenarios, the interaction between the terrain and the imposed velocity field is neglected;
• Case IV: advection tests similar to Case I except for the terrains, including a one-crest (IV-1), a five-crest (IV-2) and a ten-crest (IV-3) terrain;
• Case V: an advection test similar to Case II except that the tracer moves along the terrain surface.

The results of these tests and the corresponding illustrations are provided in Section 4.1.1–Section 4.1.5, respectively.

4.1.1. Case I: Horizontal Advection Test

Following [4], a horizontal advection test is employed over a Schär-type five-crest wavelike terrain (Figure 2a). The five-crest wavelike terrain is defined as follows:

$$h_5\left(x\right)={\cos }^2\left[\frac{\pi }{\lambda }\left(x-150\right)\right]h^{*}\left(x\right)$$

(17)

where

$$h^{*}\left(x\right)=\{\begin{array}{ll}{h}_0{\cos }^2\left[\frac{\pi }{2a}\left(x-150\right)\right],& \left|x-150\right|\le a\\ 0,& \left|x-150\right|>a\end{array}$$

(18)

where h₀ = 3 km is the maximum height of the terrain, a = 25 km is the half-width of the terrain, and λ = 8 km.

The constant 2D velocity field u(x,z) has Cartesian velocity components

$$u\left(z\right)=u_0\{\begin{array}{cc}1,\hfill & \hfill z_2\le z\\ {\sin }^2\left(\frac{\pi }{2}\frac{z-z_1}{z_2-z_1}\right),\hfill & \hfill z_1\le z\le z_2\\ 0,\hfill & \hfill z\le z_1\end{array}$$

(19)

$$w\left(x,z\right)=0$$

(20)

where u₀ = 10 m/s, z₁ = 3 km and z₂ = 4 km. The results are compared using cut-cell, TF and adaptive meshes at the beginning (t = 0 s), middle (t = 5000 s) and end (t = 10,000 s) of the integration period (Figure 3). The E_a distributions and l₂ error norms of each mesh are shown in Figure 4 and listed in

Table 1. l₂ error norms at t = 10,000 s from the five groups of 2D advection tests presented in Section 4.1.1–Section 4.1.5.

Test Cases		Cut-Cell Mesh		TF Mesh		Adaptive Mesh
		Middle	End	Middle	End	Middle	End
Case I		3.06 × 10−2	3.96 × 10−2	1.09 × 10−1	1.69 × 10−1	1.03 × 10−3	1.46 × 10−3
Case II		1.12 × 10−1	1.51 × 10−1	3.82 × 10−2	4.50 × 10−2	3.49 × 10−3	3.07 × 10−3
Case III	III-1	1.56 × 10−2	1.78 × 10−2	5.62 × 10−2	7.01 × 10−2	1.44 × 10−3	2.34 × 10−3
	III-2	1.56 × 10−2	1.78 × 10−2	5.15 × 10−2	7.01 × 10−2	1.09 × 10−3	1.57 × 10−3
	III-3	9.13 × 10−2	1.15 × 10−1	2.70 × 10−1	3.67 × 10−1	4.44 × 10−3	3.48 × 10−3
Case IV	IV-1	3.05 × 10−2	3.99 × 10−2	3.77 × 10−2	5.59 × 10−2	1.03 × 10−3	1.43 × 10−3
	IV-2	3.06 × 10−2	3.96 × 10−2	1.09 × 10−1	1.69 × 10−1	1.03 × 10−3	1.46 × 10−3
	IV-3	3.03 × 10−2	3.97 × 10−2	3.72 × 10−1	5.40 × 10−1	1.05 × 10−3	1.45 × 10−3
Case V		/	/	4.70 × 10−2	4.36 × 10−2	3.87 × 10−3	4.12 × 10−3

.

The worst solution is found with the TF mesh (Figure 3b); it contains large numerical errors and dispersions across the peaks of the terrain due to the misalignment of the velocity field with the VLs. Furthermore, the tracer cannot recover its original shape at the end of the integration period, leading to the corresponding max(E_a) = 0.67 and l₂ = 0.169 (Figure 4b and the black dashed-curve in Figure 4d). There are some numerical artifacts at the periphery of the tracer in the case of the TF mesh (Figure 3b) as well as at the periphery of the tracer of the cut-cell mesh (Figure 3a). Compared with those of the TF mesh, the errors of the cut-cell mesh are reduced by one order of magnitude, with a max(E_a) = 0.076 and l₂ = 0.0396 (Figure 4a and the blue dotted-curve in Figure 4d). As is evident from Figure 3c, the use of the adaptive mesh results in the lowest numerical dispersion and the least spurious tracer deformation. The adaptive mesh simulations achieve the most accurate results, with a max(E_a) = 0.0026 and l₂ = 0.00146.

A comparison of the absolute E_a distributions at the end of the integration period (Figure 4a–c) reveals that the cut-cell and TF meshes exhibit error extrema at both the center and at the periphery of the tracer, while the adaptive mesh displays an error maximum only in the forward direction of the tracer. However, as the simulation evolves, the direction of the mesh adaptation becomes more difficult to predict in a stream-wise non-uniform velocity field than in a linear horizontal velocity field. Consequently, the capability of the adaptive mesh in a nonlinear velocity field must be further assessed, as reported in the next section.

4.1.2. Case II: Terrain-Following Advection Test

In this stream-wise non-uniform advection test (Figure 2b), a TF velocity field [5] is considered that is given by

$$u\left(x,z\right)=-\frac{\partial \mathsf{\Phi }}{\partial z}=u_0\frac{H_{top}}{H_{top}-h_5\left(x\right)}$$

(21)

$$w\left(x,z\right)=\frac{\partial \mathsf{\Phi }}{\partial x}=u_0{H}_{top}\frac{\partial h_5}{\partial x}\frac{H_{top}-z}{{\left[H_{top}-h_5\left(x\right)\right]}^2}$$

(22)

where u₀ = 10 m/s, $\mathsf{\Phi }=-u_0{H}_{top}\frac{z-h_5\left(x\right)}{H_{top}-h_5\left(x\right)}$ is the stream function of the velocity field, and H_top = 25 km is the height of the top of the computational domain.

The results for the tracer concentrations at t = 0, 5000 and 10,000 s are shown in Figure 5, and the corresponding errors are illustrated in Figure 6 and Table 1. The shape of the tracer is distorted by the TF velocity field as the tracer moves over the terrain, after which it almost recovers its original shape at the end of the integration period. However, at the end of the integration period, the misalignment of the TF velocity field with the horizontal VLs induces larger advection errors in the case of the cut-cell mesh. Consequently, severe dispersion is readily visible (Figure 5a). In the case of the TF mesh, the tracer is only slightly dispersed, as is demonstrated by the absolute E_a distributions and l₂ error norms of O(10⁻²) (Figure 5b and the black dash-curve in Figure 6d), which are one order of magnitude less than those in the case of the cut-cell mesh (blue dotted-curve in Figure 6d). The E_a pattern in the case of the adaptive mesh resembles that of the TF velocity field (Figure 5c). Although extreme errors mainly appear at the center and periphery of the tracer in the adaptive mesh modelling scheme due to the variations in the forward tracer direction in the TF velocity field, the absolute E_a distributions and l₂ error norms are significantly reduced to O(10⁻³) (Figure 6c and the red solid-curve in Figure 6d). Therefore, the adaptive mesh modelling approach produces a more accurate solution than the cut-cell and TF meshes.

4.1.3. Case III: Advection Tests with Varying Velocity Directions

From the horizontal and TF advection tests performed above, the misalignment of the velocity field with the VLs has been found to induce more advection errors with the cut-cell and TF meshes than with the adaptive mesh. However, in the real atmosphere, the advection direction is not always that of a stratified fluid, especially in the vicinity of terrain. Building on the two advection tests described above, the effects of different velocity directions are investigated in detail in this section. Here, two time-invariant velocity fields with different non-horizontal directions and a wavelike velocity field (Figure 2c–e, respectively) are employed to evaluate the performances of the cut-cell, TF and adaptive meshes.

The time-invariant velocity fields are given by

$$u\left(x,z\right)=u_0\cos \varphi$$

(23)

$$w\left(x,z\right)=u_0\sin \varphi$$

(24)

where u₀ = 10 m/s and φ represents the angle of intersection between the velocity direction and the horizontal direction. For φ ≠ 0, the velocity penetrates the terrain surface. This test does not represent a physical advection velocity field, but it provides a numerical test for assessing the scalar transport in a non-horizontal advection velocity field. The wavelike velocity field is given by Equations (17), (21) and (22) with λ = 16 km. The vertical profiles of the velocity fields and the routes of the corresponding analytical tracer concentrations q are displayed in Figure 2c–e.

In these cases, the directions of the velocity fields are no longer aligned with VLs when the cut-cell or TF mesh is used. Consequently, these tests determine whether the adaptive mesh modelling accuracy is sensitive to the velocity direction. The direction of the first velocity field is east by north at a bearing of 2 degrees to the horizontal direction, meaning that φ = 2 degrees in Equations (23) and (24). The direction of the second velocity field is east by south at a bearing of 2 degrees to the horizontal direction, meaning that φ = −2 degrees in Equations (23) and (24). To ensure that the tracer exists throughout the entire computational domain, z₀ = 4.5 km and z₀ = 7.5 km are used in Equation (14) for the first and second velocity fields, respectively. The third velocity field is a wavelike velocity field which has three crests and two valleys and a slope gradient that is larger than that of the TF velocity field.

A comparison of the absolute E_a distributions s of the analytical tracer concentrations among the three meshes reveals that the TF mesh produces the worst solutions because it suffers from the largest misalignments of the corresponding velocity fields with respect to the TF VLs. The E_a distributions and l₂ error norms of the adaptive mesh models are one or two order(s) of magnitude less than those of the cut-cell and TF meshes for these velocity fields (3rd to 5th rows in Table 1 and Figure 7). Moreover, compared with the results for the linear fields, the advection errors obtained with the cut-cell and TF meshes for the non-stratified wavelike velocity field are one order of magnitude larger. However, the adaptive mesh modelling scheme still exhibits relatively small absolute errors and l₂ errors of O(10⁻³), without any unexpected degradation in accuracy regardless of the velocity direction (3rd–5th rows in Table 1 and Figure 7). Figure 7 also reveals that large error values appear at both ends in the forward direction of the tracer concentration q in the case of the adaptive mesh. Generally, the adaptive mesh modelling scheme produces better predictions than either the cut-cell mesh or the TF mesh, and the l₂ error norms in the adaptive mesh case are less sensitive to the direction of the velocity field.

4.1.4. Case IV: Advection Tests with Different Types of Terrains

Zou et al. [27] verified that increasingly higher terrain crests, which result in steeper terrain profiles (Figure 8), can induce larger advection errors in TF meshes. However, the proposed adaptive meshes are affected only by the given tagging levels and the Hessian matrix of specific variables (e.g., the tracer concentration in this study). Here, a one-crest terrain, a five-crest Schär-type wavelike terrain and a ten-crest wavelike terrain are used to investigate the effects of complex terrains on adaptive meshes. The one-crest bell-shaped terrain is defined as Equation (17). Here, λ = 32 km is used for the one-crest bell-shaped terrain while λ = 8 km is used for the five-crest Shär-type wavelike terrain h₅(x). The Zou-type ten-crest wavelike terrain [27] is defined as follows:

$$h_{10}\left(x\right)={\sin }^2\left[\frac{n\pi }{2a}\left(x-150\right)\right]h^{**}\left(x\right)$$

(25)

where

$$h^{**}\left(x\right)=\{\begin{array}{ll}{h}_0{\sin }^2\left[\frac{\pi }{a}\left(x-150\right)\right],& \left|x-150\right|\le a\\ 0,& \left|x-150\right|>a\end{array}$$

(26)

where h₀ = 3 km is the maximum height of the terrain, a = 25 km is the half-width of the terrain, and n = 10 is the number of terrain peaks. The slope angles ($\alpha =\arctan \frac{dh}{dx}$) of the three terrains are shown in Figure 8. Because the direction of the velocity field is horizontal in this section, the velocity field coincides with the VLs of the cut-cell mesh while the intersection angle between the velocity field and VLs at the bottom in the TF mesh is equal to α. With the increasing complexity of the terrains, the slope of the corresponding VLs in the vertical direction increases leading to the severe misalignment of the velocity fields with the VLs. The velocity field u is the same as the constant velocity field in Equations (19) and (20) in Case I. The absolute E_a distributions and l₂ error norms are shown in Figure 9 and Table 1, respectively.

First, the maximum absolute E_a values and l₂ error norms at t = 10,000 s appear in the case of the TF mesh because the misalignments of the velocity fields with respect to the VLs are largest in this case (Figure 8). More accurate results are obtained with the cut-cell mesh, due to an improved alignment of the flow with the grid. The adaptive mesh achieves the highest simulation accuracy, with errors that are one or even two order(s) of magnitude less than those displayed by the cut-cell and TF meshes (6th–8th rows in Table 1). Second, the advection errors obtained with the TF mesh increase as the slope of the VLs increases (Figure 8 and middle column in Figure 9). In contrast, as the complexity of the underlying terrain increases, the errors in the area over which the tracer passes in the cases of the cut-cell and adaptive meshes are minimally affected by changes in the terrain (left and right columns in Figure 9), thereby revealing that the errors resulting from the use of these meshes are less sensitive to the terrain. Third, compared with the results presented in Section 4.1.1, the order of and magnitude of the adaptive mesh modelling is the same as in these tests. Therefore, using adaptive meshes produces more accurate results in advection simulations than using cut-cell and TF meshes in cases with variable terrain, which has a lower effect on simulations carried out with adaptive and cut-cell meshes when the tracer is above the peak of the terrain in this case.

4.1.5. Case V: An Advection Test along the Terrain Surface

In this section, the tracer concentration q at t = 0 s is moved down to lie tangential to the surface of the terrain to test the effects of the lower boundary and the underlying terrain on the adaptive mesh. The setup for this case is described by Equation (14), except that z₀ = 3 km is used to set the tracer concentration q at t = 0 s closer to the surface of the terrain. The stratified TF velocity field described by Equations (21) and (22) is utilized to ensure the advection of the tracer along the terrain throughout the entire simulation. In this case, the lower boundary of the tracer is entirely tangential to the surface of the terrain.

Shaw and Weller [5] has demonstrated that performing horizontal advection simulations with cut-cells in the vicinity of the terrain is not as good as that with the TF mesh. Therefore, the effects of adaptive mesh are only compared with TF meshes hereafter. The TF and adaptive mesh results are compared in the middle (t = 5000 s) of the integration period. The absolute E_a distributions obtained using both meshes are illustrated in Figure 10a,b, and the l₂ error norms throughout the entire integration period are shown in Figure 10c. In the TF mesh, large errors appear throughout the region where the tracer is located. It is evident that the use of an adaptive mesh substantially reduces these advection errors throughout most of the tracer region, and large errors appear only in the center and at the periphery of the tracer. In Figure 10c the l₂ error norms in the adaptive mesh are also one order of magnitude smaller than those in the TF mesh. Therefore, the adaptive and TF meshes are both capable of representing the lower boundary of the computational domain when the tracer is tangential to the terrain surface because the use of a TF mesh ensures that the direction of the tracer flux is tangential to the surface of the terrain, while the adaptive mesh dynamically adapts to resolve complex terrain and flow features as the flow evolves. In other words, both the TF and adaptive mesh schemes based on the DG-FEM can produce physical results near a solid wall.

4.2. Mesh Refinement and Central Processing Unit (CPU) Cost

The relationship between the mesh refinement and the corresponding computational costs among cut-cell, TF and adaptive meshes are presented in this section. By giving a set of refinement tagging criterion levels $\widehat{\epsilon }$ in Equation (1), the horizontal advection test of Case IV is redone. The refinement tagging criterion levels of the adaptive mesh is set to 0.1, 0.05, 0.01, 0.005 and 0.001 in Equation (1) to investigate its effect on the mesh resolutions (the minimum mesh size), the computational cost (the runtime) and the node numbers resolving the tracer which is linked to the number of terrain peaks. The resolutions of cut-cell and TF meshes are set to 0.125 km, 0.25 km, 0.5 km and 1.0 km for the comparison to the adaptive mesh. In this study, 1-crest, 5-crest and 10-crest terrains (i.e., one wave, 5-waves and 10-waves) are used to represent the terrain complexity. The results are shown in Figure 11.

With the given tagging criterion levels in the adaptive meshes, the accuracy of the results can be achieved for all types of terrains with similar node numbers, the runtimes (Figure 11a) and the minimum mesh sizes (Figure 11b). This result is achieved because the mesh is adapted according to the tracer features passing over the terrains such that a high-resolution dense mesh is used only in areas of high tracer concentrations and a coarse mesh is used throughout the rest of the domain. The temporal and spatial adaptation of meshes ensures that the number of nodes and the corresponding computational costs used for all three types remain almost the same. However, for the cut-cell and TF meshes, the number of nodes should be linked to the wavenumber (Cases IV and V, see details in [27,28]). Although the area of the tracer is small compared to the computational domain, high-resolution meshes are required to achieve the increasing accuracy over the entire domain due to the fixed meshes. This increases the total number of the nodes during the simulation period, which leads to a significant increase in the computational cost. Therefore, the CPU time required in both the cut-cell and TF methods is more than 10 times as that in the adaptive mesh modelling when the given tagging criterion level is smaller than O(10⁻¹).

The proportion (hereinafter referred to as the proportion for simplicity) between the node numbers used for resolving the tracer and each wavelength is used for investigating the effect of the terrain complexity (Figure 11c). It can be seen that to achieve the same tagging criterion level, the proportion is increased with the increasing of the wavenumbers (representing the terrain complexity) in these three meshes. Compared to the cut-cell and TF meshes, the proportion required in adaptive mesh modelling is smaller in all the types of terrains. It reveals that for a given node numbers per wavelength, adaptive mesh modelling requires a less number of nodes for resolving the tracers than other two meshes, thus saving the CPU time. It is also noted that with the increased wavenumber, the corresponding growth rate of the proportion in the adaptive mesh is smaller than that in the cut-cell and TF meshes.

Therefore, to achieve the same tagging levels, the adaptive mesh requires fewer nodes, less runtime and lower proportion for the advection simulations through the entire computational domain, thus reducing the computational cost.

5. Discussion and Conclusions

Advection errors pose a significant computational problem in atmospheric modelling over steep terrains. Such errors are induced primarily by the misalignment of the flow direction with respect to the VLs, distortions of the grid and the complexity of the underlying terrain. To reduce these errors and to avoid wasting computational resources, this study introduces adaptive meshes that adapt to variations in the tracer concentration as the flow evolves and become sparser at locations where the gradients of variable solutions are relatively low. Through a series of 2D idealized advection experiments, the ability of these adaptive meshes to reduce advection errors was investigated.

To evaluate the performance of adaptive mesh modelling, five velocity fields and three different terrain types were employed. In the 2D simulations, adaptive meshes achieved higher accuracy in all advection tests. This is substantiated by the absolute values of the absolute errors and l₂ error norms being approximately one to two order(s) of magnitude smaller than those obtained by using cut-cell or TF meshes. Moreover, the adaptive meshes resulted in the lowest numerical dispersion above the peaks in the terrain profiles. Additionally, compared with the errors in the cut-cell and TF meshes, the errors in the adaptive meshes were less sensitive to both the direction of the velocity field and the complexity of the terrain. Finally, the adaptive meshes reduced the advection errors and were able to represent the underlying terrain without severe numerical dispersions where the tracer passed in close proximity to the surface of the terrain tangentially. All these results demonstrate that the time-dependent spatial density of an adaptive mesh ensures the accuracy and effectiveness of the advection simulation and that the accuracy of such a mesh is comparable with that of a cut-cell mesh of higher uniform spatial density (see details in [5]). Finally, the computational cost is analyzed. To achieve the same tagging criterion levels, the adaptive mesh requires fewer nodes, smaller minimum mesh sizes, less runtime and lower proportion between the node numbers used for resolving the tracer and each wavelength for the advection simulations than other mesh techniques, thus reducing the computational cost.