Shallow water model

The shallow water model describes the evolution of a 2D flow described by its velocity and an interface height that conceptually represents pressure. A divergent flow affects the interface height which in turn can impose a pressure gradient force onto the flow. The dynamics include advection, forces, dissipation, and continuity.

The following description of the shallow water model largely follows the idealized models with spectral dynamics developed at the Geophysical Fluid Dynamics Laboratory[^1]: The Shallow Water Equations[^2].

Shallow water equations

The shallow water equations of velocity   and interface height (i.e. the deviation from the fluid's rest height ) are, formulated in terms of relative vorticity   , divergence   

We denote time , Coriolis parameter , hyperdiffusion   ( is the hyperdiffusion order, see Horizontal diffusion), gravitational acceleration , dynamic layer thickness , and a forcing for the interface height . In the shallow water model the dynamics layer thickness is

that is, the layer thickness at rest plus the interface height minus orography . We also add various possible forcing terms   which can be defined to force the respective variables, see Extending SpeedyWeather.

In the shallow water system the flow can be described through or which are related through the stream function and the velocity potential (which is zero in the Barotropic vorticity equation).

With being the rotated gradient operator, in cartesian coordinates :   . See Derivatives in spherical coordinates for further details. Especially because the inversion of the Laplacian and the gradients of can be computed in a single pass, see U, V from vorticity and divergence.

The divergence/curl of the vorticity flux   are combined with the divergence/curl of the forcing vector , as

equivalently to how this is done in the Barotropic vorticity equation with   .

Algorithm

0. Start with initial conditions of relative vorticity , divergence , and interface height in spectral space and transform this model state to grid-point space:

• Invert the Laplacian of to obtain the stream function in spectral space

• Invert the Laplacian of to obtain the velocity potential in spectral space

• obtain velocities    from  

• Transform velocities , to grid-point space

• Unscale the factor to obtain

• Transform , , to in grid-point space

Now loop over

1. Compute the forcing (or drag) terms

2. Multiply with   in grid-point space

3. Add     and    

4. Transform these vector components to spectral space ,

5. Compute the curl of in spectral space and add to forcing of

6. Compute the divergence of in spectral space and add to forcing of divergence

7. Compute the kinetic energy   and transform to spectral space

8. Add to the kinetic energy the "geopotential" in spectral space to obtain the Bernoulli potential

9. Take the Laplacian of the Bernoulli potential and subtract from the divergence tendency

10. Compute the volume fluxes in grid-point space via    

11. Transform to spectral space and take the divergence for  . Add to forcing for .

12. Correct the tendencies following the semi-implicit time integration to prevent fast gravity waves from causing numerical instabilities

13. Compute the horizontal diffusion based on the tendencies

14. Compute a leapfrog time step as described in Time integration with a Robert-Asselin and Williams filter

15. Transform the new spectral state of , , to grid-point as described in 0.

16. Possibly do some output

17. Repeat from 1.

Semi-implicit time integration

Probably the biggest advantage of a spectral model is its ability to solve (parts of) the equations implicitly a low computational cost. The reason is that a linear operator can be easily inverted in spectral space, removing the necessity to solve large equation systems. An operation like   in grid-point space is costly because it requires a global communication, coupling all grid points. In spectral space is a diagonal operator, meaning that there is no communication between harmonics and its inversion is therefore easily done on a mode-by-mode basis of the harmonics.

This can be made use of when facing time stepping constraints with explicit schemes, where ridiculously small time steps to resolve fast waves would otherwise result in a horribly slow simulation. In the shallow water system there are gravity waves that propagate at a wave speed of (typically 300m/s), which, in order to not violate the CFL criterion for explicit time stepping, would need to be resolved. Therefore, treating the terms that are responsible for gravity waves implicitly would remove that time stepping constraint and allows us to run the simulation at the time step needed to resolve the advective motion of the atmosphere, which is usually one or two orders of magnitude longer than gravity waves.

In the following we will describe how the semi implicit time integration can be combined with the Leapfrog time stepping and the Robert-Asselin and Williams filter for a large increase in numerical stability with gravity waves. Let be the model state of all prognostic variables at time step , the leapfrog time stepping is then

with the right-hand side operator evaluated at the current time step . Now the idea is to split the terms in into non-linear terms that are evaluated explicitly in and into the linear terms , solved implicitly, that are responsible for the gravity waves. Linearization happens around a state of rest without orography.

We could already assume to evaluate at  , but in fact, we can introduce   so that for   we use   (i.e. explicit), for   it is centred implicit  , and for   a fully backwards scheme evaluated at  .

Let   be the tendency we need for the Leapfrog time stepping. Introducing   we have

because is a linear operator. This is done so that we can solve for by inverting , but let us gather the other terms as first.

For the shallow water equations we will only make use of the last formulation, meaning we first evaluate the whole right-hand side at the current time step as we would do with fully explicit time stepping but then add the implicit terms   afterwards to move those terms from to  . Note that we could also directly evaluate the implicit terms at   as it is suggested in the previous formulation  , the result would be the same. But in general it can be more efficient to do it one or the other way, and in fact it is also possible to combine both ways. This will be discussed in the semi-implicit time stepping for the primitive equations.

We can now implicitly solve for by

So what is ? In the shallow water system the gravity waves are caused by

which is a linearization of the equations around a state of rest with uniform constant layer thickness  . The continuity equation with the term, for example, is linearized to   . The divergence and continuity equations can now be written following the    formulation from above as a coupled system (The vorticity equation is zero for the linear gravity wave equation in the shallow water equations, hence no semi-implicit correction has to be made to the vorticity tendency).

with

Inserting the second equation into the first, we can first solve for , and then for . Reminder that we do this in spectral space to every harmonic independently, so the Laplace operator     takes the form of its eigenvalue   (normalized to unit sphere, as are the scaled shallow water equations) and its inversion is therefore just the inversion of this scalar.

Where the last formulation just makes it clear that    is the operator to be inverted. is then obtained via insertion as written above. Equivalently, by adding a superscript for every degree of the spherical harmonics, we have

The idea of the semi-implicit time stepping is now as follows:

1. Evaluate the right-hand side explicitly at time step to obtain the explicit, preliminary tendencies (and without a need for semi-implicit correction)

2. Move the implicit terms from to   when calculating

3. Solve for , the new, corrected tendency for divergence.

4. With obtain , the new, corrected tendency for .

5. Apply horizontal diffusion as a correction to as outlined in Horizontal diffusion.

6. Leapfrog with tendencies that have been corrected for both semi-implicit and diffusion.

Some notes on the semi-implicit time stepping

• The inversion of the semi-implicit time stepping depends on , that means every time the time step changes, the inversion has to be recalculated.

• You may choose   to dampen gravity waves but initialization shocks still usually kick off many gravity waves that propagate around the sphere for many days.

• With increasing   these waves are also slowed down, such that for   they quickly disappear in several hours.

• Using the scaled shallow water equations the time step has to be the scaled time step   which is divided by the radius . Then we use the normalized eigenvalues   which also omit the scaling, see scaled shallow water equations for more details.

Scaled shallow water equations

Similar to the scaled barotropic vorticity equations, SpeedyWeather.jl scales in the shallow water equations. The vorticity and the divergence equation are scaled with , the radius of the sphere squared, but the continuity equation is scaled with . We also combine the vorticity flux and forcing into a single divergence/curl operation as mentioned in Shallow water equations above

As in the scaled barotropic vorticity equations, one needs to scale the time step, the Coriolis force, the forcing and the diffusion coefficient, but then enjoys the luxury of working with dimensionless gradient operators. As before, SpeedyWeather.jl will scale vorticity and divergence just before the model integration starts and unscale them upon completion and for output. In the semi-implicit time integration we solve an equation that also has to be scaled. It is with radius squared scaling (because it is the tendency for the divergence equation which is also scaled with )

As is only scaled with we have

The normalizes the Laplace operator in the numerator, but using the scaled we also scale (which is convenient, because the time step within is the one we use anyway). The denominator does not actually change because   as is scaled with , but the Laplace operator with . So overall we just have to use the scaled time step and normalized eigenvalues for .

References

[^1]: Geophysical Fluid Dynamics Laboratory, Idealized models with spectral dynamics

[^2]: Geophysical Fluid Dynamics Laboratory, The Shallow Water Equations.