In this continuation paper to Starius (2018), referred to as (S-E-P) below, we consider numerical methods for the shallow water equations on the sphere, using the Equator-Pole (E-P) grid system. We make use of a reduced grid system for an annular band around the equator, which is divided into segments, each covered with a lat-lon net. As in (S-E-P) the equatorial grid is complemented by two orthogonal polar grids obtained by suitably modified stereographic coordinates. The whole grid system is optimised by minimising a global measure of uniformity, namely the maximum of all deviation factors on the sphere, where as in (S-E-P), the deviation factor for a net rectangle is the greatest ratio of its side-lengths.
The methods of this paper give an accuracy in the total mass, which is considerably higher than needed for NWP. This was already shown for the smooth Rossby-Haurwitz wave example for the method in (S-E-P). Here it is also verified, with and without reductions, for non-smooth solutions by considering the Cosine bell no. 1 and the mountain problem no. 5 in Williamson et al. (1992), the former with second order and the latter with first order discontinuous derivatives in their solutions. We have also integrated the mountain problem for a time period of 5 years, with maximal relative errors in the total mass of order For the total energy the corresponding quantities were of order
If K denotes the number of segments on a hemisphere, then there are 2 K-1 segments on the whole sphere, since the two segments closest to the equator are merged together. We note that K = 1 is the case considered in (S-E-P). By reductions the uniformity is improved and the total number of grid points can be decreased by 20% or more, with no substantial change in accuracy. For high resolution also the computing time will be reduced by a similar percentage. High order one-dimensional centred interpolation formulas are used to connect the discretisations on adjacent segments. The kind of deviations from uniformity that arise by changing the number of grid points on parallels do not cause any problems, at least not for the methods described in this paper.
We now briefly comment on some other methods considered in the literature. Our main sources are the two overview papers Staniforth and Thuburn (2012) and Williamson (2007). There are many other grid techniques than ours for the sphere, among which the most popular are the cubed sphere, the icosahedral and the Yin-Yang grid systems. The paper (S-E-P) contains a more detailed overview than the one given below.
Cubed sphere grids (Sadourny, 1972) consist of quadrilaterals and have been used for finite volume methods in, e.g. (Putman and Lin, 2007; Ullrich et al., 2010), for discontinuous Galerkin methods in, e.g. (Nair et al., 2005; Bao et al., 2014) and for finite difference methods in, e.g. (Ronchi et al., 1996). Some drawbacks are lack of orthogonality and deviation from uniformity, particularly at the 8 corner-points of the cube, where angles up to instead of the ideal occur.
Icosahedral grids (Williamson, 1968) are generally triangular and have been used for finite volume methods in, e.g. (Pudykiewics, 2011; Chen et al., 2014) and for discontinuous Galerkin methods in, e.g. (Giraldo, 2006; Läuter et al., 2008). At the 12 corner-points of the icosahedron there are angles up to instead of the ideal
The more recently developed Yin-Yang grid, in Kageyama and Sato (2004), is an overset grid system consisting of two congruent orthogonal subgrids based on spherical coordinates and having maximal deviation factor This grid has been used by, e.g. Qaddouri and Lee (2011) and for a full 3d model by Allen and Zerroukat (2016). Both this grid system and ours use overset and orthogonal grids, which leads to an obvious similarity between the two techniques.
The great importance of the axis of the Earth for the dynamics of the weather is clear. However, the Yin-Yang grid does not have good rotational symmetry properties relative to this axis, not even near the poles, for details see, e.g. Section 1 in (S-E-P). At the end of Section 3.4, we compare the Yin-Yang grids and ours with respect to the total number of grid points, both with the same number on the equator.
The paper is organised as follows. In Section 2, we introduce the set of governing equations, some notation, and also some coordinate systems. This section is essentially a summary of Section 2 in (S-E-P).
Section 3 is devoted to reduced grid systems for annular bands around the equator. Two methods, I and II, to divide the band into segments and then to connect them by interpolation, are considered in Sections 3.1 and 3.2. Method I is quite general and easy to implement but requires more interpolation coefficients than method II. However, method II has limitations in its applicability, which is discussed in Section 3.2. In Section 3.3, two different measures of uniformity, used to generate suitable grid systems for the sphere, are compared. In Section 3.4, an E-P grid system is constructed with square net rectangles at mid-latitudes, which is often required for NWP models.
In Section 4, we investigate our E-P grid systems numerically, both with and without reductions. Centred finite differences of high order are used for the spatial discretisations and the explicit fourth order Runge-Kutta method for the time integration. First some notation for various errors is introduced and then, in Section 4.1, smoothing by addition of a ’hyper-diffusion’ term is considered. For a discretisation method of order 2p a smoothing term of the same order is used, in contrast to the method in (S-E-P) where the order is always equal to 4. In Section 4.2, the rotation of the Cosine bell is considered and for this scalar advection example both the total mass and energy are obtained with much higher accuracy than for the height of the fluid. In Section 4.3, our results for an example by Läuter et al. (2005) show that the reduction techniques work quite well, the accuracy is very much the same for K = 1,2,3,4, particularly for the 6th order method. In Section 4.4, we consider zonal flow over an isolated mountain, example no. 5 in Williamson et al. (1992), with non-smooth solutions. The errors in the total mass shown in Table 6 are clearly small enough for NWP. We mention that no ‘corrections’ or ’fixers’ of any kind have been used.
In (S-E-P) further numerical examples are given, and also investigations of grid imprinting and comparisons between formal and computational orders of discretisation errors are considered.
Coordinate systems and governing differential equations
In this section, we introduce some notation and give a frame of reference for the paper. It is essentially a summary of Section 2 in (S-E-P).
Spherical coordinates and the shallow water equations
The transformation between Cartesian and spherical coordinates on the sphere is
Standard and modified stereographic coordinates
Consider first the standard stereographic coordinates. The relation between these and spherical coordinates can be expressed as
Since the velocities in spherical coordinates, are not uniquely defined at the poles, we use instead the stereographic velocities Observe the use of the minus sign in the definition of U, intended to get nicer formulas, namely
Stereographic projection is only suitable for grid generation close to the poles. Here we will use a modification of (3), already recommended in Section 2.3 of (S-E-P), so that the lat-lon and our new modified stereographic grid points will coincide on the meridian corresponding to λ = 0, and automatically also on the meridians provided the number of grid points on parallels in the last segment of the reduced grid is divisible by 4. The assumption for λ = 0 leads to new coordinates related to by
The modified coordinates give rise to more uniform polar grids and will also simplify the connection between these grids and the last segments of the equatorial grid system. The shallow water equations expressed in the coordinates defined by (3) and (5), and the stereographic velocities U and V given in (4), are written down in full detail in Section 2.3 of (S-E-P).
Reduced latitude-longitude grid systems covering an equatorial band
In this Section, we consider the possibility to utilise a reduced grid system for an annular band around the equator, cf. Gates and Riegel (1963), Kurihara (1965), Tolstykh and Shashkin (2012), and Starius (2014). We also mention the reduced Gaussian grids for spectral models in Hortal and Simmons (1991), reported to decrease the computing time by about 20%, and used by ECMWF in their operational model. Further, the papers Layton (2002) and Li (2018) are worth mentioning in this context. The band is divided into segments, each covered with a lat-lon net. Let n be the given number of grid points on the equator and K the number of segments on a hemisphere. Here we choose to have invariance under reflection in the equator, where the closest two segments are merged together, which implies that the total number of segments is 2 K-1 and the total number of subgrids is 2 K + 1. We recall that K = 1 is the case studied in (S-E-P).
We will later define latitudes such that the i:th segment Si is bounded by and θi, where Also the number of grid points ni on a parallel in Si has to be defined and satisfy and the longitudinal steps are defined by where We note that and when only one longitudinal step is used, for example for K = 1, it is called
We now discretise the segments introduced above. Assume that all the latitudes θi are multiples of a latitudinal step Consider for each i, with the parallels where and let be the corresponding lat-lon net. Note that we have let and both contain the parallel corresponding to θi, which will ease the presentation in Sections 3.1 and 3.2. When the discrete segments have been defined, we omit the parallel θi from i = 1,…K-1, cf. Fig. 1.
In this section, we will often use the fact that for a parallel with latitude θ and angular longitudinal step the distance along the parallel between two consecutive grid points is
Method I: reduced latitude-longitude grid system
Here the number of segments K will be given by the user. The deviation factor at the equator is where will be the same for all discretised segments and determined later. We consider the segment and require that the distances between grid points on the parallel are equal to that is where Further, the maximal deviation factor for the segment is which is obviously minimised, with respect to when the two deviation factors in the above are equal. We now have the equalities
Now can be evaluated and for practical reasons modified by defining the integer m to be rounded to the nearest integer and then redefining With known we now round the latitudes θi, to their closest multiples of Finally we set round it to the nearest integer, and define
In Table 1, the overlap between the last segment in the reduced grid and the polar grid is 3, which only affects the number of parallels in this segment. An increase in the number of segments leads to a decrease in the total number of discretisation points and in the maximal deviation factor. In Fig. 1, a reduced grid with n = 160 and K = 3 is illustrated. More realistic resolutions for NWP are for instance n = 1000, K = 10, which gives p = 20.7% and and n = 2000, K = 10, which gives p = 20.8% and We recall that the above procedure requires only one-dimensional interpolation and that the coefficients in the formulas should be computed and stored during preprocessing. The number of sequences of interpolation coefficients is which for high resolution is negligible compared to other storage requirements. A less general method using considerably fewer interpolation formulas is briefly described in Section 3.2.
A way to avoid scalability problems for the couplings between the reduced grid and the northern polar grid, say, is to let the former reach a high latitude. Since θK increases slowly with K (for instance K = 10 gives only ) we will now consider a modification of method I, which is suitable for high latitudes. We choose K and determine the first K segments and the corresponding maximal deviation factor σK, and then complement this by the desired number of additional segments, each optimised with respect to σK. More specifically, additional latitudes can be determined by the equation and then also where the ni have been rounded in a suitable way.
We conclude this subsection with two examples. For initial K = 3, we have and the corresponding latitudes are: 34.0, 46.8, 55.1, 61.9268, 67.0595, 71.1657, 74.4912, 77.2050, 79.4305, 81.2615. If we want the reduced grid to reach the polar circles, then we add two segments, so that K = 5, and reach latitude The maximal deviation factor for the whole reduced grid system is still σ3. In a second example, with initial K = 1 and two segments have to be added, so that K = 3, for the reduced grid to pass the polar circles. There is obviously a trade-off between uniformity and the possibility to use a small number of segments.
Method II: reduced latitude-longitude grid system
This method is essentially applicable with a small number of segments and special choices of n. However, its range of applicability can be extended in ways described below. The notation and assumptions of Section 3.1 will be used here, and only a brief outline of the method will be given.
The method is based on the following simple idea. In a consecutive manner we replace q longitudinal mesh intervals in segment with q−1 intervals in the following segment As for method I we require that where and further that each ni, with i < K, is divisible by q. This implies that Further, together with (7) leads to where we write σ instead of σK, since σ is now independent of K. We choose q and determine σ and where is the width of the reduced grid, by using
We now determine θi by until Further, we can choose K from and then redefine For low resolution it might be necessary to replace K by a smaller value in order to get a reasonable number of parallels in the segment Since is known we can use -1, and then define
The above works only if n is divisible by where K has our desired value, in which case the number of sequences of interpolation coefficients required is only q−1, provided symmetry is used to reduce the number of sequences whenever possible.
If we get a number with not divisible by q, then will not correspond to grid points in Since we now cannot have equidistant grid points in all our discretisations, we switch to a modification of method II, which allows that. From segment onward, we define for and then round them to integers and define for This transition to method I is very simple but requires more sequences of interpolation coefficients.
For high resolution, we indicate an improvement of the above modification of method II. If is rounded to an integer divisible by q, then the interpolation between and can be done with the q−1 interpolation formulas. We still use method I for the interpolation between and This can be repeated by choosing divisible by q, using method I between and and then the q−1 interpolation formulas between and etc.
We conclude this subsection by considering an example with q = 6, which preliminarily gives We have and the following sequence of latitudes determining the segments: 33.6, 46.0, 54.6, 61.2, 66.3, 70.4, 73.8, 76.6, 78.8, 80.7. Thus according to the above we choose K = 3, redefine and obtain a new If n is divisible by 62, then method II works perfectly well for K = 3. However, if we want K > 3 and n is not divisible by 63, we can use a modification as described above for the segments following
Comparison between two measures of uniformity
The measure of uniformity most commonly referred to in the literature is the ratio of maximum to minimum grid length, which we call γ in this section. Our main interest here is to investigate which measure of uniformity is best suited for grid generation. Our maximal deviation factor has already been used for this purpose, and we therefore now turn our focus to the ratio γ by studying two examples.
Example (i): Let Sd be a segment on the northern hemisphere bounded by the latitudes and consisting of a lat-lon net with the longitudinal step We define σ as the deviation factor for net rectangles at θ1, that is where is the latitudinal step, not yet determined. The ratio of maximum to minimum grid length on Sd is
The minimum of γ is which is attained for all The result is unsatisfactory since we want the grid to be uniquely defined.
For comparison we mention that the minimum of the maximal deviation factor for this segment is which is attained if and only if We get square net rectangles at the latitude θ for which implies
This formula will be used in the next Section 3.4. ⊗
Example (ii): In this example the ratio γ will be minimised over the whole sphere. We choose the E-P grid system with K = 1, which we optimised in (S-E-P) by using the maximal deviation factor. Now instead we minimise γ. Let the equatorial grid cover the band which will be determined later. Our northern polar grid will go down to and the corresponding map factor is According to (12) in (S-E-P) the greatest and smallest grid lengths in our polar grids are and respectively, where For our polar grids the ratio γ and the maximal deviation factor are the same.
The ratio γ for the whole sphere can now be written as
where We now consider the equation with a unique root in the interval Elementary computations show that
We look at two cases, namely
In both cases with equality if and only if σ = 1 and This result was also reported in (S-E-P) but without mentioning the condition σ = 1, which we regret.
The sole reason why lat-lon coordinates are not suitable for the whole sphere is the convergence of meridians when we are moving towards a pole. Therefore, we find it inappropriate to use square net rectangles at the equator instead of rectangles with longer longitudinal than latitudinal sides. ⊗
The examples (i) and (ii) above convincingly show that the maximal deviation factor is a better tool for generation of suitable grid systems on the sphere, than the ratio of maximum to minimum grid length.
Square net rectangles at mid-latitudes
Grids for NWP are generally constructed to be uniform (having square net rectangles) at mid-latitudes, a condition we will study in this section. We will also compare the total number of grid points used in the Yin-Yang grids and ours, both with the same number of grid points on the equator, and by considering the mid-latitude condition above.
For our ’mid-latitude method’ we use notation with bars to distinguish the variables from those used for methods I and II. In order to avoid over-determined systems of equations we replace equation (8) by (12) and use the left part of (7) to get which together with (12) gives We then determine and by solving the equations
We conclude this subsection with an example, in Table 2, where two reduced lat-lon grids are considered. Choose n = 3600 and K = 10 in method I, which corresponds to on the equator and about 11 km between neibouring grid points everywhere on the sphere.
For the Yin-Yang grid we let the Yin subgrid have square net rectangles at If the Yang grid is congruent with the Yin grid the total number of grid points will be which for n = 3600 is about This grid system only satisfies the mid-latitude condition in a crude sense. By replacing this grid by one of ours in Table 2 we get a reduction in the number of grid point of about 34%. Moreover, we get good rotational symmetry properties relative to the axis of the Earth, much better uniformity, and perfect possibilities to satisfy the mid-latitude condition in suitable ways.
Provided we accept less uniformity than in Table 2, a smaller number of segments can be used, e.g. K = 3 or 5. If we also want the reduced grid to reach high latitudes it can be complemented by additional segments by the procedure described at the end of Section 3.1.
For a full lat-lon grid, satisfying the mid-latitude condition and with n = 3600, the total number of grid points is about
In this section, we investigate the E-P grid system numerically, both with and without reductions in the equatorial grid. Test examples no. 1 and 5 from Williamson et al. (1992) and time dependent flow from Läuter et al. (2005) are considered, cf. also (S-E-P) where other numerical examples and investigations, such as grid imprinting, are studied. We recall that the number of segments in the equatorial grid is 2 K-1, that the total number of grids in the E-P system is 2 K + 1, and that K = 1 corresponds to the grid system used in (S-E-P). Method I, intended for reduction and connection between segments, has been used in all the tables below. However, in Sections 4.3 and 4.4 we give some short comparisons between methods I and II.
The spatial discretisation is obtained by using centred finite difference approximation of orders 6 or 4 on non-staggered grids, for which formulas can be found in Section 2.4 of (S-E-P). The classical explicit fourth order Runge-Kutta method is used for the time integration, which requires interpolation between grids for each of the four stages in a time step.
Explicit time integration seems to work quite well for the shallow water equations. However, for realistic NWP, 3d models are used and with vertical spacing much smaller than the horizontal ones, which would lead to very small time steps for explicit methods. This is why HEVI(horizontal explicit and vertical implicit) schemes are of interest in the atmospheric community.
In the tables below, 2p is the approximation order of the spatial discretisations of the underlying differential equations. The order of the interpolation formulas connecting the polar and equatorial grids is 2p + 1, and between discretised segments the order used is 2p + 2. The overlapping number is denoted by cf. Section 3.3 in (S-E-P).
The geopotential is defined as where g is the acceleration of gravity and H the height of the fluid. We will often use the relative error norm introduced in Williamson et al. (1992). Let denote this norm after d days. For the vector where T is the integration time in days, we define a mean value by
The global conservation properties for mass and energy will also be studied below. The integral formulas can be found, e.g. in Starius (2014) and their numerical computation is described in (S-E-P). The relative errors for the total mass after each day are stored in the vector with mean value defined by (14). The quantity for the total energy is defined analogously.
Smoothing is essential for the centred finite difference methods to work properly and will be achieved by adding a so called ‘hyper-diffusion’ term to each equation. Since ‘hyper-diffusion’ has no known physical meaning, we consider our smoothing a purely numerical stabilisation process.
Let Lh be the five-point operator associated with minus the Laplacian and defined by
To the first, second, and third equation in (2), we add the terms and respectively. In the equations for the polar grids, U and V are used instead of u and v. Since our grid systems are almost uniform, it is natural to use Lh for the whole sphere. We observe that for any smooth function In the temporal discretisations the smoothing terms will be multiplied by The smoothing constant will often be determined by minimising or We have used the MATLAB® routine fminbnd, for the numerical minimisation.
In an earlier work we have used Lh as a five-point difference approximation to minus the Laplacian in spherical coordinates, without any improvement for almost uniform grids.
We point out that the optimisations in this paper, are essential in order to ease comparisons between different methods.
Rotation of the Cosine bell
In this example we consider a solid body rotation of the Cosine bell, according to test problem no. 1 from Williamson et al. (1992). The advecting wind is given by
The E-P grid system will be used for K equal to 1 and 3 and the smoothing parameter c will be determined by minimising and The results are given in Table 3, first for and then for
Despite the solution having discontinuous second order derivatives, the accuracy of H is considerably higher for the case than for The actual approximation orders can still be fairly equal. The accuracy for the total mass must be caused by cancellation of errors since it is much higher than for H. In this specific example, the same can be said about the total energy. Minimisation of leads to slightly more accurate values of H, for K = 1 than for K = 3, but the values of are much smaller for the smoother case K = 3, which also require fewer grid points.
For the rest of this subsection we only consider the most important case The minimisation of leads to much greater values of c than for the other alternatives, which results in greater values for and see for example the case K = 1. For this reason we reject the use of for the determination of c and note that works better for this example if we focus on the other errors and and not on the values of c themselves, which are much smaller compared to the ones corresponding to
Time dependent flow, the Läuter example
We make some comments concerning Table 4. For order the errors decrease somewhat from the case to 200. This indicates that there is no real cancellation between the spatial and temporal local discretisation errors. For order there seems to be cancellation, and more so for than for For the accuracy obtained for K = 1,2,3,4 is very much the same, particularly for For each case the error grows very close to linearly in the number of days.
The third row of Table 4 will be if method II is used instead of method I. For this very smooth example there is no real difference in accuracy between the two methods.
In Table 5, is minimised in order to investigate whether it may be used instead of since the latter quantity is not easily available for realistic NWP. The values of c in the second to fourth rows in Table 5 are much smaller than those in Table 4, but the corresponding errors do not differ very much, except for the two cases K = 3, 4 and T = 15 days. However, the minimisation of might still give some insight in certain situations. In the next section, we will look at another way to determine values of c.
The solutions (15) are very smooth and have very small variation, which explains the high accuracy for the geopotential, and for the total mass and energy cancellation of errors has further increased the accuracy, see Table 5.
Zonal flow over an isolated mountain
Test problem no. 5 in Williamson et al. (1992) is concerned with zonal flow impinging on a mountain on the sphere, with a circular base. The centre of the base is denoted by and the initial conditions are
It is fairly easy to experimentally determine a ’stability interval’
In Table 6, the conservation properties are about the same for and which is natural since they are caused by cancellation of errors in the height of the fluid rather than by accuracy. We recall that q is the number of longitudinal mesh intervals that is replaced by q−1 intervals in the following segment, cf. Section 3.2. Only q−1 interpolation formulas are required and repeated periodically. This can be seen as a more regular way to interpolate between segments than the one of method I, which has resulted in greater accuracy for the total mass. The values of c in Table 6 have been obtained by using some tests.
In order to further investigate the reliability of our Equator-Pole grid systems, we have integrated the mountain problem for 5 years, using methods of order 6. In Table 7, the maximal relative errors in magnitude of the total mass appear in the fourth column and the corresponding quantity for the total energy in the last one, both after 5 years. The values of c are taken from Table 6.
The Equator-Pole grid system introduced in (S-E-P), for the solution of the shallow water equations on the sphere, consists of one equatorial and two polar grids. Here we generalise this by permitting the equatorial grid to be reduced, that is divided into segments. We give two methods of reduction, called method I and method II, in Sections 3.1 and 3.2, respectively. The reduction technique seems to work quite well, as can easiest be seen in Table 4.
By using reductions the total number of grid points on the sphere can be decreased by 20% or more, and for high resolution also the computing time by a similar percentage. The change in accuracy for the height of the fluid is negligible, and for the total mass the accuracy is increasing with increased uniformity of the grid system, cf. Table 3. High uniformity might have additional advantages, for example in conjunction with comprehensive wave propagation. The different wave speeds computed will, among other things, depend on the uniformity of the grid system used and its resolution. Low uniformity or resolution might lead to interaction between waves with no analogue in the continuous model.
We have demonstrated high accuracy in the total mass for fairly simple problems, both with and without smooth solutions. For realistic NWP models the situation can be more complicated, e.g. water can contribute to the atmosphere’s mass and give rise to far less smooth problems than those considered here. We mention that Allen and Zerroukat (2015) have derived a two-dimensional generalisation of the shallow water equations, which includes moist.
We use high order centred finite difference approximations of first order spatial derivatives, in the shallow water equations. This requires only a fraction of the work per grid point compared with some other methods of the same order. This fraction decreases rapidly with the order of approximation. Further, one can very easily increase the order of approximation by increasing the order of approximation for first order derivatives. For methods of order 4, 6, 8 and 10 the corresponding spatial stencils contain only 9, 13, 17 and 21 grid points, respectively. The higher the order of the spatial approximation, up to a bound depending on the resolution, the better treatment of fronts in the solutions, see for instance Table 2 in (S-E-P). For the time integration, we have used the explicit fourth order Runge-Kutta method.
The Equator-Pole grids presented here only focus on uniformity, except in Section 3.4, and are mainly intended to be basic grid systems. To make them more flexible, for instance when taking jet streams and orography into account, they need to be modified, perhaps by using overlapping techniques.