Barotropic vorticity model

The barotropic vorticity model describes the evolution of a 2D non-divergent flow with velocity components   through self-advection, forces and dissipation. Due to the non-divergent nature of the flow, it can be described by (the vertical component) of the relative vorticity   .

The dynamical core presented here to solve the barotropic vorticity equations largely follows the idealized models with spectral dynamics developed at the Geophysical Fluid Dynamics Laboratory[^1]: A barotropic vorticity model[^2].

Many concepts of the Shallow water model and the Primitive equation model are similar, such that for example horizontal diffusion and the Time integration are only explained here.

Barotropic vorticity equation

The barotropic vorticity equation is the prognostic equation that describes the time evolution of relative vorticity with advection, Coriolis force, forcing and diffusion in a single global layer on the sphere.

We denote time, velocity vector  , Coriolis parameter , and hyperdiffusion   ( is the hyperdiffusion order, see Horizontal diffusion). We also include possible forcing terms   which act on the vorticity and/or on the zonal velocity and the meridional velocity and hence the curl   is a tendency for relative vorticity . See Extending SpeedyWeather how to define these.

Starting with some relative vorticity , the Laplacian is inverted to obtain the stream function

The zonal velocity and meridional velocity are then the (negative) meridional gradient and zonal gradient of

which is described in Derivatives in spherical coordinates. Using and we can then advect the absolute vorticity  . In order to avoid to calculate both the curl and the divergence of a flux we rewrite the barotropic vorticity equation as

with    the rotated velocity vector, because    . This is the form that is solved in the BarotropicModel, as outlined in the following section.

Algorithm

We briefly outline the algorithm that SpeedyWeather.jl uses in order to integrate the barotropic vorticity equation. As an initial step

0. Start with initial conditions of in spectral space and transform this model state to grid-point space:

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

• obtain zonal velocity through a Meridional derivative

• obtain meridional velocity through a Zonal derivative

• Transform zonal and meridional velocity , 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, add to to accumulate the tendency of

6. Compute the horizontal diffusion based on that tendency

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

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

9. Possibly do some output

10. Repeat from 1.

Horizontal diffusion

In SpeedyWeather.jl we use hyperdiffusion through an -th power Laplacian   (hyper when  ) which can be implemented as a multiplication of the spectral coefficients with    (see spectral Laplacian). It is therefore computationally not more expensive to apply hyperdiffusion over diffusion as the    can be precomputed. Note the sign change   here is such that the dissipative nature of the diffusion operator is retained for odd and even.

In SpeedyWeather.jl the diffusion is applied implicitly. For that, consider a leapfrog scheme with time step of variable to obtain from time steps   and , the next time step  

with being some tendency evaluated from . Now we want to add a diffusion term   with coefficient , which however, is implicitly calculated from , then

As the application of   is, for every spectral mode, equivalent to a multiplication of a constant, we can rewrite this to

and expand the numerator to

Hence the diffusion can be applied implicitly by updating the tendency as

which only depends on . Now let    be the explicit part and     the implicit part. Both parts can be precomputed and are    the implicit part. Both parts can be precomputed and are only an element-wise multiplication in spectral space. For every spectral harmonic we do

Hence 2 multiplications and 1 subtraction with precomputed constants. However, we will normalize the (hyper-)Laplacians as described in the following. This also will take care of the alternating sign such that the diffusion operation is dissipative regardless the power .

Normalization of diffusion

In physics, the Laplace operator is often used to represent diffusion due to viscosity in a fluid or diffusion that needs to be added to retain numerical stability. In both cases, the coefficient is of units and the full operator reads as with units  . This motivates us to normalize the Laplace operator by a constant of units and the coefficient by its inverse such that it becomes a damping timescale of unit . Given the application in spectral space we decide to normalize by the largest eigenvalue   such that all entries in the discrete spectral Laplace operator are in . This also has the effect that the alternating sign drops out, such that higher wavenumbers are always dampened and not amplified. The normalized coefficient    (always positive) is therefore reinterpreted as the (inverse) time scale at which the highest wavenumber is dampened to zero due to diffusion. Together we have

and the hyper-Laplacian of power follows as

and the implicit part is accordingly   . Note that the diffusion time scale is then also scaled by the radius, see next section.

Radius scaling

Similar to a non-dimensionalization of the equations, SpeedyWeather.jl scales the barotropic vorticity equation with to obtain normalized gradient operators as follows. A scaling for vorticity and stream function is used that is

This is also convenient as vorticity is often in the atmosphere, but the stream function more like and so this scaling brings both closer to 1 with a typical radius of the Earth of 6371km. The inversion of the Laplacians in order to obtain from therefore becomes

where the dimensionless gradients simply omit the scaling with ,  . The Barotropic vorticity equation scaled with is

with

•  , the scaled time

•  , the velocity vector (no scaling applied)

•  , the scaled Coriolis parameter

•  , the scaled forcing vector

•  , the scaled diffusion coefficient , which itself is normalized to a damping time scale, see Normalization of diffusion.

So scaling with the radius squared means we can use dimensionless operators, however, this comes at the cost of needing to deal with both a time step in seconds as well as a scaled time step in seconds per meter, which can be confusing. Furthermore, some constants like Coriolis or the diffusion coefficient need to be scaled too during initialization, which may be confusing too because values are not what users expect them to be. SpeedyWeather.jl follows the logic that the scaling to the prognostic variables is only applied just before the time integration and variables are unscaled for output and after the time integration finished. That way, the scaling is hidden as much as possible from the user. In hopefully many other cases it is clearly denoted that a variable or constant is scaled.

References

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

[^2]: Geophysical Fluid Dynamics Laboratory, The barotropic vorticity equation.