2.1 Spherical centroidal Voronoi tessellation

We start by defining the spherical Voronoi concepts that will be important for the local refinement technique discussed in this work.

Given a set of distinct points on the sphere that we call generators, $\mathit{\{}{\mathit{x}}_i{\mathit{\}}}_{i=\mathrm{1}}^n\subset {\mathcal{S}}^{\mathrm{2}}$, we define the ith Voronoi region as

where d(x,y) is the geodesic distance on the sphere, $\forall \mathit{x},\mathit{y}\in {\mathcal{S}}^{\mathrm{2}}$. The sets $\mathit{\{}\mathrm{\Omega }{\mathit{\}}}_{i=\mathrm{1}}^n$ are called spherical Voronoi tessellation (SVT) and decompose the sphere into disjoint convex spherical polygons to generate a spherical grid (Ju et al., 2011).

With the Voronoi tessellation, we can associate a dual grid by connecting nearest-neighbor generators. The geodesic segment connecting neighbors is orthogonal to the sharing edge between these two generators. This orthogonality is important for the numerical schemes that we shall present later. Under the assumption that no four near neighbors are co-circular with empty an interior, these dual cells all consist of triangles and form a Delaunay triangulation (Okabe et al., 2000). In this duality relation, each Voronoi region corresponds to a unique dual-cell vertex, each Voronoi vertex is associated with a unique dual cell and each Voronoi edge is associated with a unique dual edge.

An example of SVT is the primal grid of the icosahedral grid, the pentagonal–hexagonal grid. The icosahedral grid is obtained by gnomonically projecting the edges of an icosahedral inscribed within a sphere. This process generates 12 vertices and 20 triangles. Each of these triangles can be divided into four new triangles and we can proceed in this manner recursively. Given a level of recursion $l=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\cdots }$, it can be shown that the icosahedral grid obtained by recursion has $N_l=\mathrm{10}\cdot {\mathrm{2}}^{\mathrm{2}l}+\mathrm{2}$ vertices, and consequently, the dual pentagonal–hexagonal grid has N_l cells. The number l is also called the grid level.

A particular case of SVT that we are interested in is the spherical centroidal Voronoi tessellation (SCVT). Given a density function $\mathit{\rho }:{\mathcal{S}}^{\mathrm{2}}\to ]\mathrm{0},\mathrm{\infty }[$, for each Voronoi region Ω, the center of mass with respect to ρ is given by

For this definition of the center of mass, it does not necessarily hold that ${\mathit{x}}^{\ast }\in {\mathcal{S}}^{\mathrm{2}}$. Therefore, following Du et al. (2003), we introduce the concept of a spherical center of mass, which may be calculated using radial projection as

An SCVT is an SVT such that each generator is the center of mass of the corresponding Voronoi region with respect to a given density function ρ. As we can notice, SCVT generates spherical grids that we call SCVT grids. An algorithm that creates SCVT from initial random generators is the iterative Lloyd's method (Du et al., 1999, 2003), which iteratively pushes the Voronoi generators to the Voronoi mass centroid position, and is the method employed in this work.

Given two Voronoi cells Ω_i and Ω_j, we denote their diameter by h_i and h_j, respectively. It is possible to show that (Du et al., 2003; Ju et al., 2011)

This estimation plays a key role in this work because it allows us to build grids with local refinement.

We show examples of SCVTs for ρ=1 generated using Lloyd's method in Fig. 1 where we use the pentagonal–hexagonal grid as an initial guess in Lloyd's method. The corresponding dual grid is showed in Fig. 2. We use N_l generators for $l=\mathrm{2},\mathrm{3},\mathrm{4}$. The dual grid of an SCVT grid may be seen as an optimized version of the icosahedral grid (Miura and Kimoto, 2005).

Figure 1SCVT grids using a uniform density function for 12 (a), 42 (b) and 162 (c) generators.

Download

Figure 2Dual (triangular) grids of the grids shown in Fig. 1 for 12 (a), 42 (b) and 162 (c) generators.

Download

Next, we present how to build a grid with local refinement on the Andes Mountains based on Eq. (4), with the density function presented in Ju et al. (2011) and the topography data from Amante and Eakins (2009).

2.2 Locally refined SCVT

We represent the ℝ³ coordinates of a point x∈𝒮² in spherical coordinates as

and define the parameters $({\mathit{\varphi }}_{\mathrm{c}},{\mathit{\lambda }}_{\mathrm{c}})=\left(-\frac{\mathit{\pi }}{\mathrm{9}},-\frac{\mathit{\pi }}{\mathrm{3}}\right)$ and ${\mathit{x}}_{\mathrm{c}}=(\sin {\mathit{\varphi }}_{\mathrm{c}}\cos {\mathit{\lambda }}_{\mathrm{c}},\sin {\mathit{\varphi }}_{\mathrm{c}}\sin {\mathit{\lambda }}_{\mathrm{c}},\cos {\mathit{\varphi }}_{\mathrm{c}})$. We also consider real values γ, α and ε that will be set later. Denoting the Euclidean norm as $\Vert \cdot \Vert$, we define the distance from x_c by $d(\mathit{x},{\mathit{x}}_{\mathrm{c}})=\Vert \mathit{x}-{\mathit{x}}_{\mathrm{c}}\Vert$ for x∈𝒮² and introduce the following auxiliary function:

A density function that will allow us to refine based on two criteria, in our case the topography (Andes Mountains) and the continent (South America), is given by

where b represents the topography (in our case, the Andes region), in this work retrieved from Amante and Eakins (2009) and normalized for the interval [0,1], γ≥1 and $\mathit{\mu }\in [\mathrm{0},\mathrm{1}]$. Note that for μ=1, the density function will refine only South America without capturing the Andes' sharp topography. For μ=0, the density function will refine only the Andes' sharp topography. Therefore, for a value $\mathit{\mu }\in ]\mathrm{0},\mathrm{1}[$, the generated grid will capture both the Andes' sharp topography and the South American continent and will have a transition to the coarse global grid. The closer the parameter μ is to 1, the less the topography will be represented in the refinement. And the closer the parameter μ is to 0, the more the sharp topography will be represented and the lower the refinement will be on the rest of the South American continent. With α, μ and s(x) properly specified, the parameter γ will approximately represent the ratio between coarse- and fine-grid-resolution areas, which follows from Eq. (4). We will provide further details at the end of this section.

For the numerical schemes that will be adopted in this work, sometimes termed mimetic finite-volume schemes, grid quality plays a major role. On staggered (C-type, using the meteorology terminology) grids, different variables are placed on different positions of the grid, and operators, differential or not, act on taking variables from one place to another. On highly distorted grids, this can lead to loss of accuracy of some operators (see, for instance, Peixoto and Barros, 2013) or even inadequate physical representation of a model field. One of the key properties of interest is the position of the circumcenters of the Delaunay triangles, which are also the Voronoi cell vertices. Delaunay grids have the so-called well-centered property if all circumcenters lie in the interior of their respective triangles (see, for instance, Engwirda, 2017). Here, we will consider grid smoothing techniques to ensure grids with the well-centered property, as we will discuss below.

The ETOPO data (Amante and Eakins, 2009) used in the function b are defined in a latitude–longitude grid with 720×1440 points, which is enough resolution to represent the Andes' topography well for our experiments. We employed bilinear interpolation to compute the function b in SCVT grid points that are not in the latitude–longitude grid. In order to guarantee grids with smooth transitions between the high- and low-resolution regions, we applied a smoothing technique on ETOPO data based on the Jacobi method. For each index (i,j) representing a point on the latitude–longitude grid, $\mathrm{0}\le i\le \mathrm{720}$, $\mathrm{0}\le j\le \mathrm{1440}$, we apply the classic Jacobi iteration process:

where $k=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\cdots }$ and $b_{ij}^{\mathrm{0}}$ is the initial Andes topography data from ETOPO. The number of iterations was empirically defined as 500 steps, ensuring that the basic shape of the Andes Range is still well represented in the topography dataset. More importantly, this smoothing later ensures that the generated SCVT primal grids all have circumcenters inside the triangles.

We also applied a smoothing filter in the function s, considering it defined in the latitude–longitude grid with 720×1440 points. We used a moving average method by replacing the value of s at a point in the latitude–longitude grid considering the mean of s in a box with 55×55 points centered at the considered point. This box size has been shown to guarantee that the generated grids have circumcenters inside the triangles.

Figure 3SCVT grid with local refinement on the Andes using 10 242 generators (a). Cell diameters in kilometers for the grid with (b) 10 242 and (c) 163 842 generators considering γ=3, $\mathit{\alpha }=\frac{\mathrm{7}\mathit{\pi }}{\mathrm{45}}$, $\mathit{\epsilon }=\frac{\mathit{\pi }}{\mathrm{12}}$ and μ=0.8.

We set the parameters to γ=3, $\mathit{\alpha }=\frac{\mathrm{7}\mathit{\pi }}{\mathrm{45}}$, $\mathit{\epsilon }=\frac{\mathit{\pi }}{\mathrm{12}}$ and μ=0.8. The parameter ε, which represents the width of the circular refinement region, and the parameter γ were defined as in Ju et al. (2011). The parameter α was defined in such a way to guarantee that the circular refinement region contains the South American continent, and μ was empirically chosen aiming to represent the Andes' shape well. In Fig. 3, we show the grid generated with 10 242 nodes after applying Lloyd's method. The diameters of cells are also shown in Fig. 3, including the diameters of the grid generated with 163 842 nodes. The diameter of a Voronoi cell is estimated as the average distance from its generator to its neighbors. As we can notice, the ratio of the diameters of the cells on the Andes and the cells in the coarser-grid region is approximately 3. We can derive this from Eqs. (4) and (7). Indeed, observe that if x_i∈𝒮² is a Voronoi generator such that $d({\mathit{x}}_{\mathrm{c}},{\mathit{x}}_i)>\mathit{\alpha }+\mathit{\epsilon }$, i.e., x_i is in the uniform-resolution region, we have $s\left({\mathit{x}}_i\right)=b\left({\mathit{x}}_i\right)=\mathrm{0}$ and therefore $\mathit{\rho }\left({\mathit{x}}_i\right)=\frac{\mathrm{1}}{{\mathit{\gamma }}^{\mathrm{4}}}$. On the other hand, notice that if x_j∈𝒮² is a Voronoi generator and a point in the Andes region, we have s(x_j)=1 and b(x_j)≈1, and therefore ρ(x_j)≈1 using Eq. (7). At last, using Eq. (4), we can see that $\frac{h_i}{h_j}\approx \frac{\mathrm{1}}{\mathit{\gamma }}$, where h_i and h_j are the diameters of the Voronoi regions i and j, respectively.

Table 1Minimum and maximum diameters of the Voronoi cells for the uniform SCVT grid (UR) and for the variable-resolution (VR) grid that captures the Andes well. The VR grids considers the parameters γ=3, $\mathit{\alpha }=\frac{\mathrm{7}\mathit{\pi }}{\mathrm{45}}$, $\mathit{\epsilon }=\frac{\mathit{\pi }}{\mathrm{12}}$ and μ=0.8.

Download Print Version | Download XLSX

We generate the variable-resolution SCVT grids using the density function defined in Eq. (7) using $N_l=\mathrm{10}\cdot {\mathrm{2}}^{\mathrm{2}l}+\mathrm{2}$ generators for $l=\mathrm{1},\mathrm{2},\mathrm{\cdots },\mathrm{7}$, and we refer to these grids as VR1, VR2, ⋯, VR7, respectively. These generated grids and the uniform SCVT grids will be used in our numerical experiments in Sect. 4. The uniform-resolution SCVT grids were generated from grid level 1 up to 9, and we refer to these grids as UR1 to UR9. We show the minimum and maximum diameters of the cells for the uniform- and variable-resolution SCVT grids in Table 1. We also generated locally refined grids with a 1 : 6 ratio of diameters between the coarse grid and cells in the Andes region with γ=5, $\mathit{\alpha }=\frac{\mathrm{7}\mathit{\pi }}{\mathrm{45}}$, $\mathit{\epsilon }=\frac{\mathit{\pi }}{\mathrm{12}}$ and μ=0.8. These grids will be analyzed in our numerical experiments only in the Sect. 4.2.1 aiming to investigate the impact of increasing the resolution of the refined region.

3.1 Shallow-water model

As is usual in the development of numerical methods for global atmospheric models, we first evaluate it considering the shallow-water equations on the sphere. These equations represent important aspects of the atmosphere, such as Rossby waves, the Coriolis effect and geostrophic adjustment, and they have the advantage of being two-dimensional; therefore, they are computationally cheaper to run at high resolutions.

Table 2Variables and parameters for the shallow-water model.

Download Print Version | Download XLSX

We denote a point on the sphere by x and a time instant by t≥0. Then, we consider the nonlinear shallow-water equations on the sphere written in the vector-invariant form (Ringler et al., 2010):

where $h=h(\mathit{x},t)$ denotes the total fluid depth, $\mathit{u}=\mathit{u}(\mathit{x},t)$ is the velocity vector tangent to the sphere for all t≥0 and all other terms are defined in Table 2. The operators ∇ and ∇⋅ are assumed to be the horizontal gradient and divergence on the sphere, respectively. These equations define an initial value problem that describes the evolution of a thin layer of a two-dimensional fluid on the sphere.

Figure 4Variable positions for a C-grid staggering.

Download

The finite-volume scheme adopted here was proposed in Thuburn et al. (2009) and Ringler et al. (2010), also denoted as TRSK, and considers a grid with an orthogonal dual grid and a C-grid staggering of the variables. We consider three positions in the C-grid: primal-cell centers (or dual-cell vertices), primal-cell vertices (or dual-cell centers), and the intersection point between a primal edge and dual edge that we call edge points (Fig. 4).

Table 3Grid notation following Ringler et al. (2010).

Download Print Version | Download XLSX

We summarize the notation needed for TRSK in Table 3. In TRSK, the height h is stored at primal centers and is denoted by $h_i\left(t\right)=h({\mathit{x}}_i,t)$. Given an edge e, we define a normal vector n_e to the edge pointing to the direction such that the normal component of the velocity given by $u_e\left(t\right)=\mathit{u}({\mathit{x}}_e,t)\cdot {\mathit{n}}_e$ is positive. The normal component of the velocity is stored at edge points. Our prognostic variables in TRSK are h_i and u_e for each primal center i and edge e. Taking the dot product of Eq. (9) with n_e, we get

where $[\mathit{F}{]}_e=\mathit{F}({\mathit{x}}_e,t)\cdot {\mathit{n}}_e$ and $[G{]}_i=G({\mathit{x}}_i,t)$ for any vector field F and scalar field G. Each term on the right-hand side of Eqs. (11) and (12) is approximated using the spatial discrete operators from Ringler et al. (2010). After discretizing each term for all primal cells and edges, we get an ordinary differential equation system. In this work, we employed a fourth-order Runge–Kutta scheme to solve the ordinary differential equation system since we are aiming to analyze only the impact of the local refinement on the horizontal discretization. At last, TRSK has been known for its mimetic properties; i.e., the discrete equations mimic some of the continuous equations, such as

1. mass conservation,

2. the Coriolis term not being an energy source or sink $(u_e{u}_e^{⟂}=\mathrm{0})$,

3. total energy being conserved within time error truncation, and

4. preservation of stationary geostrophic modes on the f sphere.

The TRSK discretization is also known to be very low-order accurate and may even lack consistency on nonregular grids (Peixoto and Barros, 2013; Peixoto, 2016). However, due to its mimetic properties and flexibility, it is still adopted in some modern global atmospheric models such as MPAS (Skamarock et al., 2012) and DYNAMICO (Dubos et al., 2015).

3.2 Moist shallow-water model

Even though the shallow-water model retains keys properties of the atmosphere, it lacks the representation of sub-scale physical processes. Therefore, we consider here an intermediate step of a two-dimensional moist shallow-water model, avoiding the computational expense of a three-dimensional primitive equation model while being able to account for water-related sub-grid processes.

Table 4Variables and parameters for the moist shallow-water model (Zerroukat and Allen, 2015).

Download Print Version | Download XLSX

We follow the moist shallow-water model proposed by Zerroukat and Allen (2015), wherein moist dynamics are included in the classic shallow-water model. This model is two-dimensional, which is computationally cheaper than three-dimensional models at high resolutions, and is derived from the Boussinesq approximation of the three-dimensional Euler equations. The model can be summarized in the following equations:

where θ is the temperature, $k=\mathrm{1},\mathrm{2},\mathrm{3}$, q¹=q_v is the water vapor state, q²=q_c is the cloud state and q³=q_r is the rain state. We are using the flux form for the advection equations for numerical purposes that we shall see later. The moisture variables are tracers that are advected in this system. The other variables are the same as in the shallow-water model (Table 2). The source terms of the advection equations define a three-state moist physics, wherein vapor may be converted to clouds, clouds are allowed to evaporate and clouds may be converted to rain. These three-state moist physics depend on a saturation function, which is given by Zerroukat and Allen (2015) as

and the initial vapor state is given by $q_{\mathrm{v}}(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=G(\mathit{\lambda },\mathit{\varphi })q_{\mathrm{sat}}\left(\mathit{\theta }\right)$, where λ denotes the longitude and q₀ is a constant chosen in order to guarantee that the initial maximum value of q_v is 0.02 and G is a function between 0 and 1. The additional variables and parameters to the moist shallow-water model are shown in Table 4. Note that $S_{q_{\mathrm{v}}}+S_{q_{\mathrm{c}}}+S_{q_{\mathrm{r}}}=\mathrm{0}$; therefore, the integral of $h(q_{\mathrm{v}}+q_{\mathrm{c}}+q_{\mathrm{r}})$ is conserved in this model.

In this model, we consider the prognostic variables u, h, hθ and hq^k, $k=\mathrm{1},\mathrm{2},\mathrm{3}$. Similarly to the shallow-water model, we store h, hθ and hq^k at Voronoi centers, and the normal component of u is stored at the edges. We write equations similar to Eqs. (11) and (12) for the moist shallow-water equations and discretize the right-hand side using TRSK operators.

The divergence term in Eqs. (15) and (16) for hθ and hq^k may be computed using the divergence discretization of TRSK, and this explains why we used the advection equation in flux form. The source terms for Eqs. (15)–(16) are straightforward to calculate once they depend only on θ and q^k, and thus we only need to divide these variables by h. Equation (13) has a source term and it needs to be evaluated at edge points. This source may be calculated using TRSK interpolation and gradient operators. We point out that MPAS solves the tracer equations such as Eqs. (15)–(16) using a high-order advection scheme (Skamarock and Gassmann, 2011). Here, we consider the TRSK scheme to solve tracer equations aiming to analyze how TRSK simulates the formation of cloud and rain. Furthermore, as we shall see in Sect. 4.3, a monotonic (positivity-preserving) filter is employed in our simulations.

4.1 Hyperdiffusion

As we will show in the next sections, the variable-resolution grid triggers unstable modes of the numerical method and spurious waves are generated, destabilizing the numerical solution. While undesirable for the shallow-water equations, artificial numerical diffusion can be used to mitigate the problem. Here we discuss the numerical diffusion mechanisms that will be explored in this study. To ensure selectiveness of the diffusion mechanism, we will employ fourth-order diffusion (hyperdiffusion), but results with second-order diffusion may be found in an earlier version of the work, available from the journal open discussion material.

The numerical hyperdiffusion was included in the momentum equation (Eq. 11), following Skamarock et al. (2012) and Jablonowski and Williamson (2011), based on the following vector identity for the vector Laplacian:

Evaluating the above identity on an edge e and taking the dot product with corresponding normal vector n_e, we have

The first term on the right-hand side of Eq. (19) is calculated using the gradient at edges in the normal direction using the divergence given at Voronoi centers. The second term of the right-hand side of Eq. (19) is calculated as ∇ζ(x_e)⋅t_e (Skamarock et al., 2012), where ${\mathit{t}}_e=\mathit{k}\times {\mathit{n}}_e$, and is obtained as the gradient at edges in the tangent direction using the relative vorticity ζ given at primal vertex positions. The fourth-order hyperdiffusion ∇⁴ is obtained using Eq. (19) and the relations ∇²=Δ and ${\mathrm{\nabla }}^{\mathrm{4}}\mathit{u}={\mathrm{\nabla }}^{\mathrm{2}}\left({\mathrm{\nabla }}^{\mathrm{2}}\mathit{u}\right)=\mathrm{\Delta }\left(\mathrm{\Delta }\mathit{u}\right)$.

In order to introduce a variable coefficient of hyperdiffusion, we follow Klemp (2017) using the formula

where K₂ denotes the variable diffusion coefficient. This formula (20) is not exact. However, as Klemp (2017) shows, this formulation ensures a dissipative behavior. This formulation allows us to discretize the Eq. (20) just as we discretized Eq. (19). Given a variable hyperdiffusion coefficient K₄, we may obtain a fourth-order filter using Eq. (20) recursively with $K_{\mathrm{2}}=\sqrt{K_{\mathrm{4}}}$, as explained by Klemp (2017). After the hyperdiffusion in computed, we add it on the right-hand side of the momentum equation (Eq. 11) considering a negative sign to ensure the dissipative behavior the hyperdiffusion operator.

We have examined different possibilities of hyperdiffusion mechanisms. First, we considered constant and cell size (diameter-based) hyperdiffusion coefficients in Sect. 4.2.1. The diameter-based hyperdiffusion coefficient (K₄) used in this work follows Zarzycki et al. (2014) and, in a Voronoi cell i, is given by $K_{max}\cdot {\left(\frac{h_i}{h_{max}}\right)}^p$, where h_i is the diameter of the ith Voronoi cell, h_max is the maximum diameter of all cells and p=log ₂10. With this choice of p, the hyperdiffusion coefficient reduces 10 times when the cell diameter is divided by 2.

As will be shown in the numerical results, the dominant factor for larger errors in the numerical solution is the cell geometry. With the variable-resolution grid, the occurrence of more distorted cells is almost inevitable, and these are also responsible for triggering the numerical instabilities. Peixoto and Barros (2013) showed that one can detect such problematic cells following an “alignment index”, which measures how well the cell edges are aligned. Perfectly aligned cells have edges 2×2 parallel and of the same length, ensuring second-order accuracy of typical finite-volume operators. Peixoto (2016) shows how the grid alignment can affect the solution of the shallow-water equations, and Peixoto and Barros (2014) used the alignment index to build vector reconstructions that use different schemes depending on the cell geometry. Following similar reasoning, here we propose a variable hyperdiffusion that adapts the coefficient based on this alignment index. We therefore employ less diffusion in better-aligned cells, so we can be more selective in imposing the diffusion mechanism where it is really needed.

For the alignment-based hyperdiffusion, we considered a variable coefficient K₄, which, in a Voronoi cell, is given by a smoothed alignment index (normalized between 0 and 1) multiplied by a value K_max for each cell. The alignment index is defined in Peixoto and Barros (2013), and the value of 0 indicates a perfectly aligned cell; therefore, on almost perfectly aligned cells the diffusion coefficient will be close to zero. The value K_max indicates that the maximum of the fourth-order hyperdiffusion coefficient allowed is employed in the most ill-aligned cells. To avoid abrupt transitions in the diffusion mechanism, here the alignment index is smoothed by replacing its cell value by the average of its neighbors (including its own value). The smoothed alignment on edges and vertex positions may be obtained using TRSK interpolation operators.

The values of K_max used in this work were chosen after an investigation of the geostrophically balanced flow test case from Williamson et al. (1992), looking for the minimum coefficient that could minimize the error of the normal component of the velocity. We describe this experiment in detail in Sect. 4.2.1.

We highlight the fact that other (different) diffusion mechanisms were investigated in this work but are not shown here. Second-order Laplace diffusion was applied in the early stages of this work, with adequate stabilization of the model but with greater impacts on accuracy and energy conservation. Vorticity base stabilization mechanisms (Weller, 2012) were also explored, but these generally failed to stabilize the model (see specific comments in Sect. 4.2.3).

4.2 Shallow-water model

4.2.1 Steady geostrophic flow

This test case is proposed in Williamson et al. (1992) and is known as test case 2 (hereafter, TC2). The initial conditions are given by

$$\begin{array}{}\text{(21)}& {\displaystyle }h(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=h_{\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{g}}\left(a\mathrm{\Omega }{u}_{\mathrm{0}}+{\displaystyle \frac{u_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}}}\right){\sin }^{\mathrm{2}}\mathit{\varphi },\text{(22)}& {\displaystyle }u(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=u_{\mathrm{0}}\cos \mathit{\varphi },\text{(23)}& {\displaystyle }v(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=\mathrm{0}.\end{array}$$

The bottom topography is equal to zero. The functions defined above satisfy the shallow-water equations and represent a balanced geostrophic flow. Thus, the initial condition should remain constant. We define the parameters $h_{\mathrm{0}}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{3}}$ and $u_{\mathrm{0}}=\frac{\mathrm{2}\mathit{\pi }a}{\mathrm{12}\mathrm{d}}$, where a is the Earth's radius.

Figure 5Steady geostrophic flow test case: time evolution of relative error in the L_∞ norm for the height field (a) and the normal component of the velocity field (b). The blue lines indicate the errors for the UR7 grid, while the orange lines show the error for the VR7 grid with refinement based on topography, and the green lines show the error for the same VR7 grid but adding numerical fourth-order hyperdiffusion based on the alignment index with $K_{max}={\mathrm{10}}^{\mathrm{12}}$ m⁴ s⁻¹ in the shallow-water model.

Download

In Fig. 5 we show the evolution of the relative error in norm L_∞ for the grids of level 7 (UR7 and VR7). This result shows that the error for the normal component of the velocity field grows from early stages of the integration for the VR7 grid without numerical hyperdiffusion. This error triggers errors in the height field that start to grow significantly after day 15. This behavior does not occur in the UR7 grid. A similar error evolution as shown in Fig. 5 may be obtained for the VR6 grid (not shown). There is therefore a clear need for artificial numerical diffusion mechanisms for stabilization purposes.

Figure 6(a) Height field error using a constant hyperdiffusion coefficient and a (b) diameter-based hyperdiffusion coefficient for the steady geostrophic flow test case considering a VR6 grid at day 30 and $K_{max}={\mathrm{10}}^{\mathrm{13}}$ m⁴ s⁻¹ (see Sect. 4.1). (c) Alignment index (Peixoto and Barros, 2013) for the VR6 grid.

To decide on the best hyperdiffusion strategy, we now examine the different formulations: a constant coefficient, a diameter-based variable coefficient and an alignment-based variable coefficient. For the variable-resolution grid of level 6 (VR6), we show in Fig. 6 the height field error after 30 d considering constant and diameter-based hyperdiffusion coefficients considering $K_{max}={\mathrm{10}}^{\mathrm{13}}$ m⁴ s⁻¹. The error pattern of both solutions (Fig. 6a and b) has a clear correspondence with the alignment index (Fig. 6c). However, we note that the error patterns are not necessarily related to the grid resolution, indicating that, while diameter-based hyperdiffusion may be well justified from a physical point of view, it may not be ideal for stabilization requirements originating from grid irregularities.

Figure 7Steady geostrophic flow test case showing relative errors after 30 d for the VR6 grid using the different hyperdiffusion methods and different values of K_max. (a) L_∞ error for the height field. (b) L_∞ error for the normal component of the velocity field. (c) L₂ error for the height field. (d) L₂ error for the normal component of the velocity field. The blue lines consider a constant hyperdiffusion coefficient, the orange lines consider an alignment-based coefficient and green lines consider a diameter-based coefficient.

Download

In Fig. 7 we present the relative errors after 30 d for the height and velocity fields in L_∞ and L₂ norms considering the VR6 grid with respect to different hyperdiffusion coefficient choices (K_max) for the constant, diameter-based and alignment-based hyperdiffusion described in Sect. 4.1. We can notice from this figure that the error of the velocity field requires more hyperdiffusion to be stabilized for all the employed methods. Also, the error of the diameter-based coefficient requires higher values of K_max to stabilize the error of the height field. The stabilization behavior of the alignment-based method is more similar to the behavior of the constant coefficient, while imposing diffusion only where necessary. This motivated us to move forward with the alignment-based hyperviscosity for further experiments.

Figure 8Hyperdiffusion coefficient (in log 10 scale) based on the smoothed alignment index for each cell in the VR7 grid (a) and L_∞ error of the normal component of the velocity field for different values of K_max considering the steady geostrophic flow test case for the VR7 grid (b).

For a variable-resolution grid of level 7 (VR7), an optimal choice for the hyperdiffusion coefficient may be found at around $K_{max}={\mathrm{10}}^{\mathrm{12}}$ m⁴ s⁻¹, as shown in Fig. 8b, where we also show the distribution of the coefficient on the sphere (Fig. 8a). Note that large areas of the globe are dominated by coefficients orders of magnitude smaller than K_max (ranging around 10⁹ to 10¹⁰ m⁴ s⁻¹). With this approach, we see in Fig. 5 how the model is now stable in the period of 30 d, with the error evolution in the VR7 grid with numerical hyperdiffusion having behavior similar to the UR7 grid.

Figure 9Steady geostrophic flow test case – plots for the VR7 grid: (a) height field error for each grid cell after 30 d (without diffusion). (b) Height field error for each grid cell after 30 d using numerical fourth-order hyperdiffusion based on the alignment index in the shallow-water model. (c) Grid alignment index (Peixoto and Barros, 2013). Plots for the UR7 grid: (d) height field error for each grid cell after 30 d. (e) Grid alignment index.

Figure 9 shows the error after 30 d for the height field in the VR7 grid without (Fig. 9a) and with alignment-based hyperdiffusion (Fig. 9b). We also show the alignment index (Fig. 9c). From Fig. 9a, we can see that the error for the height field propagates similarly as gravity waves, although spurious, when no numerical fourth-order hyperdiffusion is employed. From Fig. 9b and c, it is clear that the error obtained for the height field is greater in the cells that are ill-aligned, and the same holds for the UR7 grid, as Fig. 9d–e show. This result is in agreement with the results from Peixoto (2016), wherein it is shown that the divergence discretization accuracy used in TRSK is related to the alignment index.

Following these analyses, we will consider in the next sections the hyperdiffusion always based on alignment and adopt for VR7 a coefficient $K_{max}={\mathrm{10}}^{\mathrm{12}}$ m⁴ s⁻¹. For VR6 we will adopt $K_{max}={\mathrm{10}}^{\mathrm{13}}$ m⁴ s⁻¹, and for all coarser grids (VR1 to VR5), $K_{max}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{13}}$ m⁴ s⁻¹ was enough to stabilize the model.

In this test case, the total energy is conserved with precision 10⁻¹⁰ for the UR7 grid and 10⁻⁹ for the VR7 grid without fourth-order hyperdiffusion and with precision 10⁻⁸ considering fourth-order hyperdiffusion on the VR7 grid. Overall, we can see that the addition of fourth-order hyperdiffusion did not deteriorate the total energy conservation of TRSK, since this dissipation method is a scale-selective damping mechanism (Jablonowski and Williamson, 2011). We point out that the second-order diffusion ∇² is also able to stabilize the errors in this test case. However, the operator ∇² is not as selective as ∇⁴ and leads to deterioration of total energy conservation. Indeed, we re-run this test for the VR7 grid using ∇² with constant diffusion coefficient equal to 8.2×10³ m² s⁻¹, as in spectral methods with truncation level T85 (Jablonowski and Williamson, 2011), and we observe that the errors remained stable and the total energy was conserved with precision 10⁻⁴ (not shown).

Figure 10Steady geostrophic flow test case: convergence of relative error in the L_∞ norm for the height field (a) and the normal component of the velocity field (b) with respect to errors at day 30. Panels (c) and (d) show the same as (a) and (b), respectively, but for the L₂ norm. The blue lines indicate the errors for the uniform SCVT grids, while the orange lines show the error for the variable-resolution grids with refinement based on topography, and the green lines show the error for the same variable-resolution grids but adding numerical fourth-order hyperdiffusion based on the alignment index in the shallow-water model.

Download

Furthermore, Fig. 10 shows that adding fourth-order hyperdiffusion in the variable-resolution SCVT grid simulations is needed to achieve better convergence in both L_∞ and L₂ norms. We can see that the error does not converge to zero for both height and velocity fields, even when we add the fourth-order hyperdiffusion. Nevertheless, lack of convergence is also observed in the uniform-resolution SCVT grids although in the maximum norm, as is shown in Peixoto (2016). We point out that Fig. 10 may suggest that the grids VR1 to VR5 do not develop spurious gravity waves as in the grids VR6 and VR7. Nevertheless, they form such waves but for longer times of integrations.

Figure 11(a) Diameters in kilometers for the variable-resolution grid with 1 : 6 coarse-to-fine cell ratio. The grid level is equal to 7. (b) Steady geostrophic flow test case showing the time evolution of relative error in the L_∞ norm for the height field considering the variable-resolution grid with a 1 : 6 coarse-to-fine cell ratio. The blue lines indicate the errors for the uniform SCVT grids, while the orange lines show the error for the variable-resolution grids with refinement based on topography with 1 : 6 coarse-to-fine cell ratio, and the green lines show the error for the same variable-resolution grids but adding numerical fourth-order hyperdiffusion based on the alignment index with $K_{max}={\mathrm{10}}^{\mathrm{12}}$ m⁴ s⁻¹ in the shallow-water model.

In order to analyze the impact of increasing the resolution of the refined region, we consider the 1 : 6 ratio of fine-to-coarse cell diameter grid as discussed in Sect. 2.2. In Fig. 11a we illustrate the diameters (in km) of the 1 : 6 variable-resolution grid with a grid level equal to 7. We ran the steady geostrophic flow for this grid. In Fig. 11b we show the time evolution of the relative error in the L_∞ norm. The behavior is similar to the one depicted in Fig. 7a; however, the error starts to grow at an earlier time. Also, the error after 30 d is greater for the 1 : 6 variable-resolution grid. Therefore, we can conclude that the error starts to develop earlier and is larger when we increase the ratio between fine and coarse cell diameters, which is in agreement with Liu and Yang (2017), wherein the authors investigate the impact of SCVT grids with a circular refinement region considering the classical shallow-water model. Their work suggests that ratios between high- and low-resolution cell diameters bigger than 1 : 3 produce unacceptable errors, recommending only 1 : 2 and 1 : 3 ratios. Finally, we observe that using the same amount of fourth-order hyperdiffusion based on the alignment index stabilizes the error, as the green line in Fig. 11b shows.

4.2.2 Flow over the Andes

This test case is an adaption of the flow over a mountain test case proposed by Williamson et al. (1992), wherein we replace the bottom topography b with the smooth Andes topography described in Sect. 2.2. The velocity field is initialized as in the previous test case of the steady geostrophic flow. The initial height field is then given by

$$\begin{array}{}\text{(24)}& h(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=h_{\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{g}}\left(a\mathrm{\Omega }{u}_{\mathrm{0}}+{\displaystyle \frac{u_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}}}\right){\sin }^{\mathrm{2}}\mathit{\varphi }-b(\mathit{\lambda },\mathit{\varphi }).\end{array}$$

We adopted h₀=5400 and normalized the bottom topography to lie in the interval [0,2000]. Since the exact solution for this test case is not available, we use as a reference solution the solution yielded by ENDGame, a semi-Lagrangian dynamical core model developed by Thuburn et al. (2010), in a latitude–longitude grid of 512×1024 points. We interpolate the values from the latitude–longitude grid to the grid points needed in SCVT grids using cubic interpolation.

Figure 12Flow over the Andes test case for the shallow-water model: height field cell-wise error plots for the VR7 grid after 1 h (a), 2 d (b), 3 d (c) and 30 d (d). We consider the ENDGame solution as a reference solution. We show the height field for the VR7 grid in (e) and the reference height field (f) after 30 d.

We ran this test case for 30 d on the uniform- and variable-resolution SCVT grids of level 7 (UR7 and VR7, respectively). In Fig. 12, we show the error of the height field considering the ENDGame solution as a reference. We can notice that the error is initially concentrated on the Andes Mountains (Fig. 12a), and it is propagated to the other regions of the grid (Fig. 12b–c) and finally on the whole grid (Fig. 12d). A similar result is obtained in Tian (2020), wherein the author considers this test case using a cosine bell topography instead of the smoothed Andes topography; it is shown that this test case generates errors that propagate as gravity waves, and these errors can be larger in the high-resolution grid region.

Figure 13Flow over the Andes test case for the shallow-water model: (a) time evolution of relative error in the L_∞ norm for the normal component of the velocity field for the UR7 grid (blue line) and for the VR7 with (green line) and without (orange line) numerical fourth-order hyperdiffusion based on the alignment index. (b) Cell-wise error for the normal component of the velocity field at day 30 for the VR7 grid without fourth-order hyperdiffusion. We consider the ENDGame solution as a reference solution.

Our simulation shows that, even though the Andes are represented more accurately in the VR7 grid, the height field errors are larger in the Andes in the first days of integration. Despite that, the magnitude of the error does not increase as in the steady geostrophic flow test case. However, Fig. 13a illustrates how the L_∞ error for the velocity field increases with time. Even though they increase, they do not trigger larger errors as in the steady geostrophic flow test case. Figure 13b shows that the normal component of the velocity error is greater on the Africa continent, which is a coarse-grid region. This velocity field error in the VR7 grid is mitigated when we add fourth-order hyperdiffusion based on the alignment index, as in Sect. 4.2.1, in the shallow-water model, as Fig. 13a shows. Again, the addition of fourth-order hyperdiffusion allowed us to have a similar temporal error evolution in both variable- and uniform-resolution SCVT grids, breaking the total energy conservation, which is conserved with precision 10⁻¹⁰ without fourth-order hyperdiffusion and 10⁻⁸ with fourth-order hyperdiffusion.

4.2.3 Barotropic unstable zonal jet

Now, we shall analyze the test case proposed by Galewsky et al. (2004). The initial fields are given by

$$\begin{array}{}\text{(25)}& {\displaystyle }u(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=\left\{\begin{array}{ll}\mathrm{0}& \text{if }\mathit{\varphi }\le {\mathit{\varphi }}_{\mathrm{0}}\\ \frac{u_{\text{max}}}{e_n}\exp \frac{\mathrm{1}}{(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}})(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{1}})}& \text{if }{\mathit{\varphi }}_{\mathrm{0}}<\mathit{\varphi }\le {\mathit{\varphi }}_{\mathrm{1}},\\ \mathrm{0}& \mathit{\varphi }\ge {\mathit{\varphi }}_{\mathrm{1}},\end{array}\right.\text{(26)}& {\displaystyle }v(\mathit{\lambda },\mathit{\varphi },\mathrm{0})=\mathrm{0},\text{(27)}& {\displaystyle }\begin{array}{rl}h(\mathit{\lambda },\mathit{\varphi },\mathrm{0})& =h_{\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{g}}\stackrel{\mathit{\varphi }}{\int }au(\mathit{\lambda },{\mathit{\varphi }}^{\prime },\mathrm{0})\\ & \times \left(f+{\displaystyle \frac{\tan {\mathit{\varphi }}^{\prime }}{a}}u(\mathit{\lambda },{\mathit{\varphi }}^{\prime },\mathrm{0})\right)\mathrm{d}{\mathit{\varphi }}^{\prime },\end{array}\end{array}$$

where we define ${\mathit{\varphi }}_{\mathrm{0}}=-\mathrm{5}$^∘, ${\mathit{\varphi }}_{\mathrm{1}}=-\mathrm{45}$^∘, u_max=80 and $e_n=\exp -\frac{\mathrm{4}}{({\mathit{\varphi }}_{\mathrm{0}}-{\mathit{\varphi }}_{\mathrm{1}}{)}^{\mathrm{2}}}$. These parameters define the jet in the Southern Hemisphere. At last, we add the following perturbation in the height field in order to trigger the instability:

$$\begin{array}{}\text{(28)}& \widehat{h}(\mathit{\lambda },\mathit{\varphi })=\widehat{h}\cos \left(\mathit{\varphi }\right)e^{-{\left(\frac{\mathit{\lambda }}{\mathit{\alpha }}\right)}^{\mathrm{2}}}{e}^{-{\left(\frac{\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{2}}}{\mathit{\beta }}\right)}^{\mathrm{2}}}\end{array}$$

where $\widehat{h}=\mathrm{120}$, $\mathit{\alpha }=\frac{\mathrm{1}}{\mathrm{3}}$, $\mathit{\beta }=\frac{\mathrm{1}}{\mathrm{15}}$ and ${\mathit{\varphi }}_{\mathrm{2}}=-\mathrm{25}$^∘. The bottom topography is set to zero.

Figure 14Barotropic unstable zonal jet for the shallow-water model: potential vorticity for the VR7 grid at (a) day 6 and (b) day 7 without fourth-order hyperdiffusion. Panel (c) shows the potential vorticity at day 7 with fourth-order hyperdiffusion in the shallow-water model for the VR7 grid. In (d) and (e) we depict the potential vorticity at day 7 for the UR7 and UR8 grids, respectively, as references.

We ran this test for 7 d for the VR7, UR7 and UR8 grids. In Fig. 14 we illustrate the potential vorticity. On day 6, in the VR7 grid, numerical noise is generated in a few grid cells on the topography-based refined region, as is shown in Fig. 14a. This noise propagates to other grid cells at day 7 (Fig. 14b) and its magnitude is reduced. These results show that the performance of the VR7 grid is problematic for this test case. For this reason, we ran this test again employing the fourth-order hyperdiffusion based on the alignment index, as in Sect. 4.2.1, in the shallow-water model. As Fig. 14c shows, the numerical noise in the VR7 grid is mitigated when we add numerical hyperdiffusion. In Fig. 14d and e, we depict the potential vorticity at day 7 for the UR7 and UR8 grids, respectively, for which no numerical fourth-order hyperdiffusion was employed. When we compare the potential vorticity at day 7 for the grids VR7 and UR7, both of level 7, with the UR8 grid, we can observe that the potential vorticity in the VR7 grid (Fig. 14c), considering numerical hyperdiffusion, better represents the vortex formed near the Andes region than in the uniform-resolution grid with the same number of cells (Fig. 14d); i.e., the vortex formed near the Andes region in the VR7 grid with numerical fourth-order hyperdiffusion is more similar to the vortex formed in the same area in the UR8 grid (Fig. 14e) compared to the vortex generated in the UR7 grid. This simulation shows that the variable-resolution grid has a good performance in this test case when we add numerical hyperdiffusion in the shallow-water model compared to the uniform grid with the same number of grid cells, showing a computational advantage of the variable-resolution grid. We do, however, highlight the fact that since the test case is dynamically unstable, exact agreement of uniform- and variable-resolution grid solutions is not necessarily expected, but one may notice the similar qualitative representation of the vortex formation in both solutions.

Figure 14b shows a grid-scale checker-boarding pattern in the potential vorticity. TRSK is known for having grid-scale oscillations in the potential vorticity, and some potential stabilization methods have been proposed in the literature to tackle this problem, such as the anticipated potential vorticity method (APVM) (Ringler et al., 2010) and the continuous, linear-upwind stabilized transport scheme (CLUST) (Weller, 2012). We re-run the barotropic unstable jet test case for the VR7 grid using APVM and CLUST; however, our results show that these methods do not help to mitigate the potential vorticity grid-scale oscillations.

4.2.4 Matsuno baroclinic wave

This test case has been proposed by Shamir et al. (2019) and considers analytical wave solutions to the linearized shallow-water equations on the beta plane that approximate analytical solutions to the shallow-water equations on the sphere for low gravity wave speed and are taken as initial and reference solutions for this test case. It is recommended to analyze two types of waves: eastward inertia–gravity waves (EIGWs hereafter) and Rossby waves (RWs hereafter). We follow Shamir et al. (2019) by analyzing both waves in this work. The mean atmosphere depth is set to H=30 m, the wavenumber is defined as k=5 and the mode number is defined as n=1. The EIGW period is T=1.9 d and the RW period is T=18.5 d. We recall that the time frequency is defined by $\mathit{\omega }=\frac{\mathrm{2}\mathit{\pi }}{T}$. For these parameters, the analytical solutions are expected to be waves that propagate in the zonal direction, while the EIGW propagates in the eastward direction and the RW propagates in the westward direction. In order to guarantee that waves pass through the whole Andes topography-based locally refined region, we applied a three-dimensional rotation of 21^∘ in the x−z plane in the VR7 grid. We integrate the shallow-water equations for 100 wave periods for each EIGW and RW, as is suggested in Shamir et al. (2019), for both rotated VR7 and UR7 grids to investigate any local refinement impact on the solution.

Figure 15Matsuno baroclinic wave test case: simulated h^′ for the EIGW in the rotated VR7 grid (a) and the analytical h^′ (b) after 100 EIGW periods (190 d). In (c) and (d) we depict the simulated and analytical h^′ for the Rossby wave in the rotated VR7 grid after 100 RW periods (1850 d).

Download

We define the variable $h^{\prime }=h-H$ in order to analyze the mean atmosphere depth perturbation and for plot reasons. In Fig. 15 we depict the h^′ field obtained for the EIGW (Fig. 15a) and the analytical h^′ (Fig. 15b) in the VR7 grid. The same results for the RW are shown Fig. 15c and d. For both EIGW and RW, the error caused is due to phase error, which is a common error for many numerical schemes, and similar results are obtained for the UR7 grid. Despite the waves propagating along the locally refined region, the obtained solutions preserved the wave shape even for a long time of integration.

Figure 16EIGW latitude–time Hovmöller diagrams for the simulated h^′ in the rotated VR7 (a) and UR7 (b) grids, obtained fixing the longitude at $\mathit{\lambda }=-\mathrm{65}$^∘, and EIGW longitude–time Hovmöller diagrams for the simulated h^′ in the rotated VR7 (c) and UR7 (d) grids obtained fixing the latitude at ϕ=0^∘. T denotes the period of the EIGW (1.9 d). The white dashed line has slope $\frac{k}{\mathit{\omega }}$, where ω and k are the EIGW time frequency and wavenumber, respectively.

Download

Figure 17Same as Fig. 16 for the RW.

Download

We also depict Hovmöller diagrams for h^′ considering the EIGW and RW in Figs. 16 and 17, respectively, for both rotated VR7 and UR7 grids. The latitude–time Hovmöller diagram is obtained by intersecting h^′ at $\mathit{\lambda }=-\mathrm{65}$^∘, and the longitude–time Hovmöller diagram is obtained by intersecting h^′ at ϕ=0^∘. As Fig. 16a and b show, the VR7 and UR7 grids have the same phase error for h^′ in the EIGW test case. From Fig. 16c–d, comparing the black lines with the white line, we can conclude that the obtained EIGW is propagating in agreement with the analytical speed for both grids since the analytical solutions propagate zonally; the longitude–time Hovmöller diagrams are expected to be straight lines with slope $\frac{k}{\mathit{\omega }}$ (Shamir et al., 2019). The same analysis holds for the RW Hovmöller diagrams depicted in Fig. 17. Therefore, the local refinement introduces no impact on the obtained solution for this test case as we observed in the previous test cases.

We remark that these test cases proposed by Shamir et al. (2019), while providing valuable insight into wave propagation properties, seem not to be very sensitive to different numerical methods. We expect that most numerical methods would perform well on these tests. This was the case for our experiments on both uniform and locally refined grids and was also noticed in Brachet and Croisille (2021) for a cubed-sphere-based model.

4.3 Moist shallow-water model

We now focus our analyses of the variable-resolution grids on the moist shallow-water model framework (Zerroukat and Allen, 2015) described in Sect. 3.2. We consider two test cases. The first test case is suggested by Zerroukat and Allen (2015) and is similar to the flow over a mountain test case (Williamson et al., 1992). The second test case is similar to the barotropic unstable jet test case (Galewsky et al., 2004). This test case was analyzed in Ferguson et al. (2019); however, in their work, the rain formed is considered precipitated and is removed from the model. In our analysis, we consider the rain to be advected as a tracer following Zerroukat and Allen (2015). In this section, on all test cases that use variable-resolution grids, we have added the alignment-based hyperdiffusion in the momentum equation (Eq. 13) just as in Sect. 4.1, aiming to mitigate the variable-resolution numerical noise observed in Sect. 4.2. We considered the maximum value of hyperdiffusion as described in Sect. 4.2.1. As we pointed out in Sect. 3.2, we discretize the tracer equations using the TRSK divergence operator, which is a centered scheme. Such centered schemes are known for having inherent grid-scale oscillations that are not necessarily removed using the dissipation mechanism proposed in this work. Despite that, this analysis will allow us to understand how the discrete divergence is affected by grid distortions.

The vapor, cloud and rain variables in the moist shallow-water model must assume only non-negative values. For this purpose, we employed a monotonic filter to ensure that the moisture variables assume only non-negative values. At each time step of the model integration, we apply the monotonic filter for each tracer (vapor, cloud and rain), which is implemented as an iterative process. At each iteration of the filter, the tracer mass is computed for each grid cell. The grid cells with negative tracer mass are set to zero mass and the negative mass is distributed equally in the neighbor cells. This iterative process is repeated until the maximum magnitude of the negative values is small enough. After that, all the remaining mass negative values are set to zero, and the tracer value for a cell is computed as the value of its mass divided by its grid cell area. This process, even though it does not conserve the tracer mass, is shown to avoid negative mass, and the total mass of $h(q_{\mathrm{v}}+q_{\mathrm{c}}+q_{\mathrm{r}})$ relative to variation is less than 10⁻⁵ for the simulations that we shall see in the remaining of this section.

4.3.1 Flow over a mountain

This test case is proposed in Zerroukat and Allen (2015). We only slightly changed the initial conditions by translating the initial fields into 30^∘ longitude in order to ensure cloud and rain formation on the South American continent. The initial height and velocity are defined as in the flow over a mountain test case for the classical shallow-water model described in Sect. 4.2.2. The initial temperature and vapor are respectively given by

$$\begin{array}{}\text{(29)}& {\displaystyle }\begin{array}{rl}\mathit{\theta }(\mathit{\lambda },\mathit{\varphi },\mathrm{0})& =F\left({\mathit{\theta }}^{\mathrm{SP}},(\mathrm{1}-{\mathit{\mu }}_{\mathrm{1}}){\mathit{\theta }}^{\mathrm{EQ}},{\mathit{\theta }}^{\mathrm{NP}},\mathit{\varphi }\right)\\ & +{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}^{\mathrm{EQ}}\cos \mathit{\varphi }\sin \left(\mathit{\lambda }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}\right),\end{array}\text{(30)}& {\displaystyle }{q}_{\mathrm{v}}(\mathit{\lambda },\mathit{\varphi },\mathrm{0})={\mathit{\mu }}_{\mathrm{2}}{q}_{\mathrm{sat}}\left(b,h(\mathit{\lambda },\mathit{\varphi },\mathrm{0}),\mathit{\theta }(\mathit{\lambda },\mathit{\varphi },\mathrm{0})\right),\end{array}$$

where F is the function given by

$$\begin{array}{}\text{(31)}& \begin{array}{rl}F(f_{\mathrm{1}},f_{\mathrm{2}},f_{\mathrm{3}},\mathit{\varphi })& ={\displaystyle \frac{\mathrm{2}}{{\mathit{\pi }}^{\mathrm{2}}}}\left[\mathit{\varphi }\left(\mathit{\varphi }-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)f_{\mathrm{1}}\right.\\ & -\mathrm{2}\mathit{\varphi }\left(\mathit{\varphi }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)\left(\mathit{\varphi }-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)f_{\mathrm{2}}\\ & \left.+\mathit{\varphi }\left(\mathit{\varphi }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)f_{\mathrm{3}}\right].\end{array}\end{array}$$

The parameters are set as ${\mathit{\theta }}^{\mathrm{NP}}=-\mathrm{40}\mathit{\epsilon }$, θ^EQ=40ε, ${\mathit{\theta }}^{\mathrm{SP}}=-\mathrm{20}\mathit{\epsilon }$, $\mathit{\epsilon }=\frac{\mathrm{1}}{\mathrm{300}}$, μ₁=0.05 and μ₂=0.98. The function q_sat is described in Sect. 3.2, Eq. (17). The topography is given by

$$\begin{array}{}\text{(32)}& b(\mathit{\lambda },\mathit{\varphi })=\left\{\begin{array}{ll}\mathrm{2000}\left(\mathrm{1}-\frac{r}{r_{\text{max}}}\right)& \text{if }r\le r_{\text{max}}\\ \mathrm{0}& \text{otherwise.}\end{array}\right.\end{array}$$

Here, $r=\sqrt{(\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{0}}{)}^{\mathrm{2}}+(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}{)}^{\mathrm{2}}}$, $r_{\max }=\frac{\mathit{\pi }}{\mathrm{9}}$, ${\mathit{\lambda }}_{\mathrm{0}}=\frac{\mathrm{2}\mathit{\pi }}{\mathrm{3}}$ and ${\mathit{\varphi }}_{\mathrm{0}}=-\frac{\mathit{\pi }}{\mathrm{6}}$. The initial fields are shown in Fig. 18. The initial cloud and rain are set as zero.

Figure 18Flow over a mountain test case for the moist shallow-water model: initial conditions for the vapor (a) and temperature (b) fields, as well as the considered topography (c).

We ran this test for 30 d on the uniform-resolution SCVT grids using the grids UR6 to UR9 and the VR7 grid. In Fig. 19 we present the results for the cloud field in the south of the continent of South America at day 30 for the grids considered in our simulations. Figure 20 illustrates the same for the rain field. We employed numerical fourth-order hyperdiffusion based on the alignment index in the VR7 grid simulation as in Sect. 4.2.1.

Figure 19Flow over a mountain test case for the moist shallow-water model: cloud field at day 30 for the grids UR6 (a), UR7 (b), UR8 (c), UR9 (d) and VR7 (e) as well as the cell diameters in kilometers for the VR7 grid (f). The UR6, UR7, UR8 and UR9 grids have an average resolution of 120, 60, 30 and 15 km, respectively.

Figure 20Flow over a mountain test case for the moist shallow-water model: rain field at day 30 for the grids UR6 (a), UR7 (b), UR8 (c), UR9 (d) and VR7 (e) as well as the cell diameters in kilometers for the VR7 grid (f). The UR6, UR7, UR8 and UR9 grids have an average resolution of 120, 60, 30 and 15 km, respectively.

From Figs. 19a–d and 20a–d we can notice that both cloud and rain fields are converging in the uniform-resolution SCVT grids. For the variable-resolution SCVT grid with local refinement on the Andes, the cloud field is better represented in the refined region when we compare Fig. 19e and b with Fig. 19d. This illustrates a potential benefit of the locally refined grid, as the resulting clouds are well represented in the refined region while having lower computational cost with respect to the globally high-resolution uniform grid. However, when we look at the rain field in Fig. 20, we can see that the rain in the VR7 grid (Fig. 20e) is more similar to the rain in the UR7 (Fig. 20b) than in the UR8 (Fig. 20c) or UR9 (Fig. 20d) when we look at the topography-based refined region. Notice that the formed rain in the topography-based refined region that is not formed in the UR8 and UR9 grids is present in the UR6 and UR7 grids. This indicates that the low-resolution region of the variable-resolution grid is affecting the observed rain in the high-resolution region.

4.3.2 Barotropic unstable zonal jet

In this test case, we initialize the height and velocity fields as in the unstable barotropic jet for the classical shallow-water model (Sect. 4.2.3). The temperature and vapor are set as in the previous test (Eqs. 29 and 30) with μ₁=0. This test considers the bottom topography to be zero on the whole sphere. A perturbation in the height field is also added, and we consider the longitude to be coordinated and translated using the transformation ${\mathit{\lambda }}^{\prime }=\mathit{\lambda }-\mathrm{0.93}\mathit{\pi }$ only for the perturbation since the other fields do not depend on λ.

Figure 21Unstable barotropic jet test case for the moist shallow-water model. (a–c) Cloud field at day 7 for the UR8, UR9 and VR7 grids, respectively. (d–f) Cloud field at day 8 for the UR8, UR9 and VR7 grids, respectively. (g) Grid cell diameters in kilometers for the VR7 grid. The UR8 and UR9 grids have an average resolution of 30 and 15 km, respectively.

Figure 22Unstable barotropic jet test case for the moist shallow-water model. (a–c) Rain field at day 7 for the UR8, UR9 and VR7 grids, respectively. (d–f) Rain field at day 8 for the UR8, UR9 and VR7 grids, respectively. (g) Grid cell diameters in kilometers for the VR7 grid. The UR8 and UR9 grids have an average resolution of 30 and 15 km, respectively.

We ran this test for 8 d in the UR8, UR9 and VR7 grids, whereby we consider numerical fourth-order hyperdiffusion only for the VR7 grid just as the previous test case. Figures 21 and 22 show cloud and rain formation for these grids at days 7 and 8. The rain at day 7 in the variable-resolution grid (Fig. 22c) has a higher magnitude when compared to the UR8 (Fig. 22a) and UR9 grids (Fig. 22b). Nevertheless, this is expected since the region of the maximum magnitude of the rain is in the transition zone between the low- and high-resolution grid. A similar analysis applies for the cloud field at day 7.

At day 8, the rain is generated on the Andes in the UR8 (Fig. 22d) and UR9 grids (Fig. 22e). However, this rain formation is not properly generated in the Andes region (Fig. 22f), being misplaced westward of the Andes, even though the Andes are better represented in our variable-resolution grid. This illustrates how the local refinement was not able to solve the rain field in this case and influenced the large-scale dynamics. Once again, a similar influence is observed for the cloud field on day 8.