A high-order conservative collocation scheme and its application to global shallow water equations

Introduction Conclusions References Tables Figures


Introduction
A recent trend in developing global models for atmospheric and oceanic general circulations is the increasing use of the high-order schemes that make use of local reconstructions and have the so-called spectral convergence.Among many others are those reported in Giraldo et al. (2002), Thomas and Loft (2005), Giraldo and Warburton (2005), Nair et al. (2005a, b), Taylor and Fournier (2010) and Blaise and St-Cyr (2012).Two major advantages that make these models attractive are (1) they can reach the targeted numerical accuracy more quickly by increasing the number of degrees of freedom (DOFs) (or unknowns), and (2) they can be more computationally intensive with respect to the data communications in parallel processing (Dennis et al., 2012).
The discontinuous Galerkin (DG) (Cockburn et al., 2000;Hesthaven and Warburton, 2008) and spectral element (SE) (Patera, 1984;Karniadakis and Sherwin, 2005) methods are the widely used frameworks in this context.A more general formulation, the so-called flux reconstruction (FR), was presented in Huynh (2007) which covers a wide spectrum of nodal type schemes, including the DG and SE as the special cases.A FR scheme solves the values at the solution points located within each grid element, and the volume-integrated values, which are the weighted summation of the solutions, can be numerically conserved.We recently proposed a class of local high-order schemes, named multi-moment schemes, which were used to develop the accurate shallow-water models on different spherical grids (Chen and Xiao, 2008;Li et al., 2008;Ii and Xiao, 2010;Chen et al., 2014b).By introducing a multi-moment concept, we showed in Xiao et al. (2013) that the flux reconstruction can be implemented in a more flexible way, and other new schemes can be generated by properly choosing different types of constraint conditions.
In this paper, we introduce a new scheme which is different from the existing nodal DG and SE methods under the FR framework.The scheme, the so-called Gauss-Legendrepoint-based conservative collocation (GLPCC) method, is a kind of collocation method that solves the governing equations of differential form at the solution points, and is very simple and easy to follow.The Fourier analysis and the numerical tests show that the present scheme has the same super convergence property as the DG method.A global shallowwater equation (SWE) model has been developed by implementing the three-point GLPCC scheme on a cubed-sphere grid.The model has been verified by the benchmark tests.The numerical results show the fifth-order accuracy of the present global SWE model.All the numerical outputs look favorably comparable to other existing methods.
The rest of this paper is organized as follows.In Sect.2, the numerical formulations in a one-dimensional case are described in detail.The extension of the proposed scheme to a global shallow-water model on a cubed-sphere grid is discussed in Sect.3. In Sect.4, several widely used benchmark tests are solved by the proposed model to verify its performance in comparison with other existing models.Finally, the Conclusion is given in Sect. 5.

Scheme in one-dimensional scalar case
The first-order scalar hyperbolic conservation law in one dimension is solved in this subsection: where q is a dependent variable and f a flux function.The computational domain, x ∈ [x l , x r ], is divided into I elements with the grid spacing of The computational variables (unknowns) are defined at several solution points within each element, e.g., within element C i the point values, q im (m = 1, 2, . .., M), are defined at the solution points (x im ).High-order schemes can be built by increasing the number of the solution points.In this paper, we describe the GLPCC scheme that has three solution points for each grid element (M = 3).The configuration of local degrees of freedom is shown in Fig. 1 by the hollow circles.To achieve the best accuracy, the DOFs are arranged at Gauss-Legendre points in this study: where x i is the center of the element The unknowns are updated by applying the differentialform governing equations (Eq. 1) at solution points as As a result, the key task left is to evaluate the derivatives of the flux function, which is realized by reconstructing the piecewise polynomial for flux function, F i (x), over each element.Once the reconstructed flux function is obtained, the derivative of flux function is approximated by In Huynh (2007), FR is formulated by two correction functions which assure the continuity at the two cell boundaries and collocate with the so-called primary Lagrange reconstruction at their zero points.Therefore, the existing nodal type schemes can be recast under the FR framework with different correction functions.In Xiao et al. (2013), a more general FR framework was proposed by introducing multimoment constraint conditions including nodal values, firstorder derivatives and even second-order derivatives to determine the flux reconstruction.Here, we will develop a new method to reconstruct the flux function, which is more straightforward and simpler compared with the methods discussed in either Huynh (2007) or Xiao et al. (2013).
We assume that the reconstructed flux function over the ith element, F i (x), has the form of where the coefficients, c i0 , c i1 , . .., c i4 , are determined by a collocation method, which meets five constraint conditions specified at five constraint points (shown in Fig. 1 by the solid circles) as where f i± 1 2 are the values of flux function at the cell boundaries.
In Eq. ( 6), f (q im ) are calculated by three known DOFs at solution points.The values of flux function at the boundaries are obtained by solving the Riemann problems with the values of dependent variables interpolated separately from two adjacent elements.Considering the interface at x i− 1 2 , we get two values of flux function from elements C i−1 and C i as where Q i (x) is a spatial reconstruction for the dependent variable based on local DOFs, having the form of where the Lagrange basis function at the boundary is obtained by an approximate Riemann solver as where a = f q avg i− 1 2 with f (q) = ∂f (q) ∂q being the characteristic speed.A simple averaging q is used in the present paper.
Based on the Riemann solver at cell boundaries, the proposed scheme is essentially an upwind type method.As a result, the inherent numerical dissipation is included and stabilizes the numerical solutions.We did not use any extra artificial viscosity in the shallow-water model for the numerical tests presented in the paper.
It is easy to show that the proposed scheme is conservative in terms of the volume-integrated average of each element: where the weights w im are obtained by integrating the Lagrange basis function as and are exactly the same as those in Gaussian quadrature of degree 5.
A direct proof of this observation is obtained by integrating Eq. (3) over the grid element, yielding the following conservative formulation: where x i q i is the total mass within the element C i .
With the above spatial discretization, the Runge-Kutta method is used to solve the following semi-discrete equation (ODE): where D represents the spatial discretization and q * is the dependent variables known at time t = t * .A fifth-order Runge-Kutta scheme (Fehlberg, 1958) is adopted in the numerical tests to examine the convergence rate: where In other cases, a third-order scheme (Shu, 1988) is adopted to reduce the computational cost, which does not noticeably degrade the numerical accuracy since the truncation errors of the spatial discretization are usually dominant.It is written as where

Spectral analysis and convergence test
We conduct the spectral analysis (Huynh, 2007;Xiao et al., 2013) to theoretically study the performance of the GLPCC scheme by considering the following linear equation This linear equation is discretized on an uniform grid with x = 1.Since the advection speed is positive, the spatial discretization for the three DOFs defined in element C i involves the six DOFs within elements C i and C i−1 and can be written as the following linear combination as where b ms and b ms are the coefficients for the DOFs within elements C i−1 and C i , respectively, which can be obtained by applying the proposed scheme to governing equation Eq. ( 18) in element C i as and with the parameter .
With a wave solution of q (x, t) = e Iω(x+t) (I = √ −1), we have Above spatial discretization can be simplified as (B ms q is ) and Considering the all of DOFs in element C i , a matrix-form spatial discretization formulation is obtained as where q i = [q i1 , q i2 , q i3 ] T and the components of the 3 × 3 matrix B are coefficients B ms (m = 1 to 3, s = 1 to 3).With the wave solution, the exact expression for the spatial discretization of Eq. ( 18) is The numerical property of the proposed scheme can be examined by analyzing the eigenvalues of matrix B in Eq. ( 24).Truncation errors of the spatial discretization are computed by comparing the principal eigenvalues of matrix B and its exact solution −Iω, and the convergence rate can be approximately estimated by the errors at two different wave numbers.The results are shown in Table 1 and the fifthorder accuracy is achieved.The spectrum of B is shown in Fig. 2. A scheme achieves better numerical performance when the hollow circles become closer to the imaginary axis; furthermore, the maximum of spectral radius determines the largest available Courant-Friedrichs-Lewy (CFL) number, i.e., a larger spectral radius corresponding to a smaller available CFL number.Numerical dispersion and dissipation relations dominated by the principal eigenvalues are shown in Fig. 3. Numerical properties of several schemes were analyzed in Xiao et al. (2013), shown in their Fig. 1 for spectra and Fig. 2 for numerical dispassion and dispersion relations.We conduct a comparison between DG3 (three-point discontinuous Galerkin scheme; Huynh, 2007), MCV5 (fifthorder multi-moment constrained finite volume scheme; Ii and Xiao, 2009) and the proposed scheme since these three schemes have the fifth-order accuracy and can be derived by FR framework using different constraint conditions for spatial reconstruction of flux functions.As detailed in Huynh (2007), the DG3 scheme uses the Radau polynomial as the correction functions to derive the flux reconstruction which assure the continuity of the numerical fluxes computed from Riemann solvers at the cell boundaries.The MCV5 scheme can be derived by a general framework for flux reconstruction using multi-moments proposed in Xiao et al. (2013)  uses constraint conditions on the point values, first-and second-order derivatives of flux functions at the cell boundaries where Riemann solvers in terms of derivatives of the flux function are required.Compared with the DG3 scheme, the proposed scheme is easier to be implemented and thus has less computational overheads.Though the MCV5 scheme gives better spectra (eigenvalues are closer to imaginary) than the DG3 scheme and the present scheme, it adopts more DOFs under the same grid spacing, i.e., 4I + 1 DOFs for MCV5 and 3I DOFs for DG3 and the present scheme, where I is the total number of elements.Both MCV5 and the present scheme show slightly higher numerical frequency in the high wave number regime, which is commonly observed in other spectral-convergence schemes, such as DG.Considering the results of the spectral analysis, the proposed scheme is a very competitive framework to build high-order schemes compared with existing advanced methods.Advection of a smooth sine wave is then computed by the GLPCC scheme on a series of refined uniform grids to numerically checking the converge rate.The test case is specified by solving Eq. ( 18) with initial condition q(x, 0) = sin(2π x) and periodical boundary condition over x ∈ [0, 1].
A CFL number of 0.1 is adopted in this example.Normalized l 1 , l 2 and l ∞ errors and corresponding convergence rate are given in Table 2. Again, the fifth-order convergence is obtained, which agrees with the conclusion in the above spectral analysis.

Extension to system of equations
The proposed scheme is then extended to a hyperbolic system with L equations in one dimension, which is written as where q is the vector of dependent variables and f the vector of flux functions.Above formulations can be directly applied to each equation of the hyperbolic system, except that the Riemann problem, which is required at the cell boundaries between different elements to determine the values of flux functions, is solved for a coupled system of equations.
For a hyperbolic system of equations, the approximate Riemann solver used at interface x i− 1 2 is obtained by rewriting Eq. ( 9) as where the vectors by applying the formulations designed for scalar case to each component of the vector.In this paper, we use a simple approximate Riemann solver, the local Lax-Friedrichs (LLF) solver, where a is reduced to a positive real number as where λ l (l = 1 to L) are eigenvalues of matrix A q with A (q) = ∂f (q) ∂q and q 3 Global shallow-water model on cubed-sphere grid

Cubed-sphere grid
The cubed-sphere grid (Sadourny, 1972), shown in Fig. 4, is obtained by projecting an inscribed cube onto a sphere.As a result, the surface of a sphere is divided into six identical patches and six identical curvilinear coordinates are then constructed.Two kinds of projections are adopted to construct the local curvilinear coordinates, i.e., gnomonic and conformal projections (Rancic et al., 1996).Considering the analytic projection relations and more uniform grid spacing, the equiangular gnomonic projection is adopted in the present study.For the transformation laws and the projection relations, one can refer to Nair et al. (2005a, b) for details.Furthermore, a side effect of this choice is that the discontinuous coordinates are found along the boundary edges between adjacent patches.In Chen and Xiao (2008), we have shown that the compact stencils for the spatial reconstructions through using local DOFs are beneficial to suppress the extra numerical errors due to the discontinuous coordinates.

Global shallow-water model
The local curvilinear coordinate system (ξ, η) is shown in Fig. 5, where P is a point on sphere surface, and P is corresponding point on the cube surface through a gnomonic projection.λ and θ represent the longitude and latitude.α and β are central angles spanning from − π 4 to π 4 for each patch.Local coordinates are defined by ξ = Rα and η = Rβ where R is the radius of the Earth.
To build a high-order global model, the governing equations are rewritten onto the general curvilinear coordinates.As a result, the numerical schemes developed for Cartesian grid are straightforwardly applied in the computational space.The shallow-water equations are recast on each spherical patch in flux form as where dependent variables are q = √ Gh, u, v T with water depth h, covariant velocity vector (u, v) and Jacobian of transformation T in η direction with gravitational acceleration g, height of the bottom mountain h s and contravariant velocity vector ( u, v); source term where Here, taking √ Gh as the model variable assures the global conservation of total mass, and the total height is used in the flux term.Consequently, the proposed model can easily deal with the topographic source term in a balanced way (Xing and Shu, 2005).
The numerical formulations for a two-dimensional scheme are easily obtained under the present framework by implementing the one-dimensional GLPCC formulations in ξ and η directions respectively as are discretized along the grid lines in ξ and η directions.
We describe the numerical procedure in ξ direction here as follows.In η direction, similar procedure is adopted for spatial discretization by simply exchanging e and ξ with f and η.Considering three DOFs, i.e., q ij 1nk , q ij 2nk and q ij 3nk , along the nth row (n = 1 to 3) of element on patch k (defined at solution points denoted by the hollow circles in Fig. 6), we have the task to discretize the following equations: As in a one-dimensional case, a fourth-order polynomial E ij nk (x) is built for spatial reconstructions of flux functions e to calculate the derivative of e with regard to ξ as where E (ξ ) can be obtained by applying the constraint conditions at five constraint points (solid circles in Fig. 6) along the nth row of element C ij k , which are pointwise values of flux functions e including three from DOFs directly and other two by solving Riemann problems along the nth rows of the adjacent elements.
The LLF approximate Riemann solver is adopted.It means that the parameter a in Eq. ( 27) reads a = | u| + G 11 gh.Details of solving the Riemann problem in a global shallowwater model using governing equations Eq. ( 29) can be referred to in Nair et al. (2005b).
How to set up the boundary conditions along the twelve patch boundaries is a key problem to construct a global model on cubed-sphere grid.With enough information from the adjacent patch, above numerical formulations can be applied on each patch independently.In the present study, the values of dependent variables are required to be interpolated from the grid lines in the adjacent patch, for example, as shown in Fig. 7 for the boundary edge between patch 1 and patch 4. When we solve the Riemann problem is obtained by interpolation along the grid line P P 1 .Whereas, need to be interpolated from the DOFs defined along grid line P 4 P on patch 4. Since the coordinates on patch 1 and patch 4 are discontinuous at point P , the values of the covariant velocity vector on the coordinate system on patch 4 should be projected to coordinate system on patch 1 and the values of the scalar can be adopted directly.In comparison with our previous study (Chen and Xiao, 2008), in present the study we solve the Riemann problem at patch boundary only in the direction perpendicular to the edge.The parameter a in Eq. ( 27) is determined by the contravariant velocity component perpendicular to the edge and the water depth, which is exactly the same in two adjacent coordinate systems, since the water depth is a scalar independent of the coordinate system and the basis vector perpendicular to the edge is continuous between adjacent patches.As a result, solving the Riemann problem obtains the same result wherever the numerical procedure is conducted on patch 1 or patch 4. So, no additional corrections are required and the global conservation is guaranteed automatically.

Numerical tests
Representative benchmark tests, three from Williamson's standard test cases (Williamson et al., 1992) and one introduced in Galewsky et al. (2004)  water model.All measurements of errors are defined following Williamson et al. (1992).

Williamson's standard case 2: steady-state geostrophic flow
A balanced initial condition is specified in this case by using a height field as where gh 0 = 2.94 × 10 4 , u 0 = 2π R/ (12 days) and the parameter γ represents the angle between the rotation axis and polar axis of the Earth, and a velocity field (velocity components in longitude-latitude grid u λ and u θ ) as As a result, both height and velocity fields should keep unchanging during integration.Additionally, the height field in this test case is considerably smooth.Thus, we run this test on a series of refined grids to check the convergence rate of GLPCC global model.The results of l 1 , l 2 and l ∞ errors and convergence rates are given in Table 3.After extending the proposed high-order scheme to the spheric geometry through the application of the cubed-sphere grid, the original fifthorder accuracy as shown in one-dimensional simulations and spectral analysis is preserved in this test.Numerical results of height fields and absolute errors are shown in Fig. 8 for tests on grid G 12 , which means there are 12 elements in both ξ and η directions on every patch, in the different flow directions, i.e., γ = 0 and γ = π 4 .Compared with our former global model on a cubed sphere, the present model is more accurate in this test.On grid G 20 (240 DOFs along the equator), the normalized errors are l 1 = 1.278 × 10 −7 , l 2 = 2.008 × 10 −7 and l ∞ = 8.045 × 10 −7 , which are almost 1 order of magnitude smaller than those on grid 32×32×6 (with similar number of DOFs; 256 DOFs along the equator) in Chen and Xiao (2008).The influence of patch boundaries on the numerical results can be found in the plots of the absolute errors.The distributions of absolute errors can reflect the locations of patch boundaries, especially in the flow with γ = 0.

Williamson's standard case 5: zonal flow over an isolated mountain
The total height and velocity field in this case is same as the above case 2 with γ = 0, except h 0 = 5960 m and u 0 = 20 m s −1 .A bottom mountain is specified as where . This test is adopted to check the performance of a shallowwater model to deal with a topographic source term.We run this test on a series of refined grids G 6 , G 12 , G 24 and G 48 .Numerical results of height fields are shown in Fig. 9 for the total height field of the test on grids G 12 at day 5, 10 and 15, which agree well with the spectral transform solutions on T213 grid (Jakob-Chien et al., 1995).Furthermore, the oscillations occurring at the boundary of the bottom mountain observed in spectral transform solutions are completely removed through a numerical treatment which balances the numerical flux and topographic source term (Chen and Xiao, 2008).The numerical results on finer grids are not depicted here since they are visibly identical to the results shown in Fig. 9. Present model assures the rigorous conservation of the total mass as shown in Fig. 10.The conservation errors of total energy and enstrophy are of particular interest for evaluating the numerical dissipation of the model.As shown in Fig. 11, the conservation errors for total energy (left panel) and potential enstrophy (right panel) of tests on a series of refined grids are checked.As in the above case, to compare with our former fourth-order model this test case is checked on grid G 20 having the similar DOFs as the former 32×32×6 grid.The conservation errors are −9.288× 10 −7 for total energy and −1.388 × 10 −5 for potential enstrophy and much smaller than those by fourth-order model in Chen and Xiao (2008).99 00 99 00 9 9 0 0 9 9 0 0 99 00 99 00 9900 990 0 1 0 1 0 0 1 0 1 0 0 1 0 1 0 0 10100 1 0 1 0 0 1 0 1 0 0 10 10 0 10 10 0 10 30 0 1 0 3 0 0 1 0 3 0 0 1 0 3 0 0 1 0 3 0 0  high-order schemes are always preferred to better capture the evolution of small scales.The spectral transform solution on fine T213 grid given by Jakob-Chien et al. (1995) is widely accepted as the reference solution to this test due to its good capability to reproduce the behavior of small scales.Numerical results of height fields by the GLPCC model are shown in Fig. 12 for tests on grids G 12 and G 24 at day 7 and 14.At day 7, no obvious difference is observed between the solutions on different grids and both agree well with the reference solution.At day 14, obvious differences are found on different grids.Eight circles of 8500 m exist in the results on the coarser grid G 12 , which are also found in the spectral transform solution on the T42 grid, but not in the results on the finer grid G 24 by the GLPCC model and the spectral transform solutions on the T63 and T213 grids.Additionally, the contour lines of 8100 m exist in spectral transform solution on the T213 grid, but not in present results and spectral transform solutions on the T42 and T63 grids.According to the analysis in Thuburn and Li (2000), this is due to the less inherent numerical viscosity on finer grids.As in case 5, total mass is conserved to the machine precision as shown in Fig. 13 and the conservation errors for total energy and potential enstrophy are given in Fig. 14 for tests with different resolutions.Total energy error of −6.131 × 10 −6 and potential enstrophy error of −1.032 × 10 −3 are obtained by the present model running on grid G 20 , which are smaller than those obtained by our fourth-order model on the 32 × 32 × 6 grid (Chen and Xiao, 2008).This test was also checked in Chen et al. (2014a) by a third-order model (see their Fig.19c  and d), where many more DOFs (9 times more than those on grid G 24 ) are adopted to obtain a result without eight circles of 8500 m at day 14.It reveals a well-accepted observation that a model of higher order converges faster to the reference solution, and should be more desirable in the atmospheric modeling.

Barotropic instability
A barotropic instability test was proposed in Galewsky et al. (2004).Two types of setups of this test are usually checked in the literature, i.e., the balanced setup and unbalanced setup.The balanced setup is same as Williamson's standard case 2, except the water depth changes with much larger gradient within a very narrow belt zone.This test is of special interest for global models on the cubed-sphere grid, since that narrow belt zone is located along the boundary edges between patch 5 and patches 1, 2, 3 and 4. Extra numerical errors near boundary edges would easily pollute the numerical results.
In practice, four-wave pattern errors may dominate the simulations on the coarse grids.For this case, we run the proposed model on a series of refined grids.By checking the convergence of the numerical results, we can figure out if the extra numerical errors generated by discontinuous coordinates can be suppressed by the proposed models with the increasing resolution.The unbalanced setup introduces a small perturbation to the height field.Thus, the balanced condition can not be preserved and the flow will evolve to a very complex pattern.Exact solution does not exist for unbalanced setup and a spectral transform solution on the T341 grid to this case given in Galewsky et al. (2004) at day 6 is adopted as reference solution.The details of setup of this test can be referred to Galewsky et al. (2004).

Balanced setup
We test the balanced setup at first.The proposed model runs on two grids with different resolutions of G 24 and G 72 .Numerical results of water depth after integrating for 5 days are shown in Fig. 15 and evolution of normalized l 1 errors of water depth of two simulations are depicted in Fig. 16.On a coarse grid with G 24 , the numerical result is dominated by four-wave pattern errors and the balanced condition can not be preserved in simulation.The accuracy is obviously improved by increasing the resolution using grid G 72 .The numerical result of height field at day 5 is visually identical to the initial condition.The improvement of the accuracy can be also proven by checking the velocity component u θ .Numerical results of u θ , which stay at zero in exact solution, vary within a range of ±31 m s −1 on grid G 24 and a much smaller range of ±0.8 m s −1 on grid G 72 .This test is more challenging for a cubed-sphere grid than other quasi-uniform spherical grids, e.g., a yin-yang grid or icosahedral grid.As shown in Fig. 16, at the very beginning of the simulation the l 1 errors increase to a magnitude of about 10 −4 on coarse grid G 24 and this character does not change on refined grid G 72 .This evolution pattern of l 1 errors are different from those of models on yin-yang and icosahedral grids, where initial startup errors also decrease on fine grids as shown in Chen et al. (2014a, Fig. 23).

Unbalanced setup
We run the unbalanced setup on a series of refined grids to check if the numerical result will converge to the reference solution on refined grids.Numerical results for relative vorticity field after integrating the proposed model for 6 days are shown in Fig. 17.Shown are the results on four grids with gradually refined resolutions of G 24 , G 48 , G 72 and G 96 .On grid G 24 , the structure of numerical result is very different from the reference solution.After refining the grid resolution, the result is improved on grid G 48 ; except for the structure in top-left corner, it looks very similar to the reference solution.On grid G 72 and G 96 , numerical results agree with the reference solution very well and there is no obvious difference between these two contour plots.Compared with the results of our former fourth-order model, the contour lines look slightly less smooth.Similar results are found in the spectral transform reference solution.Since this test contains more significant gradients in the solution, a highorder scheme might need some extra numerical dissipation to remove the noise around the large gradients.Increasing the grid solution can effectively reduce the magnitude of the oscillations as shown in the present simulation.

Conclusions
In this paper, a three-point high-order GLPCC scheme is proposed under the framework of flux reconstruction.Three local DOFs are defined within each element at Gauss-Legendre points and a super convergence of fifth order is achieved.This single-cell-based method shares advantages with the DG and SE methods, such as high-order accuracy, grid flexibility, global conservation and high scalability for parallel processing.Meanwhile, it is much simpler and easier to implement.With the application of the cubed-sphere grid, the global shallow-water model has been constructed using the GLPCC scheme.Benchmark tests are checked by using the present model, and promising results reveal that it is a potential framework to develop high-performance general circulation models for atmospheric and oceanic dynamics.
As any high-order numerical scheme, additional dissipation or limiter projection might be needed in simulations of real case applications.Because of the algorithmic similarity, the existing works on high-order limiting projection and artificial dissipation devised for DG or SE methods are applicable to the GLPCC without substantial difficulty.Also future studies should focus on designing more reliable limiting projection formulations for the GLPCC and other FR schemes, which are able to deal with discontinuities without losing the overall high-order accuracy.

Figure 1 .
Figure 1.Configuration of DOFs and constraint conditions in a onedimensional case.

Figure 2 .
Figure 2. The spectrum of the semi-discrete scheme.

Figure 3 .
Figure 3. Numerical dispersion (left) and dissipation (right) relations of the semi-discrete scheme.

Figure 6 .
Figure 6.Configuration of DOFs and constraint conditions in a twodimensional case.

Figure 7 .
Figure 7.The Riemann problem along patch boundary edge between patch 1 and 4.

Figure 8 .
Figure 8. Numerical results and absolute errors of water depth for case 2 on grid G 12 at day 5. Shown are water depth (top left) and absolute error (top right) of the flow with γ = 0 and water depth (bottom left) and absolute error (bottom right) of the flow with γ = π 4 .

Figure 9 .
Figure 9. Numerical results of total height field for case 5 on grid G 12 at day 5 (top left), day 10 (top right) and day 15 (bottom).

Figure 10 .
Figure 10.Normalized conservation error of total mass on grid G 12 for case 5.

Figure 11 .
Figure 11.Normalized conservation errors of total energy and potential enstrophy on refined grids for case 5.
's standard case 6: Rossby-Haurwitz wave The Rossby-Haurwitz wave case checks a flow field including the phenomena of a large range of scales.As a result, the www.geosci-model-dev.net/8

Figure 12 .
Figure12.Numerical results of water depth for case 6 on grid G 12 at day 7 (top left), day 14 (top right) and on grid G 24 at day 7 (bottom left) and day 14 (bottom right).

Figure 13 .
Figure 13.Normalized conservation error of total mass on grid G 12 for case 6.

Figure 14 .
Figure 14.Normalized conservation errors of total energy and potential enstrophy on refined grids for case 6.

Figure 15 .Figure 16 .
Figure 15.Numerical results of water depth for balanced setup of barotropic instability test on two grids G 24 (left) and G 72 (right).Contour lines vary from 9000 to 10 100 m.

Table 1 .
Numerical errors at two wave numbers and corresponding convergence rate.

Table 2 .
Numerical errors and convergence rates for advection of a sine wave.

Table 3 .
Numerical errors and convergence rates for case 2 of the flow with γ = π 4 .