Quadratic spline methods for the shallow water equations on the sphere: Galerkin

https://doi.org/10.1016/j.matcom.2004.10.008Get rights and content

Abstract

Currently in most global meteorological applications, low-order finite difference or finite element methods, or the spectral transform method are used. The spectral transform method, which yields high-order approximations, requires Legendre transforms. The Legendre transforms have a computational complexity of O(N3), where N is the number of subintervals in one dimension, and thus render the spectral transform method unscalable. In this study, we present an alternative numerical method for solving the shallow water equations (SWEs) on a sphere in spherical coordinates. In this implementation, the SWEs are discretized in time using the two-level semi-Lagrangian semi-implicit method, and in space on staggered grids using the quadratic spline Galerkin method. We show that, when applied to a simplified version of the SWEs, the method yields a neutrally stable solution for the meteorologically significant Rossby waves. Moreover, we demonstrate that the Helmholtz equation arising from the discretization and solution of the SWEs should be derived algebraically rather than analytically, in order for the method to be stable with respect to the Rossby waves. These results are verified numerically using Boyd’s equatorial wave equations [J.P. Boyd, Equatorial solitary waves. Part I. Rossby solitons, J. Phys. Oceanogr. 10 (1980) 1699–1717] with initial conditions chosen to generate a soliton.

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 O(N2log(N)), 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 O(N3). 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., O(N2)). 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 v 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 f=2Ωsinθ the Coriolis parameter, where R and Ω are assumed to be constant. In spherical coordinates, the SWEs are given bydudtf+utanθRv+ϕλRcosθ=0,dvdt+f+utanθRu+ϕθR=0,dϕdt+ϕRcosθ[uλ+(vcosθ)θ]=0,where the Lagrangian derivative is defined byddtt+uRcosθλ+vRθ,and the subscripts λ and θ denote the spatial derivatives in the two directions. Since u and v 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: Uucosθ/R and Vvcosθ/R. Thus, we solve the following equations, obtained by multiplying the motion equations (1), (2) and the continuity equation (3) by cosθ/R and cosθ, 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 termdUdtfV+ϕλR2=0,dVdt+fU+cosθR2ϕθ+δ=0,cosθddtlogϕ+Uλcosθ+Vθ=0,where δ(U2+V2)sinθ/cos2θ.

Along the longitude, functions are assumed periodic, whereas at the poles (θ=±π/2), homogeneous Dirichlet boundary conditions are imposed on the wind image components: U(λ,±π/2)=V(λ,±π/2)=0. The latitudinal boundary conditions on ϕ are designed to mimic the behaviour of its spherical harmonic expansions [7], and are set to be homogeneous Neumann: ϕθ(λ,±π/2)=0.

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,Δλ{0=λ0<λ1<<λNλ=2π},ΔˆλΔλ2=λˆ1<λˆ0<<λˆNλ=2π+Δλ2,where Nλ>0 is an integer and Δλ=2π/Nλ 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 λi=iΔλ for i=0,,Nλ, and λˆi=(i+1/2)Δλ for i=1,,Nλ, 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)

  • A.T. Layton et al.

    Quadratic spline methods for the shallow water equations on the sphere: collocation

    Math. Comput. Simulat.

    (2006)
  • J.R. Bates et al.

    Multiply-upstream, semi-Lagrangian advection schemes: analysis and application to a multi-level primitive equation model

    Mon. Wea. Rev.

    (1982)
  • J.P. Boyd

    Equatorial solitary waves. Part I. Rossby solitons

    J. Phys. Oceanogr.

    (1980)
  • J. Côté

    Variable resolution techniques for weather prediction

    Meteorol. Atmos. Phys.

    (1997)
  • J. Côté et al.

    The operational CMC-MRB global environment multiscale (GEM) model. Part II. Results

    Mon. Wea. Rev.

    (1998)
  • J. Côté et al.

    The operational CMC-MRB global environmental multiscale (GEM) model. Part I. Design considerations and formulation

    Mon. Wea. Rev.

    (1998)
  • J. Côté et al.

    A two-time-level semi-Lagrangian semi-implicit scheme for spectral models

    Mon. Wea. Rev.

    (1988)
  • J. Côté et al.

    An accurate and efficient finite-element global model of the shallow water equations

    Mon. Wea. Rev.

    (1990)
  • G.J. Haltiner et al.

    Numerical Prediction and Dynamic Meteorology

    (1980)
  • M. Hortal

    Aspects of the numerics of the ECMWF model

There are more references available in the full text version of this article.

Cited by (2)

  • Quadratic spline methods for the shallow water equations on the sphere: Collocation

    2006, Mathematics and Computers in Simulation
    Citation 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.

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.

View full text