A discontinuous/continuous low order finite element shallow water model on the sphere

https://doi.org/10.1016/j.jcp.2011.11.018Get rights and content

Abstract

We study the applicability of a new finite element in atmosphere and ocean modeling. The finite element under investigation combines a second order continuous representation for the scalar field with a first order discontinuous representation for the velocity field and is therefore different from continuous and discontinuous Galerkin finite element approaches. The specific choice of low order approximation spaces is attractive because it satisfies the Ladyzhenskaya–Babuska–Brezzi condition and is, at the same time, able to represent the crucially important geostrophic balance.

The finite element is used to solve the viscous and inviscid shallow water equations on a rotating sphere. We introduce the spherical geometry via a stereographic projection. The projection leads to a manageable number of additional terms, the associated scaling factors can be exactly represented by second order polynomials.

We perform numerical experiments considering steady and unsteady zonal flow, flow over topography, and an unstable zonal jet stream. For ocean applications, the wind driven Stommel gyre is simulated. The experiments are performed on icosahedral geodesic grids and analyzed with respect to convergence rates, conservation properties, and energy and enstrophy spectra. The results match quite well with results published in the literature and encourage further investigation of this type of element for three-dimensional atmosphere/ocean modeling.

Introduction

Finite element (FE) schemes are successfully employed in numerous computational fluid dynamics applications. FE methods rely on a powerful mathematical apparatus and offer a good representation of the physical fields that are approximated using sets of basis functions of chosen accuracy. FE methods are applicable to unstructured grids and can accommodate grid refinement.

Today most dynamical cores of global atmosphere and ocean models are based on finite difference, or spectral transform schemes – mostly due to the ease of implementation and computational efficiency [4], [26]. Compared to those discretization methods, FE schemes allow more flexibility with regard to unstructured or locally refined grids and can approximate physical fields with higher order polynomials. There are dynamical cores of global ocean or atmosphere models based on FE schemes using continuous low order finite elements methods [10], spectral elements [33], or discontinuous Galerkin (DG) FE methods [11], [13], [23], [19].

Global scale presents special challenges for the dynamical core of an atmospheric or ocean model. In a finite element framework, four important properties need to be considered:

  • 1.

    The core needs to enable high resolution simulations.

  • 2.

    To support a smooth representation of topography and coastlines, one should be able to use very small grid cells. The second property is different to the first one since the internal resolution inside the grid cells can vary for different discretization schemes.

  • 3.

    The core must be able to represent the geostrophic balance between the Coriolis force and the pressure gradient force.

  • 4.

    The element must satisfy the Ladyzhenskaya–Babuska–Brezzi-condition (LBB-condition) – a necessary condition for stability of the discretization scheme.

The first two properties suggest low order elements. Unfortunately, it is difficult to find a low order element that fulfills the properties three and four at the same time. Failing in one of them leads to spurious modes [21]. Therefore, continuous low order finite element methods often require stabilization schemes [10].

A hybrid FE approach that combines a second order continuous representation for the scalar field with a first order discontinuous representation for the velocity field P1DGP2 was proposed as a potential candidate to combine the four properties above. It was shown for the linear equations that the element is LBB-stable and is able to represent the geostrophic balance [6], [7]. The analysis of the P1DGP2 element in terms of a Helmholtz decomposition for the linearized shallow-water equations showed that the element has no spurious pressure or Rossby modes. Spurious modes occur in the least harmful place via inertia oscillations that do not propagate [5].

In this paper we contribute to the understanding of the P1DGP2 element. We use the element to solve the spherical shallow water equations with and without viscous dissipation and on unbounded domains as well as on domains with lateral boundaries. The solutions are analyzed with respect to their convergence, conservation properties, and energy spectra.

While there exists extensive literature on how to obtain stable finite element models when only continuous or only discontinuous field representations are used, it is not obvious how to obtain an adequate model setup for the hybrid continuous/discontinuous P1DGP2-element. As opposed to continuous FE configurations, we need to solve a Riemann problem at cell boundaries due to the discontinuous velocity field, and, for a stable model setup, we need to perform a partial integration of the nonlinear divergence term in the scalar equation. As opposed to typical DG approaches, we use the non-conservative form of the shallow water equations and non-orthogonal Lagrange polynomials as basis functions. Since the continuous scalar field is well defined on cell boundaries, our Riemann problem reduces to two dimensions for the two components of the velocity field, this is different to DG methods. Except for the Riemann solver, our FE configuration is similar to the one presented in [9] for a P1DGP2 FE model on the plane. In [9] the weak form of the equations was derived from the non-conservative shallow water equations while the three-dimensional Riemann solver for velocity and scalar fields was deduced from the conservative form.

The introduction of the spherical geometry for global FE methods requires some care. To develop the physical fields into sets of basis functions, the triangles in physical space need to be mapped onto a reference triangle on which the basis functions are defined. On the sphere, the geometry of the physical triangles is given by trigonometric functions and cannot be mapped exactly onto a reference triangle, which is typically defined on the plane. In the literature, three different approaches can be found to introduce curved manifolds, such as the sphere, to finite element models. In the first approach, the differential equations are written for the curved manifold. This is done in the case of the cubed sphere, for example [28]. In the second approach, the differential equations are formulated in the three-dimensional space and a global Cartesian coordinate system is considered. The vector fields are forced to stay on the manifold via constraints [12]. In the third approach, the vector fields are written in the local tangent basis while the fluxes and spatial operators are expressed in the three dimensional Cartesian basis [3].

We propose the use of the stereographic projection to introduce the spherical geometry to global finite element models. The projection is part of the first approach mentioned above. It has already been used for finite difference and finite volume methods ([25], [8], [17], CCSR Ocean Component Model). In the stereographic projection, the sphere is projected from one of the poles onto a plane at the opposite side of the sphere. The pole itself is mapped to infinity. Global circulation models either use a combination of two stereographic projections from each pole that are connected at the equator or two stereographic projection from each pole that are connected to a Mercator projection in tropical regions [2], [25].

For an ocean model, it is sufficient to use one projection since the pole can be placed on land. For a global atmosphere model, we use two projections from each pole that are connected at the equator. The scaling factors that appear can be represented exactly by second order polynomials. The projection leads to a manageable number of additional terms in the differential equations.

In section two, we give an overview of the model setup, including the equations of motion, the discontinuous/continuous finite element discretization, and the incorporation of the spherical geometry. In section three, we apply the model to the standard test set for global shallow water models and a Stommel gyre ocean test case. In section four, we present a summary and conclusions.

Section snippets

The viscous shallow water equations

The viscous shallow water equations can be written in several different forms which are analytically equivalent but have different properties when discretized.

One of those forms is the non-conservative formtu+u·(u)+fk×u+gh-1H·(Hν(u))=τ,th+·(Hu)=0,where u is the two dimensional velocity vector, f is the Coriolis parameter, k is the vertical unit vector, g is the gravitational acceleration, ν is the eddy viscosity, τ is a forcing term (for example, bottom friction or wind stress in ocean

Numerical results

We evaluate our model using the standard suite of global shallow water test cases for the inviscid equation (ν = 0) as well as a non-linear Stommel gyre test with boundaries for the viscous case.

Geodesic grids

For our numerical experiments, we use icosahedral geodesic grids [1], due to their quasiuniform coverage of the sphere. In principle, each spherical triangular grid could be used for simulations. The grids can be refined to the next level by bisecting the sides of the triangles; the new

Conclusions

We investigate the potential applicability of the P1DGP2 finite element pair in atmosphere and ocean modeling. We present a stable P1DGP2 discretization for the viscous and inviscid shallow water equations on the sphere. We do not observe any spurious modes that necessitate additional diffusion or stabilization schemes.

Regarding the convergence properties and kinetic energy spectra, our discretization shows the expected behavior. A comparison to the ICON shallow water model shows very promising

Acknowledgments

We thank Hui Wan and the two anonymous reviewers for their useful and constructive comments, which helped to improve earlier versions of the paper.

References (35)

  • Matthias Läuter et al.

    A discontinuous Galerkin method for the shallow water equations in spherical triangular coordinates

    Journal of Computational Physics

    (2008)
  • Matthias Läuter et al.

    Unsteady analytical solutions of the spherical shallow water equations

    Journal of Computational Physics

    (2005)
  • C. Ronchi et al.

    The ‘cubed sphere’: a new method for the solution of partial differential equations in spherical geometry

    Journal of Computational Physics

    (1996)
  • Khosro Shahbazi

    An explicit expression for the penalty parameter of the interior penalty method

    Journal of Computational Physics

    (2005)
  • Mark Taylor et al.

    The spectral element method for the shallow water equations on the sphere

    Journal of Computational Physics

    (1997)
  • Paul A. Ullrich et al.

    High-order finite-volume methods for the shallow-water equations on the sphere

    Journal of Computational and Physics

    (2010)
  • David L. Williamson et al.

    A standard test set for numerical approximations to the shallow water equations in spherical geometry

    Journal of Computational and Physics

    (1992)
  • Cited by (12)

    • Numerical Methods for the Nonlinear Shallow Water Equations

      2017, Handbook of Numerical Analysis
      Citation Excerpt :

      To remove its disadvantage of pole-singularity, other choices of meshes have been proposed and studied, including the Yin–Yang mesh, cubed-sphere mesh, among others. Well-balanced and positivity-preserving methods under these meshes, including spectral element, finite volume, continuous low order finite elements and DG methods, have also been investigated, and we refer to Taylor et al. (1997), Giraldo et al. (2002), Nair et al. (2005), Verkley (2009), Ullrich et al. (2010) and Duben et al. (2012) for recent development on numerical methods for the SWEs on the sphere. Solving the two-layer SWEs is a challenging problem due to several reasons: they are only conditionally hyperbolic; they contain nonconservative product terms; they admit steady state to be exactly preserved; their water heights should stay nonnegative, and their eigenstructure cannot be obtained in explicit form.

    • A novel Local Time Stepping algorithm for shallow water flow simulation in the discontinuous Galerkin framework

      2016, Applied Mathematical Modelling
      Citation Excerpt :

      Shallow water flow equations, also known as Saint–Venant equations, are widely used in hydraulic and river engineering problems. Many researchers, with the use of different numerical methods such as finite volume, finite difference, finite element, hybrid finite element/volume, smoothed-particle hydrodynamics, and lattice Boltzmann, have successfully solved these equations [9–15]. In recent years, the discontinuous Galerkin (DG) finite element method has gained popularity in numerical modeling of shallow water flows [16–26].

    View all citing articles on Scopus
    View full text