Development and technical paper 16 Nov 2021
Development and technical paper  16 Nov 2021
Topographybased local spherical Voronoi grid refinement on classical and moist shallowwater finitevolume models
 Instituto de Matemática e Estatística da Universidade de São Paulo, Rua do Matão, 1010 – Butantã, São Paulo – SP, 05508090, Brazil
 Instituto de Matemática e Estatística da Universidade de São Paulo, Rua do Matão, 1010 – Butantã, São Paulo – SP, 05508090, Brazil
Correspondence: Luan F. Santos (luan.santos@usp.br)
Hide author detailsCorrespondence: Luan F. Santos (luan.santos@usp.br)
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 highresolution 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 finestructure 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 finitevolume scheme employed in the MPAS dynamical core on this grid considering the nonlinear classic and moist shallowwater 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 shallowwater 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 uniformresolution SCVT grids. Fortunately, the spurious waves originate from smallscale gridrelated numerical errors and can therefore be mitigated using fourthorder hyperdiffusion. We exploit a grid geometrybased 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 variableresolution grid when compared to a respective uniformresolution grid with the same number of cells, while in other cases, grid effects can affect the cloud and rain representation.
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 lowlevel 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 lowlevel 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 highresolution 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 cubedsphere 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 outpropagating waves.
Aiming to run at very high resolutions, many recent global atmospheric models are employing quasiuniform 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 finitevolume method for the shallowwater equations proposed by Thuburn et al. (2009) and Ringler et al. (2010), hereafter named TRSK, was designed to work on arbitrary orthogonal Voronoi Cgrids. 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 variableresolution (VR) grids have no significant wave reflection detected in the transition zone between high and lowresolution 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 uniformresolution 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 MPASOcean 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 uniformresolution grid. Also, in the absence of hyperdiffusion, Zhou et al. (2020) report gridscale oscillations in a variableresolution 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 shallowwater 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 shallowwater model proposed by Zerroukat and Allen (2015) is investigated. This model includes physical processes in the shallowwater model and allows us to simulate rain and clouds. The moist shallowwater 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 shallowwater 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 shallowwater model. Conclusions are drawn in Sect. 5.
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 nearestneighbor 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 cocircular 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 dualcell 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).
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 ℝ^{3} coordinates of a point x∈𝒮^{2} 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}}=(\mathrm{sin}{\mathit{\varphi}}_{\mathrm{c}}\mathrm{cos}{\mathit{\lambda}}_{\mathrm{c}},\mathrm{sin}{\mathit{\varphi}}_{\mathrm{c}}\mathrm{sin}{\mathit{\lambda}}_{\mathrm{c}},\mathrm{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∈𝒮^{2} 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 finegridresolution 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 finitevolume schemes, grid quality plays a major role. On staggered (Ctype, 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 socalled wellcentered 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 wellcentered 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 lowresolution 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.
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 coarsergrid region is approximately 3. We can derive this from Eqs. (4) and (7). Indeed, observe that if x_{i}∈𝒮^{2} 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 uniformresolution 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}∈𝒮^{2} 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.
We generate the variableresolution 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 uniformresolution 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 variableresolution 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 Shallowwater model
As is usual in the development of numerical methods for global atmospheric models, we first evaluate it considering the shallowwater 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 twodimensional; therefore, they are computationally cheaper to run at high resolutions.
We denote a point on the sphere by x and a time instant by t≥0. Then, we consider the nonlinear shallowwater equations on the sphere written in the vectorinvariant 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 twodimensional fluid on the sphere.
The finitevolume 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 Cgrid staggering of the variables. We consider three positions in the Cgrid: primalcell centers (or dualcell vertices), primalcell vertices (or dualcell centers), and the intersection point between a primal edge and dual edge that we call edge points (Fig. 4).
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 righthand 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 fourthorder 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

mass conservation,

the Coriolis term not being an energy source or sink $({u}_{e}{u}_{e}^{\u27c2}=\mathrm{0})$,

total energy being conserved within time error truncation, and

preservation of stationary geostrophic modes on the f sphere.
The TRSK discretization is also known to be very loworder 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 shallowwater model
Even though the shallowwater model retains keys properties of the atmosphere, it lacks the representation of subscale physical processes. Therefore, we consider here an intermediate step of a twodimensional moist shallowwater model, avoiding the computational expense of a threedimensional primitive equation model while being able to account for waterrelated subgrid processes.
We follow the moist shallowwater model proposed by Zerroukat and Allen (2015), wherein moist dynamics are included in the classic shallowwater model. This model is twodimensional, which is computationally cheaper than threedimensional models at high resolutions, and is derived from the Boussinesq approximation of the threedimensional Euler equations. The model can be summarized in the following equations:
where θ is the temperature, $k=\mathrm{1},\mathrm{2},\mathrm{3}$, q^{1}=q_{v} is the water vapor state, q^{2}=q_{c} is the cloud state and q^{3}=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 shallowwater model (Table 2). The source terms of the advection equations define a threestate moist physics, wherein vapor may be converted to clouds, clouds are allowed to evaporate and clouds may be converted to rain. These threestate 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_{0} 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 shallowwater 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 shallowwater 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 shallowwater equations and discretize the righthand 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 highorder 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 (positivitypreserving) filter is employed in our simulations.
This section is dedicated to presenting the results of the TRSK method on the variableresolution grids considering two frameworks: the classical shallowwater model and the moist shallowwater model developed by Zerroukat and Allen (2015). We will consider the variableresolution 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 fourthorder 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 shallowwater model and for the moist shallowwater model in Sect. 4.2 and 4.3, respectively. The L_{2} 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 variableresolution grid triggers unstable modes of the numerical method and spurious waves are generated, destabilizing the numerical solution. While undesirable for the shallowwater 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 fourthorder diffusion (hyperdiffusion), but results with secondorder 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 righthand 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 righthand 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 fourthorder hyperdiffusion ∇^{4} is obtained using Eq. (19) and the relations ∇^{2}=Δ 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_{2} 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_{4}, we may obtain a fourthorder 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 righthand 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 (diameterbased) hyperdiffusion coefficients in Sect. 4.2.1. The diameterbased hyperdiffusion coefficient (K_{4}) 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 _{2}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 variableresolution 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 secondorder accuracy of typical finitevolume operators. Peixoto (2016) shows how the grid alignment can affect the solution of the shallowwater 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 betteraligned cells, so we can be more selective in imposing the diffusion mechanism where it is really needed.
For the alignmentbased hyperdiffusion, we considered a variable coefficient K_{4}, 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 fourthorder hyperdiffusion coefficient allowed is employed in the most illaligned 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. Secondorder 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 Shallowwater 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
The bottom topography is equal to zero. The functions defined above satisfy the shallowwater 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}\phantom{\rule{0.125em}{0ex}}\mathrm{d}}$, where a is the Earth's radius.
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.
To decide on the best hyperdiffusion strategy, we now examine the different formulations: a constant coefficient, a diameterbased variable coefficient and an alignmentbased variable coefficient. For the variableresolution grid of level 6 (VR6), we show in Fig. 6 the height field error after 30 d considering constant and diameterbased hyperdiffusion coefficients considering ${K}_{max}={\mathrm{10}}^{\mathrm{13}}$ m^{4} s^{−1}. 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 diameterbased hyperdiffusion may be well justified from a physical point of view, it may not be ideal for stabilization requirements originating from grid irregularities.
In Fig. 7 we present the relative errors after 30 d for the height and velocity fields in L_{∞} and L_{2} norms considering the VR6 grid with respect to different hyperdiffusion coefficient choices (K_{max}) for the constant, diameterbased and alignmentbased 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 diameterbased coefficient requires higher values of K_{max} to stabilize the error of the height field. The stabilization behavior of the alignmentbased 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 alignmentbased hyperviscosity for further experiments.
For a variableresolution grid of level 7 (VR7), an optimal choice for the hyperdiffusion coefficient may be found at around ${K}_{max}={\mathrm{10}}^{\mathrm{12}}$ m^{4} s^{−1}, 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^{9} to 10^{10} m^{4} s^{−1}). 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 shows the error after 30 d for the height field in the VR7 grid without (Fig. 9a) and with alignmentbased 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 fourthorder 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 illaligned, 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^{4} s^{−1}. For VR6 we will adopt ${K}_{max}={\mathrm{10}}^{\mathrm{13}}$ m^{4} s^{−1}, and for all coarser grids (VR1 to VR5), ${K}_{max}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{13}}$ m^{4} s^{−1} was enough to stabilize the model.
In this test case, the total energy is conserved with precision 10^{−10} for the UR7 grid and 10^{−9} for the VR7 grid without fourthorder hyperdiffusion and with precision 10^{−8} considering fourthorder hyperdiffusion on the VR7 grid. Overall, we can see that the addition of fourthorder hyperdiffusion did not deteriorate the total energy conservation of TRSK, since this dissipation method is a scaleselective damping mechanism (Jablonowski and Williamson, 2011). We point out that the secondorder diffusion ∇^{2} is also able to stabilize the errors in this test case. However, the operator ∇^{2} is not as selective as ∇^{4} and leads to deterioration of total energy conservation. Indeed, we rerun this test for the VR7 grid using ∇^{2} with constant diffusion coefficient equal to 8.2×10^{3} m^{2} s^{−1}, 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^{−4} (not shown).
Furthermore, Fig. 10 shows that adding fourthorder hyperdiffusion in the variableresolution SCVT grid simulations is needed to achieve better convergence in both L_{∞} and L_{2} norms. We can see that the error does not converge to zero for both height and velocity fields, even when we add the fourthorder hyperdiffusion. Nevertheless, lack of convergence is also observed in the uniformresolution 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.
In order to analyze the impact of increasing the resolution of the refined region, we consider the 1 : 6 ratio of finetocoarse cell diameter grid as discussed in Sect. 2.2. In Fig. 11a we illustrate the diameters (in km) of the 1 : 6 variableresolution 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 variableresolution 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 shallowwater model. Their work suggests that ratios between high and lowresolution 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 fourthorder 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
We adopted h_{0}=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 semiLagrangian 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.
We ran this test case for 30 d on the uniform and variableresolution 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 highresolution grid region.
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 coarsegrid region. This velocity field error in the VR7 grid is mitigated when we add fourthorder hyperdiffusion based on the alignment index, as in Sect. 4.2.1, in the shallowwater model, as Fig. 13a shows. Again, the addition of fourthorder hyperdiffusion allowed us to have a similar temporal error evolution in both variable and uniformresolution SCVT grids, breaking the total energy conservation, which is conserved with precision 10^{−10} without fourthorder hyperdiffusion and 10^{−8} with fourthorder 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
where we define ${\mathit{\varphi}}_{\mathrm{0}}=\mathrm{5}$^{∘}, ${\mathit{\varphi}}_{\mathrm{1}}=\mathrm{45}$^{∘}, u_{max}=80 and ${e}_{n}=\mathrm{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:
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.
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 topographybased 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 fourthorder hyperdiffusion based on the alignment index, as in Sect. 4.2.1, in the shallowwater 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 fourthorder 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 uniformresolution grid with the same number of cells (Fig. 14d); i.e., the vortex formed near the Andes region in the VR7 grid with numerical fourthorder 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 variableresolution grid has a good performance in this test case when we add numerical hyperdiffusion in the shallowwater model compared to the uniform grid with the same number of grid cells, showing a computational advantage of the variableresolution grid. We do, however, highlight the fact that since the test case is dynamically unstable, exact agreement of uniform and variableresolution 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 gridscale checkerboarding pattern in the potential vorticity. TRSK is known for having gridscale 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, linearupwind stabilized transport scheme (CLUST) (Weller, 2012). We rerun 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 gridscale 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 shallowwater equations on the beta plane that approximate analytical solutions to the shallowwater 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 topographybased locally refined region, we applied a threedimensional rotation of 21^{∘} in the x−z plane in the VR7 grid. We integrate the shallowwater 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.
We define the variable ${h}^{\prime}=hH$ 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.
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 cubedspherebased model.
4.3 Moist shallowwater model
We now focus our analyses of the variableresolution grids on the moist shallowwater 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 variableresolution grids, we have added the alignmentbased hyperdiffusion in the momentum equation (Eq. 13) just as in Sect. 4.1, aiming to mitigate the variableresolution 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 gridscale 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 shallowwater model must assume only nonnegative values. For this purpose, we employed a monotonic filter to ensure that the moisture variables assume only nonnegative 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^{−5} 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 shallowwater model described in Sect. 4.2.2. The initial temperature and vapor are respectively given by
where F is the function given by
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}}$, μ_{1}=0.05 and μ_{2}=0.98. The function q_{sat} is described in Sect. 3.2, Eq. (17). The topography is given by
Here, $r=\sqrt{(\mathit{\lambda}{\mathit{\lambda}}_{\mathrm{0}}{)}^{\mathrm{2}}+(\mathit{\varphi}{\mathit{\varphi}}_{\mathrm{0}}{)}^{\mathrm{2}}}$, ${r}_{\mathrm{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.
We ran this test for 30 d on the uniformresolution 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 fourthorder hyperdiffusion based on the alignment index in the VR7 grid simulation as in Sect. 4.2.1.
From Figs. 19a–d and 20a–d we can notice that both cloud and rain fields are converging in the uniformresolution SCVT grids. For the variableresolution 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 highresolution 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 topographybased refined region. Notice that the formed rain in the topographybased refined region that is not formed in the UR8 and UR9 grids is present in the UR6 and UR7 grids. This indicates that the lowresolution region of the variableresolution grid is affecting the observed rain in the highresolution 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 shallowwater model (Sect. 4.2.3). The temperature and vapor are set as in the previous test (Eqs. 29 and 30) with μ_{1}=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 λ.
We ran this test for 8 d in the UR8, UR9 and VR7 grids, whereby we consider numerical fourthorder 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 variableresolution 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 highresolution 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 variableresolution grid. This illustrates how the local refinement was not able to solve the rain field in this case and influenced the largescale dynamics. Once again, a similar influence is observed for the cloud field on day 8.
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 shallowwater model and a moist shallowwater 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 shallowwater model. Although the TRSK method is designed for arbitrary orthogonal Cgrids, 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 uniformresolution 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 fourthorder 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 diameterbased hyperdiffusion coefficients, and we could observe that the alignmentbased coefficient has behavior similar to the constant coefficient, while adding hyperdiffusion where it is needed. Therefore, we adopted the alignmentbased 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 lowresolution cell diameters, the errors started to develop at earlier times and the final error increased, but the same quantity of fourthorder 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 highresolution 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, fourthorder hyperdiffusion mitigated the problem, and a vortex near the Andes formation was better represented in the variableresolution 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 fourthorder 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 shallowwater model developed by Zerroukat and Allen (2015) includes moist variables in the classical shallowwater model, wherein conversion of vapor to cloud, cloud to vapor and cloud to rain is possible. We analyzed the TRSK method on the variableresolution 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 variableresolution 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 highresolution solution, whereby we observed a misposition of rain and clouds in the Andes region despite the employment of numerical fourthorder 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 quasiuniformresolution 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 slimshaped 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 quasiuniform 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 loworder. MPAS employs a highorder advection scheme for the tracer equations, which could reduce some of the issues encountered here. Results from highorder finitevolume schemes for tracer transport on unstructured locally refined grids are to be presented in a followup 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 threedimensional 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.
All codes used in this work are openly provided in the GitHub repository at https://github.com/pedrospeixoto/imodel, and the release version relative to results of this paper may be accessed as the Zenodo file https://doi.org/10.5281/zenodo.5667336 (Peixoto et al., 2021).
This paper uses the ETOPO1 data from Amante and Eakins (2009) (https://doi.org/10.7289/V5C8276M). The ETOPO1 data used in this paper may be also accessed as the Zenodo file https://doi.org/10.5281/zenodo.5667336 (Peixoto et al., 2021).
PSP proposed and obtained funding for the project. Both PSP and LFS worked on the method development and analysis. LFS did the implementation of the numerical methods. Both authors wrote and revised the paper.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This research has been supported by the Fundação de Amparo à Pesquisa do Estado de São Paulo (grant nos. 17/251914 and 16/184457).
This paper was edited by James Kelly and reviewed by Nicholas Kevlahan, Darren Engwirda, and one anonymous referee.
Amante, C. and Eakins, B.: ETOPO1 1 ArcMinute Global Relief Model: procedures, data sources and analysis, NOAA National Centers for Environmental Information [data set], https://doi.org/10.7289/V5C8276M, 2009. a, b, c, d
Barros, S. and Garcia, C.: A Global SemiImplicit SemiLagrangian ShallowWater Model on Locally Refined Grids, Mon. Weather Rev., 132, 53–65, https://doi.org/10.1175/15200493(2004)132<0053:AGSSSM>2.0.CO;2, 2004. a
Behrens, J.: Adaptive Atmospheric ModelingKey Techniques in Grid Generation, Data Structures, and Numerical Operations with Applications, Lecture Notes in Computational Science and Engineering, vol. 54, Springer, Berlin, Germany, 2006. a
Brachet, M. and Croisille, J.P.: Spherical shallowwater wave simulation by a cubedsphere finitedifference solver, Q. J. Roy. Meteor. Soc., 147, 786–800, https://doi.org/10.1002/qj.3946, 2021. a
Chou, S. C., Dereczynski, C., Gomes, J. L., Pesquero, J. F., Avila, A. M. H. d., Resende, N. C., Alves, L. F., RuizCárdenas, R., Souza, C. R. d., and Bustamante, J. F. F.: Tenyear seasonal climate reforecasts over South America using the Eta Regional Climate Model, An. Acad. Bras. Cienc., 92, 3, https://doi.org/10.1590/00013765202020181242, 2020. a
Du, Q., Faber, V., and Gunzburger, M.: Centroidal Voronoi Tessellations: Applications and Algorithms, SIAM Rev., 41, 637–676, https://doi.org/10.1137/S0036144599352836, 1999. a, b
Du, Q., Gunzburger, M., and Ju, L.: Constrained Centroidal Voronoi Tessellations for Surfaces, SIAM J. Sci. Comput., 24, 1488–1506, https://doi.org/10.1137/S1064827501391576, 2003. a, b, c, d
Dubos, T., Dubey, S., Tort, M., Mittal, R., Meurdesoif, Y., and Hourdin, F.: DYNAMICO1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility, Geosci. Model Dev., 8, 3131–3150, https://doi.org/10.5194/gmd831312015, 2015. a, b
Engwirda, D.: JIGSAWGEO (1.0): locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geosci. Model Dev., 10, 2117–2140, https://doi.org/10.5194/gmd1021172017, 2017. a
Ferguson, J., Jablonowski, C., and Johansen, H.: Assessing Adaptive Mesh Refinement (AMR) in a Forced ShallowWater Model with Moisture, Mon. Weather Rev., 147, 3673–3692, https://doi.org/10.1175/MWRD180392.1, 2019. a
Figueroa, S., Bonatti, J., Kubota, P., Grell, G., Morrison, H., R. M. Barros, S., Fernandez, J., RamirezGutierrez, E., Siqueira, L., Luzia, G., Silva, J., Silva, J., Pendharkar, J., Capistrano, V., Alvim, D., Enore, D., Diniz, F., Satyamurty, P., Cavalcanti, I., and Panetta, J.: The Brazilian Global Atmospheric Model (BAM): Performance for Tropical Rainfall Forecasting and Sensitivity to Convective Scheme and Horizontal Resolution, Weather Forecast., 31, 1547–1572, https://doi.org/10.1175/WAFD160062.1, 2016. a
Freitas, S. R., Panetta, J., Longo, K. M., Rodrigues, L. F., Moreira, D. S., Rosário, N. E., Silva Dias, P. L., Silva Dias, M. A. F., Souza, E. P., Freitas, E. D., Longo, M., Frassoni, A., Fazenda, A. L., Santos e Silva, C. M., Pavani, C. A. B., Eiras, D., França, D. A., Massaru, D., Silva, F. B., Santos, F. C., Pereira, G., Camponogara, G., Ferrada, G. A., Campos Velho, H. F., Menezes, I., Freire, J. L., Alonso, M. F., Gácita, M. S., Zarzur, M., Fonseca, R. M., Lima, R. S., Siqueira, R. A., Braz, R., Tomita, S., Oliveira, V., and Martins, L. D.: The Brazilian developments on the Regional Atmospheric Modeling System (BRAMS 5.2): an integrated environmental model tuned for tropical areas, Geosci. Model Dev., 10, 189–222, https://doi.org/10.5194/gmd101892017, 2017. a
Galewsky, J., Scott, R., and Polvani, L.: An initialvalue problem for testing numerical models of the global shallowwater equations, Tellus A, 56, 429–440, https://doi.org/10.3402/tellusa.v56i5.14436, 2004. a, b, c, d
Garreaud, R. D.: The Andes climate and weather, Adv. Geosci., 22, 3–11, https://doi.org/10.5194/adgeo2232009, 2009. a
Harris, L., Lin, S.J., and Tu, C.: HighResolution Climate Simulations Using GFDL HiRAM with a Stretched Global Grid, J. Climate, 29, 4293–4314, https://doi.org/10.1175/JCLID150389.1, 2016. a
Hoch, K. E., Petersen, M. R., Brus, S. R., Engwirda, D., Roberts, A. F., Rosa, K. L., and Wolfram, P. J.: MPASOcean Simulation Quality for VariableResolution North American Coastal Meshes, J. Adv. Model Earth Sy., 12, e2019MS001848, https://doi.org/10.1029/2019MS001848, 2020. a
Insel, N., Poulsen, C., and Ehlers, T.: Influence of the Andes Mountains on South American moisture transport, convection, and precipitation, Clim. Dynam., 35, 1477–1492, https://doi.org/10.1007/s0038200906371, 2010. a
Jablonowski, C. and Williamson, D.: The Pros and Cons of Diffusion, Filters and Fixers in Atmospheric General Circulation Models, in: Numerical Techniques for Global Atmospheric Models. Lecture Notes in Computational Science and Engineering, edited by: Lauritzen, P., Jablonowski, C., Taylor, M., and Nair, R., Springer, Berlin, Heidelberg, vol. 80, 381–493, https://doi.org/10.1007/9783642116407_13, 2011. a, b, c
Ju, L., Ringler, T., and Gunzburger, M.: Voronoi Tessellations and Their Application to Climate and Global Modeling, in: Numerical Techniques for Global Atmospheric Models. Lecture Notes in Computational Science and Engineering, edited by: Lauritzen, P., Jablonowski, C., Taylor, M., and Nair, R., Springer, Berlin, Heidelberg, vol. 80, 313–342, https://doi.org/10.1007/9783642116407_10, 2011. a, b, c, d, e
Junquas, C., Li, L., Vera, C., Treut, H., and Takahashi, K.: Influence of South America orography on summertime precipitation in Southeastern South America, Clim. Dynam., 46, 3941–3963, https://doi.org/10.1007/s0038201528148, 2015. a
Kevlahan, N. K.R. and Dubos, T.: WAVETRISK1.0: an adaptive wavelet hydrostatic dynamical core, Geosci. Model Dev., 12, 4901–4921, https://doi.org/10.5194/gmd1249012019, 2019. a
Klemp, J. B.: Damping Characteristics of Horizontal Laplacian Diffusion Filters, Mon. Weather Rev., 145, 4365–4379, https://doi.org/10.1175/MWRD170015.1, 2017. a, b, c
Kramer, M., Heinzeller, D., Hartmann, H., Berg, W., and Steeneveld, G.J.: Assessment of MPAS variable resolution simulations in the greyzone of convection against WRF model results and observations, Clim. Dynam., 55, 253–276, https://doi.org/10.1007/s003820184562z, 2018. a
Lean, H., Clark, P., Dixon, M., Roberts, N., Fitch, A., Forbes, R., and Halliwell, C.: Characteristics of HighResolution Versions of the Met Office Unified Model for Forecasting Convection over the United Kingdom, Mon. Weather Rev., 136, 3408–3424, https://doi.org/10.1175/2008MWR2332.1, 2008. a
Liu, Y. and Yang, T.: Impact of Local Grid Refinements of Spherical Centroidal Voronoi Tessellations for Global Atmospheric Models, Commun. Comput. Phys., 21, 1310–1324, https://doi.org/10.4208/cicp.050815.020916a, 2017. a
Miura, H. and Kimoto, M.: A Comparison of Grid Quality of Optimized Spherical Hexagonal Pentagonal Geodesic Grids, Mon. Weather Rev., 133, 2817–2833, https://doi.org/10.1175/MWR2991.1, 2005. a
Okabe, A., Boots, B., Sugihara, K., and Chiu, S.: Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, vol. 43, 2nd edn., John Wiley and Sons, https://doi.org/10.2307/2687299, 2000. a
Park, S.H., Klemp, J., and Skamarock, W.: A Comparison of Mesh Refinement in the Global MPASA and WRF Models Using an Idealized NormalMode Baroclinic Wave Simulation, Mon. Weather Rev., 142, 3614–3634, https://doi.org/10.1175/MWRD1400004.1, 2014. a
Peixoto, P.: Accuracy analysis of mimetic finite volume operators on geodesic grids and a consistent alternative, J. Comput. Phys., 310, 127–160, https://doi.org/10.1016/j.jcp.2015.12.058, 2016. a, b, c, d
Peixoto, P. and Barros, S. R. M.: Analysis of grid imprinting on geodesic spherical icosahedral grids, J. Comput. Phys., 237, 61–78, https://doi.org/10.1016/j.jcp.2012.11.041, 2013. a, b, c, d, e, f
Peixoto, P. S. and Barros, S. R.: On vector field reconstructions for semiLagrangian transport methods on geodesic staggered grids, J. Comput. Phys., 273, 185–211, 2014. a
Peixoto, P. S., Thuburn, J., and Bell, M. J.: Numerical instabilities of spherical shallowwater models considering small equivalent depths, Q. J. Roy. Meteor. Soc., 144, 156–171, 2018. a
Peixoto, P. S., Santos, L. F., and Granjeiro, J. B.: pedrospeixoto/iModel: iModel v1.0 (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.5667336, 2021. a, b
Rauscher, S. and Ringler, T.: Impact of VariableResolution Meshes on Midlatitude Baroclinic Eddies Using CAMMPASA, Mon. Weather Rev., 142, 4256–4268, https://doi.org/10.1175/MWRD1300366.1, 2014. a
Ringler, T., Thuburn, J., Klemp, J., and Skamarock, W.: A unified approach to energy conservation and potential vorticity dynamics on arbitrarily structured Cgrids, J. Comput. Phys., 229, 3065–3090, https://doi.org/10.1016/j.jcp.2009.12.007, 2010. a, b, c, d, e, f
Ringler, T., Jacobsen, D., Gunzburger, M., Ju, L., Duda, M., and Skamarock, W.: Exploring a Multiresolution Modeling Approach within the ShallowWater Equations, Mon. Weather Rev., 139, 3348–3368, https://doi.org/10.1175/MWRD1005049.1, 2011. a
Shamir, O., Yacoby, I., Ziskin Ziv, S., and Paldor, N.: The Matsuno baroclinic wave test case, Geosci. Model Dev., 12, 2181–2193, https://doi.org/10.5194/gmd1221812019, 2019. a, b, c, d, e, f
Silva, V., Kousky, V., and Higgins, W.: Daily Precipitation Statistics for South America: An Intercomparison Between NCEP Reanalyses and Observations, J. Hydrometeorol., 12, 101–117, https://doi.org/10.1175/2010JHM1303.1, 2011. a
Skamarock, W., Klemp, J., Duda, M., Fowler, L., Park, S.H., and Ringler, T.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and CGrid Staggering, Mon. Weather Rev., 140, 3090–3105, https://doi.org/10.1175/MWRD1100215.1, 2012. a, b, c, d
Skamarock, W. C. and Gassmann, A.: Conservative Transport Schemes for Spherical Geodesic Grids: HighOrder Flux Operators for ODEBased Time Integration, Mon. Weather Rev., 139, 2962–2975, https://doi.org/10.1175/MWRD1005056.1, 2011. a
Staniforth, A. and Thuburn, J.: Horizontal grids for global weather and climate prediction models: A review, Q. J. Roy. Meteor. Soc., 138, 1–26, https://doi.org/10.1002/qj.958, 2012. a
Thuburn, J., Ringler, T., Skamarock, W., and Klemp, J.: Numerical representation of geostrophic modes on arbitrarily structured Cgrids, J. Comput. Phys., 228, 8321–8335, https://doi.org/10.1016/j.jcp.2009.08.006, 2009. a, b
Thuburn, J., Zerroukat, M., Wood, N., and Staniforth, A.: Coupling a mass conserving semi Lagrangian scheme (SLICE) to a semi implicit discretization of the shallow water equations: Minimizing the dependence on a reference atmosphere, Q. J. Roy. Meteor. Soc., 136, 146–154, https://doi.org/10.1002/qj.517, 2010. a
Tian, X.: Evolutions of Errors in the Global Multiresolution Model for Prediction Across Scales – Shallow Water (MPASSW), Q. J. Roy. Meteor. Soc., 147, 382–391, https://doi.org/10.1002/qj.3923, 2020. a
Ullrich, P. A., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R., Reed, K. A., Zarzycki, C. M., Hall, D. M., Dazlich, D., Heikes, R., Konor, C., Randall, D., Dubos, T., Meurdesoif, Y., Chen, X., Harris, L., Kühnlein, C., Lee, V., Qaddouri, A., Girard, C., Giorgetta, M., Reinert, D., Klemp, J., Park, S.H., Skamarock, W., Miura, H., Ohno, T., Yoshida, R., Walko, R., Reinecke, A., and Viner, K.: DCMIP2016: a review of nonhydrostatic dynamical core design and intercomparison of participating models, Geosci. Model Dev., 10, 4477–4509, https://doi.org/10.5194/gmd1044772017, 2017. a
Weller, H.: Controlling the Computational Modes of the Arbitrarily Structured C Grid,, Mon. Weather Rev., 140, 3220–3234, https://doi.org/10.1175/MWRD1100221.1, 2012. a, b
Williamson, D., Drake, J., Hack, J., Jakob, R., and Swarztrauber, P.: A Standard Test Set for Numerical Approximations to the Shallow Water Equations in Spherical Geometry, J. Comput. Phys., 102, 211–224, https://doi.org/10.1016/S00219991(05)800166, 1992. a, b, c, d, e, f, g, h
Zarzycki, C. M., Jablonowski, C., and Taylor, M. A.: Using VariableResolution Meshes to Model Tropical Cyclones in the Community Atmosphere Model, Mon. Weather Rev., 142, 1221–1239, https://doi.org/10.1175/MWRD1300179.1, 2014. a
Zerroukat, M. and Allen, T.: A moist Boussinesq shallow water equations set for testing atmospheric models, J. Comput. Phys., 290, 55–72, https://doi.org/10.1016/j.jcp.2015.02.011, 2015. a, b, c, d, e, f, g, h, i, j, k
Zhou, Y., Zhang, Y., Li, J., Yu, R., and Liu, Z.: Configuration and evaluation of a global unstructured mesh atmospheric model (GRISTA20.9) based on the variableresolution approach, Geosci. Model Dev., 13, 6325–6348, https://doi.org/10.5194/gmd1363252020, 2020. a, b, c