Quadratic spline methods for the shallow water equations on the sphere: Galerkin
Introduction
Weather prediction is a science with a long history. Its objective is to describe and predict of the behaviour of the atmosphere, ocean water and sea ice. The accuracy of weather prediction depends on many factors, such as the accuracy of the knowledge of the state of the atmosphere at the initial time, the numerical methods applied, and the resolution used in these methods. Because weather prediction computations are time-consuming, there is interest among the scientific community in developing accurate and efficient methods for weather prediction. One way to achieve high accuracy in weather prediction computations is to consider high-order discretization methods.
The spatial discretization schemes that are commonly used in meteorological simulations are finite difference schemes, spectral schemes, and finite element schemes. For instance, the models used by the National Center for Atmospheric Research (NCAR) and the European Centre for Medium-Range Weather Forecasts (ECMWF) are based on the spectral transform method [9], [19], [20], while the model developed by the Canadian Meteorological Centre in partnership with the Meteorological Research Branch (CMC-MRB) uses a variable-resolution cell-integrated finite element scheme [4], [5]. Although the spectral transform method is popular and offers superior accuracy at a sufficiently high resolution, the efficiency of the method at high resolutions is reduced, because the cost of performing spectral transforms increases rapidly with spatial resolution. In the case of Fourier transforms in the longitudinal direction, fast Fourier transforms (FFTs) may be used and their computational cost grows as , where N is the number of spatial subintervals in one dimension. However, an efficient method for performing Legendre transforms, analogous to FFTs, has not yet been developed. Thus, the Legendre transforms in the latitudinal direction are often performed by summation and their costs escalate at a rate of . Therefore, finite element methods will likely offer a competitive alternative in the future, since, compared to the spectral transform method, their computational complexity grows at a slower rate. Provided that a sufficiently efficient solver is applied for the solution of the linear system arising from the Helmholtz problem, the computational costs of finite difference and finite element methods applied to the shallow water equations (SWEs) on the sphere increase quadratically with the number of gridpoints in one dimension (i.e., ). Moreover, finite element methods can incorporate adaptive grids and, compared to spectral methods, are more suitable for massively parallel computers.
In this paper and in a companion paper [12], we present finite element-based numerical methods for the SWEs in spherical coordinates. The SWEs, which describe the inviscid flow of a thin layer of fluid in two dimensions [8], have been used for many years by the atmospheric modeling community to test promising numerical methods for solving atmospheric and oceanic problems. Because the earth is approximately spherical, most global atmospheric models in use today are based on spherical coordinates. To define the equations on the sphere, let u and be the wind velocity components in the (longitudinal) and (latitudinal) directions, respectively, and be the geopotential. Let R be the radius of the earth, its rotational speed, and the Coriolis parameter, where R and are assumed to be constant. In spherical coordinates, the SWEs are given bywhere the Lagrangian derivative is defined byand the subscripts and denote the spatial derivatives in the two directions. Since u and are multi-valued at the poles, we adopt the approach of Côté and Staniforth [7] and compute the components of the wind images instead: and . Thus, we solve the following equations, obtained by multiplying the motion equations (1), (2) and the continuity equation (3) by and , respectively, by expressing the resulting equations in terms of U and V, and by isolating the nonlinearity in the continuity equation (3) in a logarithmic termwhere .
Along the longitude, functions are assumed periodic, whereas at the poles , homogeneous Dirichlet boundary conditions are imposed on the wind image components: . The latitudinal boundary conditions on are designed to mimic the behaviour of its spherical harmonic expansions [7], and are set to be homogeneous Neumann: .
In this paper, the SWEs are discretized in time using the two-level semi-Lagrangian semi-implicit (SLSI) method, and in space using the quadratic spline Galerkin (QSG) method. In the companion paper [12], the optimal quadratic spline collocation (OQSC) methods are used for the spatial discretization. Our numerical results show that the QSG method, when applied in conjunction with the SLSI method, is stable and non-dispersive for the meteorologically significant Rossby waves. Also, the method exhibits locally fourth-order spatial convergence, and compares favorably to the linear spline Galerkin method, which forms the basis of the method in [3].
Section snippets
Time discretization
Discretization schemes based on a semi-Lagrangian treatment of advection have generated considerable interest in the past decade for the efficient integration of atmospheric models, since they offer the promise of larger timestep size, with no loss in accuracy when compared to the Eulerian-based advection schemes, in which the timestep size is limited by more severe stability restrictions [15], [16].
A semi-Lagrangian time discretization scheme in spherical coordinates approximates the
Space discretization
To apply the QSG method to the time-discretized SWEs, we first define two uniform partitions along the longitude,where is an integer and denotes the meshsize in the -direction. The two partitions and are staggered with respect to each other. The gridpoints in and are chosen so that for , and for , respectively. Similarly, staggered partitions and are defined in the
Rossby wave stability
In this section, we perform a stability analysis for the QSG method applied to a simplified version of the SWEs and examine the conditions under which the discretized solutions are stable. We assume, for simplicity, an unstaggered grid and biperiodic boundary conditions. In [17] Staniforth and Mitchell showed that, if Cartesian coordinates are used, if the discretization grid is unstaggered, and if biperiodic boundary conditions are assumed, then the Helmholtz equation should be derived
Numerical results
We now demonstrate numerically the stability of our method with respect to the Rossby waves and study its convergence behaviour. In the companion paper [12], we compare the computational cost of the QSG method to those of the linear spline Galerkin method (on which the model in [3] is based) and three quadratic spline collocation methods; we also show that when the spatial discretization is done on the C-grid, gravity waves propagate in the proper directions.
Remarks
Spatial discretization schemes commonly used in meteorological applications are currently limited to spectral methods or low-order finite difference/finite element methods [4], [5], [9], [19], [20]. Spectral methods, which yield high-order solutions, give rise to dense matrices, therefore their parallel implementation is not scalable. In contrast, finite element methods, which give rise to sparse matrices (thus fewer global communications), have more potential for parallelism and give rise to
Acknowledgement
The authors wish to thank Dr. Steve Thomas for his helpful suggestions throughout this piece of research.
References (20)
- et al.
Quadratic spline methods for the shallow water equations on the sphere: collocation
Math. Comput. Simulat.
(2006) - et al.
Multiply-upstream, semi-Lagrangian advection schemes: analysis and application to a multi-level primitive equation model
Mon. Wea. Rev.
(1982) Equatorial solitary waves. Part I. Rossby solitons
J. Phys. Oceanogr.
(1980)Variable resolution techniques for weather prediction
Meteorol. Atmos. Phys.
(1997)- et al.
The operational CMC-MRB global environment multiscale (GEM) model. Part II. Results
Mon. Wea. Rev.
(1998) - et al.
The operational CMC-MRB global environmental multiscale (GEM) model. Part I. Design considerations and formulation
Mon. Wea. Rev.
(1998) - et al.
A two-time-level semi-Lagrangian semi-implicit scheme for spectral models
Mon. Wea. Rev.
(1988) - et al.
An accurate and efficient finite-element global model of the shallow water equations
Mon. Wea. Rev.
(1990) - et al.
Numerical Prediction and Dynamic Meteorology
(1980) Aspects of the numerics of the ECMWF model
Cited by (2)
Quadratic spline methods for the shallow water equations on the sphere: Collocation
2006, Mathematics and Computers in SimulationCitation Excerpt :In a companion paper [13], we present a finite element-based method, which discretizes the shallow water equations (SWEs) in time using the semi-Lagrangian semi-implicit (SLSI) method and in space using the quadratic spline Galerkin (QSG) method.
A Galerkin method with spherical splines for the shallow water equations on a sphere: error analysis
2015, Numerische Mathematik
- 1
Part of this work was done while the author was at the Department of Computer Science, University of Toronto and published as her doctoral thesis [18].
- 2
Supported in part by Natural Sciences and Engineering Research Council (NSERC) of Canada.