Topography-based local spherical Voronoi grid refinement on classical and moist shallow-water finite-volume models

Locally refined grids for global atmospheric models are attractive since they are expected to provide an alternative to solve local phenomena without the requirement of a global high-resolution uniform grid, whose computational cost may be prohibitive. Spherical centroidal Voronoi tessellation (SCVT), as used in the atmospheric Model for Prediction Across Scales (MPAS), allows a flexible way to build and work with local refinement. In addition, the Andes Range plays a key role in the South American weather, but it is hard to capture its fine-structure dynamics in global models. This paper describes how to generate SCVT grids that are locally refined in South America and that also capture the sharp topography of the Andes Range by defining a density function based on topography and smoothing techniques. We investigate the use of the mimetic finite-volume scheme employed in the MPAS dynamical core on this grid considering the nonlinear classic and moist shallow-water equations on the sphere. We show that the local refinement, even with very smooth transitions from different resolutions, generates spurious numerical inertia–gravity waves that may even numerically destabilize the model. In the moist shallow-water model, wherein physical processes such as precipitation and cloud formation are included, our results show that the local refinement may generate spurious rain that is not observed in uniform-resolution SCVT grids. Fortunately, the spurious waves originate from small-scale grid-related numerical errors and can therefore be mitigated using fourth-order hyperdiffusion. We exploit a grid geometry-based hyperdiffusion that is able to stabilize spurious waves and has very little impact on the total energy conservation. We show that, in some cases, the clouds are better represented in a variable-resolution grid when compared to a respective uniform-resolution grid with the same number of cells, while in other cases, grid effects can affect the cloud and rain representation.

1 Introduction

The Andes, the world's longest mountain range, located along entire South American western coast, act as an obstacle to atmospheric circulation and play a key role in the weather of South America. For instance, the South American low-level jet, formed in the Andes eastern slope, has an important influence on the Brazilian southeast, among other regions, with precipitation transporting moisture from the Amazon Basin (Garreaud, 2009). Furthermore, the Andes are essential for the formation of the South American low-level jet. Insel et al. (2010) and Junquas et al. (2015) show that the remotion of the Andes leads to a reduction of precipitation in the Brazilian southeast. Despite their length, the Andes have a narrow width of less than 200 km almost everywhere, and their elevation is 4000 m on average, but the elevation may suddenly change in a few kilometers.

The Andes' complex topography is hard to capture in a global model due to abrupt changes in elevation and resolution limitations imposed by computational costs. Furthermore, Figueroa et al. (2016) show that both the Brazilian Global Atmospheric Model (BAM) and the Global Forecast System (GFS) produce dry or wet biases in the Amazon Basin, over the Andes and in the Brazilian southeast. This bias is also observed in other models (Silva et al., 2011; Chou et al., 2020). Therefore, aiming to improve the precipitation forecast in those regions, it is desirable to better represent the shape of the Andes Range in the employed model. One alternative to overcome the problem of resolution limitations is to employ regional models (e.g., Freitas et al., 2017). A second alternative is to develop a locally refined global grid that captures the Andes Range well. The second alternative has the advantage of not requiring any artificial boundary conditions, as is usual for regional models, and has been developed for some recent global models (Ullrich et al., 2017).

Locally refined grids for global models are built aiming to solve local phenomena without requiring a global high-resolution grid for which computational cost may be prohibitive. For instance, Barros and Garcia (2004) have proposed a locally refined latitude–longitude grid with a coarse global uniform grid successively refined to a region of interest. Local refinements have also been developed for the cubed-sphere grid. For example, in Harris et al. (2016) a local refinement for the cubed sphere was developed using stretching transformations, whereby the region of interest has a resolution 7 times that in the coarse grid region. On the other hand, Lean et al. (2008) proposed nesting models with resolutions of 12, 4 and 1 km using the UK Met Office operational model at that time, with the 4 and 1 km models using two squares centered in England as the domain. Those works propose feasible ways to have a higher resolution in areas of interest. However, they are usually constrained to a rectangular geometry in the region of interest and impose interpolation constraints between different resolutions, which may affect the dynamics of in- and out-propagating waves.

Aiming to run at very high resolutions, many recent global atmospheric models are employing quasi-uniform geodesic spherical grids (Staniforth and Thuburn, 2012). The geodesic grids usually have more complex cell geometry, but they can allow more flexibility in local grid refinement. The added geometry complexity on these grids, such as those using Voronoi (hexagonal–pentagonal) cells, makes the development of numerical schemes on these grids more difficult than in regular quadrilateral grids. Nevertheless, the finite-volume method for the shallow-water equations proposed by Thuburn et al. (2009) and Ringler et al. (2010), hereafter named TRSK, was designed to work on arbitrary orthogonal Voronoi C-grids. This scheme has many desirable mimetic problems, such as mass and total energy conservation, as well as preservation of stationary geostrophic modes on the f sphere, among others. It is used in the dynamical core of some atmospheric models such as the Model for Prediction Across Scales (MPAS) (Skamarock et al., 2012) developed jointly by the National Center for Atmospheric Research (NCAR) and the Los Alamos Laboratory, as well as the model Icosahedral Dynamical Core (DYNAMICO) (Dubos et al., 2015).

TRSK applies to a variety of grids, in particular to spherical centroidal Voronoi tessellation (SCVT) grids, which are adopted in MPAS and may be thought of as an optimized version of the hexagonal–pentagonal grid generated from an icosahedral grid. The SCVT grids are orthogonal grids wherein each grid cell is a Voronoi diagram and the center of mass with respect to a given density function is the corresponding Voronoi region generator (Ju et al., 2011). This grid has the flexibility of allowing local refinements by carefully defining a density function and employing Lloyd's method (Du et al., 1999, 2003; Ringler et al., 2011). Using a circular refinement region, Park et al. (2014) show that SCVT variable-resolution (VR) grids have no significant wave reflection detected in the transition zone between high- and low-resolution regions, which occurs in the nested Weather Research and Forecasting (WRF) model analyzed in this work. The work of Kramer et al. (2018) investigates three different meteorological events in Europe and concludes that the VR grid used in MPAS global model yields results comparable to the WRF model. Furthermore, the VR grid performed better than the uniform-resolution grid in the organized convection test case and is computationally cheaper. VR grids are also feasible in the ocean component of MPAS. For instance, Hoch et al. (2020) consider grids with local refinement on the North American coast and conclude that even when the grid quality is degraded, no spurious effects are observed in the MPAS-Ocean model simulation results. On the other hand, some studies show that the local refinements in the MPAS grid may negatively impact the results of the model. For instance, Rauscher and Ringler (2014) show that the VR grid generates an eddy kinetic energy maximum in the refined region for the aquaplanet and Held–Suarez test cases. This maximum is caused by the use of the VR grid and is not observed in the uniform-resolution grid. Also, in the absence of hyperdiffusion, Zhou et al. (2020) report grid-scale oscillations in a variable-resolution grid considering baroclinic wave and tropical cyclone test cases. While MPAS allows a model that is suitable to work with grids with local refinement in arbitrary shapes by only defining a density function, wave representation and interaction in these grids are still not well understood.

In this work, we propose generating SCVT grids with local refinement based on Andes topography. Our locally refined grids will be designed to capture the Andes and the South American continent well with smooth transitions to the coarse grid regions, focusing on obtaining a grid that may allow better precipitation forecasts in South America. Also, we will perform a detailed analysis of the intrinsic advantages and problems of using TRSK scheme in such grids in two frameworks that capture many important structures of the atmosphere.

In the first framework, we will consider the nonlinear shallow-water equations on the sphere, some standard tests proposed in the literature (Williamson et al., 1992; Galewsky et al., 2004) and a more recent test proposed by Shamir et al. (2019). In the second framework, the moist shallow-water model proposed by Zerroukat and Allen (2015) is investigated. This model includes physical processes in the shallow-water model and allows us to simulate rain and clouds. The moist shallow-water model will allow us to study the impact of the local refinement on cloud and rain formation.

This paper is outlined as follows. In Sect. 2, we describe the development of the SCVT grids with local refinement on the Andes Mountains. In Sect. 3, the classical and moist shallow-water models are presented, along with their discretization schemes. In Sect. 4 we present the results obtained using TRSK on the developed grids considering the classical and moist shallow-water model. Conclusions are drawn in Sect. 5.

2 Grids

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

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{\Omega }}_i& =\mathit{\{}\mathit{x}\in {\mathcal{S}}^{\mathrm{2}}:d(\mathit{x},{\mathit{x}}_i)<d(\mathit{x},{\mathit{x}}_j),\\ & \forall i=\mathrm{1},\mathrm{\cdots },n,i\ne j\mathit{\}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(2)}& {\mathit{x}}^{\ast }={\displaystyle \frac{{\int }_{\mathrm{\Omega }}\mathit{x}\mathit{\rho }\left(\mathit{x}\right)\mathrm{d}\mathrm{\Omega }\left(\mathit{x}\right)}{{\int }_{\mathrm{\Omega }}\mathit{\rho }\left(\mathit{x}\right)\mathrm{d}\mathrm{\Omega }\left(\mathit{x}\right)}}.\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\mathit{x}}^{\mathrm{c}}={\displaystyle \frac{{\mathit{x}}^{\ast }}{\Vert {\mathit{x}}^{\ast }\Vert }}.\end{array}$$

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)

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{h_i}{h_j}}\approx {\left({\displaystyle \frac{\mathit{\rho }\left({\mathit{x}}_j\right)}{\mathit{\rho }\left({\mathit{x}}_i\right)}}\right)}^{\frac{\mathrm{1}}{\mathrm{4}}}.\end{array}$$

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 1 SCVT grids using a uniform density function for 12 (a), 42 (b) and 162 (c) generators.

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

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

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathit{x}=(\sin \mathit{\varphi }\cos \mathit{\lambda },\sin \mathit{\varphi }\sin \mathit{\lambda },\cos \mathit{\varphi }),\\ & \mathit{\varphi }\in \left[-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right],\mathit{\lambda }\in [-\mathit{\pi },\mathit{\pi }],\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(6)}& s\left(\mathit{x}\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if }d(\mathit{x},{\mathit{x}}_{\mathrm{c}})\le \mathit{\alpha },\\ \frac{\mathit{\alpha }+\mathit{\epsilon }-d(\mathit{x},{\mathit{x}}_{\mathrm{c}})}{\mathit{\epsilon }}& \text{if }\mathit{\alpha }\le d(\mathit{x},{\mathit{x}}_{\mathrm{c}})\le \mathit{\alpha }+\mathit{\epsilon },\\ \mathrm{0}& \text{otherwise.}\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathit{\rho }\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}^{\mathrm{4}}}}+\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}^{\mathrm{4}}}}\right)\left(\mathit{\mu }s\left(\mathit{x}\right)+(\mathrm{1}-\mathit{\mu })b\left(\mathit{x}\right)\right),\\ & \mathit{x}\in {\mathcal{S}}^{\mathrm{2}},\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(8)}& b_{ij}^{k+\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left(b_{i-\mathrm{1},j}^k+b_{i+\mathrm{1},j}^k+b_{i,j-\mathrm{1}}^k+b_{i,j+\mathrm{1}}^k\right),\end{array}$$

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 3 SCVT 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 1 Minimum 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.

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 Model description

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 2 Variables and parameters for the shallow-water model.

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):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}=-qh{\mathit{u}}^{⟂}-\mathrm{\nabla }B,\text{(10)}& {\displaystyle \frac{\partial h}{\partial t}}=-\mathrm{\nabla }\cdot \left(h\mathit{u}\right),\end{array}$$

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 4 Variable positions for a C-grid staggering.

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 3 Grid notation following Ringler et al. (2010).

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

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial u_e}{\partial t}}=-{\left[qh{\mathit{u}}^{⟂}\right]}_e-[\mathrm{\nabla }B{]}_e,\text{(12)}& {\displaystyle \frac{\partial h_i}{\partial t}}=-[\mathrm{\nabla }\cdot (h\mathit{u}){]}_i,\end{array}$$

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 4 Variables and parameters for the moist shallow-water model (Zerroukat and Allen, 2015).

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:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}=-qh{\mathit{u}}^{⟂}-\mathrm{\nabla }B+S_u,\text{(14)}& {\displaystyle \frac{\partial h}{\partial t}}=-\mathrm{\nabla }\cdot \left(h\mathit{u}\right),\text{(15)}& {\displaystyle \frac{\partial h\mathit{\theta }}{\partial t}}=-\mathrm{\nabla }\cdot \left(h\mathit{\theta }\mathit{u}\right)+h{S}_{\mathit{\theta }},\text{(16)}& {\displaystyle \frac{\partial h{q}^k}{\partial t}}=-\mathrm{\nabla }\cdot \left(h{q}^k\mathit{u}\right)+h{S}_q^k,\end{array}$$

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

$$\begin{array}{}\text{(17)}& q_{\mathrm{sat}}\left(\mathit{\theta }\right)={\displaystyle \frac{q_{\mathrm{0}}}{g(h+b)}}\exp \left(\mathrm{20}\mathit{\theta }\right),\end{array}$$

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 Numerical results

This section is dedicated to presenting the results of the TRSK method on the variable-resolution grids considering two frameworks: the classical shallow-water model and the moist shallow-water model developed by Zerroukat and Allen (2015). We will consider the variable-resolution grids developed in Sect. 2.2 and the uniform SCVT grids to compare the impact of introducing a local refinement. We also consider the addition of variable numerical fourth-order hyperdiffusion on the locally refined grid tests. It will be clear whether we are adding hyperdiffusion or not. We give details about the hyperdiffusion employed in the numerical experiments in Sect. 4.1, and we show the results for the shallow-water model and for the moist shallow-water model in Sect. 4.2 and 4.3, respectively. The L₂ error will be referred to as the root mean square error (RMSE hereafter) and the L_∞ error as the maximum error.

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:

$$\begin{array}{}\text{(18)}& \mathrm{\Delta }\mathit{u}=\mathrm{\nabla }(\mathrm{\nabla }\cdot \mathit{u})-\mathrm{\nabla }\times (\mathrm{\nabla }\times \mathit{u}).\end{array}$$

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

$$\begin{array}{}\text{(19)}& {\left[\mathrm{\Delta }\mathit{u}\right]}_e={\left[\mathrm{\nabla }(\mathrm{\nabla }\cdot \mathit{u})\right]}_e-{\left[\mathrm{\nabla }\times (\mathrm{\nabla }\times \mathit{u})\right]}_e.\end{array}$$

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

$$\begin{array}{}\text{(20)}& (\mathrm{\nabla }\cdot K_{\mathrm{2}}\mathrm{\nabla })\mathit{u}\to \mathrm{\nabla }(K_{\mathrm{2}}\mathrm{\nabla }\cdot \mathit{u})-\mathrm{\nabla }\times (K_{\mathrm{2}}\mathrm{\nabla }\times \mathit{u}),\end{array}$$

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 5 Steady 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.

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 7 Steady 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.

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 8 Hyperdiffusion 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 9 Steady 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 10 Steady 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.

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 12 Flow 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 13 Flow 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 14 Barotropic 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 15 Matsuno 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).

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 16 EIGW 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.

Figure 17 Same as Fig. 16 for the RW.

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 18 Flow 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 19 Flow 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 20 Flow 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 21 Unstable 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 22 Unstable 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.

5 Conclusions

In this work, we developed SCVT grids with local refinement based on topography that refines and captures the Andes Range well and is smoothly transitioned to a circular grid region centered in South America. This circular region is again smoothly transitioned to a coarse global uniform grid. We used smoothing techniques in the Andes topography data in order to generate these grids using Lloyd's method. These grids were developed in the hope to provide a better basis for weather forecasting in South America, aiming at reduced computational cost. We evaluate the TRSK method on these grids, considering the classical shallow-water model and a moist shallow-water model proposed by Zerroukat and Allen (2015), discussing the pros and cons of the developed grids in this work.

We began our analysis with the classical shallow-water model. Although the TRSK method is designed for arbitrary orthogonal C-grids, in particular SCVT grids, we showed some of its deficiencies when considering the SCVT grid with local refinement based on topography. For instance, in test case 2 from Williamson et al. (1992), the error in the normal velocity component grows from the beginning of the model integration, triggering spurious gravity waves in the height field. This test case is known for being stable for a long time of integration in the uniform-resolution SCVT grids. Therefore, our results show a negative impact on the local refinement. However, we could mitigate the error evolution of the velocity field by adding fourth-order hyperdiffusion with a coefficient based on the alignment index. The maximum of the hyperdiffusion coefficient was chosen by looking for the smaller maximum coefficient that stabilized the errors considering test case 2. We also investigated constant and diameter-based hyperdiffusion coefficients, and we could observe that the alignment-based coefficient has behavior similar to the constant coefficient, while adding hyperdiffusion where it is needed. Therefore, we adopted the alignment-based hyperdiffusion coefficient in the further simulations. After adding hyperdiffusion, spurious gravity waves were not generated in the height field and we obtained a similar evolution of the errors with time as the ones obtained in the uniform SCVT grids. Also, we observed that when we increased the ratio between high- and low-resolution cell diameters, the errors started to develop at earlier times and the final error increased, but the same quantity of fourth-order hyperdiffusion was enough to stabilize the error. For test case 5 from Williamson et al. (1992), we replace their mountain considering smooth Andes topography data. This test case showed that the height errors were concentrated over the Andes and generated gravity waves. This result is undesirable since the height errors were greater in the high-resolution region. We could notice that the error in the normal component of the velocity started to grow significantly after 23 d, but this was mitigated with the same amount of hyperdiffusion. For the unstable barotropic jet test case from Galewsky et al. (2004), we observed a lot of numerical noise in the potential vorticity on the Andes refined region. Again, fourth-order hyperdiffusion mitigated the problem, and a vortex near the Andes formation was better represented in the variable-resolution grid in comparison with a uniform grid considering the same number of cells. From these results, we can conclude that the local refinement introduced errors in the velocity field that generate spurious gravity waves, which are visible in the height field. A fourth-order hyperdiffusion is necessary to stabilize the velocity field error. These results are consistent with Zhou et al. (2020), wherein the authors show that the employment of horizontal hyperdiffusion removes numerical noise in the transition zone of a VR SCVT grid.

The moist shallow-water model developed by Zerroukat and Allen (2015) includes moist variables in the classical shallow-water model, wherein conversion of vapor to cloud, cloud to vapor and cloud to rain is possible. We analyzed the TRSK method on the variable-resolution grids considering two test cases. In the first test, similar to test case 5 from Williamson et al. (1992), we observed that the cloud field after 30 d was better represented in the variable-resolution grid of level 7 than in the uniform grid of the same level. However, for the cloud field, we notice spurious rain in the transition zone. Furthermore, spurious rain was also detectable over the Andes refined region, which we address as interference of the coarse grid. For the unstable barotropic test case with moist variables, we notice that both cloud and rain fields were not properly represented in the Andes refined region after 8 d in comparison with a high-resolution solution, whereby we observed a misposition of rain and clouds in the Andes region despite the employment of numerical fourth-order hyperdiffusion.

Our work shows that, despite its flexibility, working with general local refinements in SCVT grids using TRSK, which is also possible in the MPAS global model, requires some care. The main issues are related to spurious numerical waves that may trigger numerical instabilities, requiring numerical dissipation mechanisms to solve this problem. The requirement of dissipation is very common in most general circulation models, and TRSK has so far shown exceptionally good stability properties on quasi-uniform-resolution grids. On locally refined grids, our results show that, while requiring dissipation, the imposed hyperdiffusion had very little impact on the total energy conservation of TRSK. However, the grid may influence rain and cloud formation in the refined region with the TRSK advection scheme used here for the tracers. Overall, our results indicate that locally refined grids can be effective for short periods of time (a few days) but require additional caution for longer periods.

This choice of method for local grid refinement is one among several possibilities (Behrens, 2006). Nested local grid refinement, while requiring some additional computational effort on transition zones, can preserve more regularly shaped cells and is therefore less prone to such spurious modes encountered here. But it can be challenging to build a nesting grid to capture such a slim-shaped featured as the Andes Range. Static adaptive mesh refinement (AMR) could also be a potentially adequate solution to better represent the local feature of interest. A way to preserve the quality of cells observed in quasi-uniform grids, using static AMR to capture features of interest, is by using wavelet representation, as in Kevlahan and Dubos (2019). Overall, the most adequate way to represent the Andes Range for atmospheric modeling purposes still seems to be an open issue.

Additionally, we remark that our tracer equations were discretized here considering the same scheme as used for the continuity equation for the TRSK scheme, which is very low-order. MPAS employs a high-order advection scheme for the tracer equations, which could reduce some of the issues encountered here. Results from high-order finite-volume schemes for tracer transport on unstructured locally refined grids are to be presented in a follow-up work.

Finally, the grids developed in this work, which are available in the open repository https://github.com/pedrospeixoto/imodel (last access: 11 November 2021), can be applied directly in the MPAS model or other models that use unstructured grids, and their impact may be investigated considering the full three-dimensional model, along similar lines as Zhou et al. (2020). From one side, we expect that the unstable modes could be challenging to mitigate due to other already existing unstable modes (e.g., Peixoto et al., 2018). However, since such 3D models already employ a relatively large amount of diffusion or hyperdiffusion (for numerical or physical reasons), most of the issues encountered here could be adequately controlled in 3D.