Primitive equation model

The primitive equations are a hydrostatic approximation of the compressible Navier-Stokes equations for an ideal gas on a rotating sphere. We largely follow the idealised spectral dynamical core developed by GFDL[^GFDL1] and documented therein[^GFDL2].

The primitive equations solved by SpeedyWeather.jl for relative vorticity , divergence , logarithm of surface pressure , temperature and specific humidity are

with velocity  , rotated velocity   , Coriolis parameter , the Vertical advection operator, dry air gas constant , Virtual temperature , Geopotential , pressure and surface pressure , thermodynamic   with the heat capacity at constant pressure. Horizontal hyper diffusion of the form   with coefficient and power is added for every variable that is advected, meaning , but left out here for clarity, see Horizontal diffusion.

The parameterizations for the tendencies of from physical processes are denoted as   and are further described in the corresponding sections, see Parameterizations.

SpeedyWeather.jl implements a PrimitiveWet and a PrimitiveDry dynamical core. For a dry atmosphere, we have   and the virtual temperature   equals the temperature (often called absolute to distinguish from the virtual temperature). The terms in the primitive equations and their discretizations are discussed in the following sections.

Virtual temperature

In short: Virtual temperature

Virtual temperature is the temperature dry air would need to have to be as light as moist air. It is used in the dynamical core to include the effect of humidity on the density while replacing density through the ideal gas law with temperature.

We assume the atmosphere to be composed of two ideal gases: Dry air and water vapor. Given a specific humidity both gases mix, their pressures , ( for dry, for water vapor), and densities add in a given air parcel that has temperature . The ideal gas law then holds for both gases

with the respective specific gas constants   and   obtained from the universal gas constant divided by the molecular masses of the gas. The total pressure in the air parcel is

We ultimately want to replace the density    in the dynamical core, using the ideal gas law, with the temperature , so that we never have to calculate the density explicitly. However, in order to not deal with two densities (dry air and water vapor) we would like to replace temperature with a virtual temperature that includes the effect of humidity on the density. So, wherever we use the ideal gas law to replace density with temperature, we would use the virtual temperature, which is a function of the absolute temperature and specific humidity, instead. A higher specific humidity in an air parcel lowers the density as water vapor is lighter than dry air. Consequently, the virtual temperature of moist air is higher than its absolute temperature because warmer air is lighter too at constant pressure. We therefore think of the virtual temperature as the temperature dry air would need to have to be as light as moist air.

Starting with the last equation, with some manipulation we can write the ideal gas law as total density times a gas constant times the virtual temperature that is supposed to be a function of absolute temperature, humidity and some constants

Now we identify

as some constant that is positive for water vapor being lighter than dry air (  ) and

as the specific humidity. Given temperature and specific humidity , we can therefore calculate the virtual temperature as

For completeness we want to mention here that the above product, because it is a product of two variables has to be computed in grid-point space, see Spherical Harmonic Transform. To obtain an approximation to the virtual temperature in spectral space without expensive transforms one can linearize

with a global constant temperature , for example obtained from the    mode,   but depending on the normalization of the spherical harmonics that factor needs adjustment. We call this the linear virtual temperature which is used for the geopotential calculation, see #254.

Vertical coordinates

We start with some general considerations that apply when changing the vertical coordinate from height to something else. Let be some variable that depends on space and time. Now we want to express using some other coordinate in the vertical. Regardless of the coordinate system the value of at the to corresponding (and vice versa) has to be the same as we only want to change the coordinate, not itself.

So you can think of as a function of and as a function of . The chain rule lets us differentiate with respect to or

But for derivatives with respect to we have to apply the multi-variable chain-rule as both and depend on it. So a derivative with respect to on levels (where constant) becomes

So we first take the derivative of with respect to , but then also have to account for the fact that, at a given , depends on which is again dealt with using the univariate chain rule from above. We will make use of that for the Pressure gradient.

Sigma coordinates

The problem with pure pressure coordinates is that they are not terrain-following. For example, the 1000 hPa level in the Earth's atmosphere cuts through mountains. A flow field on such a level is therefore not continuous and one would need to deal with boundaries. Especially with spherical harmonics we need a terrain-following vertical coordinate to transform between continuous fields in grid-point space and spectral space.

SpeedyWeather.jl currently uses so-called sigma coordinates for the vertical. This coordinate system uses fraction of surface pressure in the vertical, i.e.

with   and   being the top (zero pressure) and   the surface (at surface pressure). As a consequence the vertical dimension is also indexed from top to surface.

Vertical indexing

Pressure, sigma, or hybrid coordinates in the vertical range from lowest values at the top to highest values at the surface. Consistently, we also index the vertical dimension top to surface. This means that   is the top-most layer, and   (or similar) is the layer that sits directly above the surface.

Sigma coordinates are therefore terrain-following, as   is always at surface pressure and so this level bends itself around every mountain, although the actual pressure on this level can vary. For a visualisation see #329.

One chooses levels associated with the -th layer and the pressure can be reobtained from the surface pressure

The layer thickness in terms of pressure is

which can also be expressed with the layer thickness in sigma coordinates times the surface pressure. In SpeedyWeather.jl one chooses the half levels first and then obtains the full levels through averaging

Geopotential

In the hydrostatic approximation the vertical momentum equation becomes

meaning that the (negative) vertical pressure gradient is given by the density in that layer times the gravitational acceleration. The heavier the fluid the more the pressure will increase below. Inserting the ideal gas law

with the geopotential   we can write this in terms of the logarithm of pressure

Note that we use the Virtual temperature here as we replaced the density through the ideal gas law with temperature. Given a vertical temperature profile and the (constant) surface geopotential   where is the orography, we can integrate this equation from the surface to the top to obtain on every layer . The surface is at    (see Vertical coordinates) with vertical levels. We can integrate the geopotential onto half levels as ( is the virtual temperature at layer , the subscript has been moved to be a superscript)

or onto full levels with

We use this last formula first to get from to , and then for every twice to get from to via . For the first half-level integration we use for the second .

Semi-implicit time integration: Geopotential

With the semi-implicit time integration in SpeedyWeather the Geopotential is not calculated from the spectral temperature at the current, but at the previous time step. This is because this is a linear term that we solve implicitly to avoid instabilities from gravity waves. For details see section Semi-implicit time stepping.

Surface pressure tendency

The surface pressure increases with a convergence of the flow above. Written in terms of the surface pressure directly, and not its logarithm

For discrete layers from 1 at the top to at the surface layer this can be written as

which can be thought of as a vertical integration of the pressure thickness-weighted divergence. In -coordinates with   (see Vertical coordinates) this becomes

Using the logarithm of pressure as the vertical coordinate this becomes

The second term is the divergence at layer . We introduce  , the -weighted vertical integration operator applied to some variable . This is essentially an average as  . The surface pressure tendency can then be written as

which is form used by SpeedyWeather.jl to calculate the tendency of (the logarithm of) surface pressure.

As we will have available in spectral space at the beginning of a time step, the gradient can be easily computed (see Derivatives in spherical coordinates). However, we then need to transform both gradients to grid-point space for the scalar product with the (vertically -averaged) velocity vector before transforming it back to spectral space where the tendency is needed. In general, we can do the -weighted average in spectral or in grid-point space, although it is computationally cheaper in spectral space. We therefore compute entirely in spectral space. With denoting spectral space and grid-point space (hence, and are the transforms in the respective directions) we therefore do

But note that it would also be possible to do

Meaning that we would compute the vertical average in grid-point space, subtract from the pressure gradient flux before transforming to spectral space. The same amount of transforms are performed but in the latter, the vertical averaging is done in grid-point space.

Semi-implicit time integration: Surface pressure tendency

With the semi-implicit time integration in SpeedyWeather the term is not evaluated from the spectral divergence at the current, but at the previous time step. This is because this is a linear term that we solve implicitly to avoid instabilities from gravity waves. For details see section Semi-implicit time stepping.

Vertical advection

The advection equation   for a tracer is, in flux form, for layer :

Starting from this equation in pressure layer , we can derive the advective form of the tracer transport. Dividing through by the layer thickness gives:

In sigma coordinates, the vertical mass flux can be expressed as  , where is the vertical velocity in sigma coordinates and is the surface pressure. Assuming is constant in time within the layer, and switching to sigma coordinates, we rewrite the equation as:

We now expand the horizontal divergence term using the product rule, thus:

From the continuity equation in sigma coordinates (mass conservation):

Substituting this into the equation:

Rearranging terms, we obtain:

With the reconstruction at the faces, and depending on one's choice of the advection scheme. For a second-order centered scheme, we choose    and obtain:

However, note that this scheme is dispersive and easily leads to instabilities at higher resolution, where a more advanced vertical advection scheme becomes necessary. For convenience, we may write to denote the vertical advection term , without specifying which schemes is used.

Vertical velocity

In the section Surface pressure tendency we used that the surface pressure changes with the convergence of the flow above, which derives from the conservation of mass. Similarly, the conservation of mass for layer can be expressed as (setting   in the advection equation in section Vertical advection)

Meaning that the pressure thickness of layer changes with a horizontal divergence   if not balanced by a net vertical mass flux into of the layer through the bottom and top boundaries of at  . is defined positive downward as this is the direction in which both pressure and sigma coordinates increase. The boundary conditions are   , such that there is no mass flux into the top layer from above or out of the surface layer and into the ground or ocean.

When integrating from the top down to layer we obtain the mass flux downwards out of layer

In sigma coordinates we have   with being the vertical velocity in sigma coordinates, also defined at interfaces between layers. To calculate we therefore compute

With denoting a sigma thickness-weighted vertical average as in section Surface pressure tendency. Now let be that average from   to   only and not necessarily down to the surface, as required in the equation above, then we can also write

See also Hoskins and Simmons, 1975[^HS75]. These vertical averages are the same as required by the Surface pressure tendency and in the Temperature equation, they are therefore all calculated at once, storing the partial averages and on the fly.

Pressure gradient

The pressure gradient term in the primitive equations is

with density and pressure . The gradient here is taken at constant hence the subscript. If we move to a pressure-based vertical coordinate system we will need to evaluate gradients on constant levels of pressure though, i.e. . There is, by definition, no gradient of pressure on constant levels of pressure, but we can use the chain rule (see Vertical coordinates) to rewrite this as (use only but is equivalent)

Using the hydrostatic equation    this becomes

Or, in terms of the geopotential  

which is the actual reason why we use pressure coordinates: As density also depends on the pressure the left-hand side means an implicit system when solving for pressure . To go from pressure to sigma coordinates we apply the chain rule from section Vertical coordinates again and obtain

where the last step inserts the hydrostatic equation again. With the ideal gas law, and note that we use Virtual temperature everywhere where the ideal gas law is used, but in combination with the dry gas constant

Combining the pressure in denominator and gradient to the logarithm and with   in Sigma coordinates (the logarithm of adds a constant that drops out in the gradient) we therefore have

From left to right: The pressure gradient force in -coordinates; in pressure coordinates; and in sigma coordinates. Each denoted with the respective subscript on gradients. SpeedyWeather.jl uses the latter. In sigma coordinates we may drop the subscript on gradients, but still meaning that the gradient is evaluated on a surface of our vertical coordinate. In vorticity-divergence formulation of the momentum equations the drops out in the vorticity equation (  ), but becomes a in the divergence equation, which is therefore combined with the kinetic energy term   similar as it is done in the Shallow water equations. You can think of    as the Bernoulli potential in the primitive equations. However, due to the change into sigma coordinates the surface pressure gradient also has to be accounted for. Now highlighting only the pressure gradient force, we have in total

In our vorticity-divergence formulation and with sigma coordinates.

Semi-implicit pressure gradient

With the semi-implicit time integration in SpeedyWeather.jl the pressure gradient terms are further modified as follows. See that section for details why, but here is just to mention that we need to split the terms into linear and non-linear terms. The linear terms are then evaluated at the previous time step for the implicit scheme such that we can avoid instabilities from gravity waves.

We split the (virtual) temperature into a reference vertical profile and its anomaly,   . The reference profile has to be a global constant for the spectral transform but can depend on the vertical. With this, the previous equation becomes

In the vorticity equation the term with the reference profile drops out as   , and in the divergence equation we move it into the Laplace operator. Now the linear terms are gathered with the Laplace operator and for the semi-implicit scheme we calculate both the Geopotential and the contribution to the "linear pressure gradient" at the previous time step for the semi-implicit time integration for details see therein.

Vorticity advection

Vorticity advection in the primitive equation takes the form

Meaning that we add the Coriolis parameter and the relative vorticity and multiply by the respective velocity component. While the primitive equations here are written with vorticity and divergence, we use here as other tendencies will be added and the curl and divergence are only taken once after transform into spectral space. To obtain a tendency for vorticity and divergence, we rewrite this as

with    the rotated velocity vector, see Barotropic vorticity equation.

Humidity equation

The dynamical core treats humidity as an (active) tracer, meaning that after the physical parameterizations for humidity are calculated in grid-point space, humidity is only advected with the flow. The only exception is the Virtual temperature as high levels of humidity will lower the effective density, which is why we use the virtual instead of the absolute temperature. The equation to be solved for humidity is therefore,

With denoting spectral space and grid-point space, so that and are the transforms in the respective directions. To avoid confusion with that notation, we write the tendency of humidity due to Vertical advection as . This equation is identical to a tracer equation, with denoting sources and sinks. Note that Horizontal diffusion should be applied to every advected variable.

A very similar equation is solved for (absolute) temperature as described in the following.

Temperature equation

The first law of thermodynamic states that the internal energy is increased by the heat applied minus the work done by the system. We neglect changes in chemical composition ([^Vallis], chapter 1.5). For an ideal gas, the internal energy is with the heat capacity at constant volume and temperature . The work done is , with pressure and the specific volume

For fluids we replace the differential here with the material derivative . With   and density we then have

Using the ideal gas law to replace with (we are using the Virtual temperature again), and using

we have

And further, with    the heat capacity at constant pressure,  , and using the logarithm of pressure

This is the form of the temperature equation that SpeedyWeather.jl uses. Temperature is advected through the material derivative and first term on the right-hand side represents an adiabatic conversion term describing how the temperature changes with changes in pressure. Recall that this term originated from the work term in the first law of thermodynamics. The forcing term is here identified as the physical parameterizations changing the temperature, for example radiation, and hence we will call it .

Similar to the Humidity equation we write the equation for (absolute) temperature as

is the Vertical advection of temperature. We evaluate the adiabatic conversion term completely in grid-point space following Simmons and Burridge, 1981[^SB81] Equation 3.12 and 3.13. Leaving out the for clarity, the term at level is

with

In sigma coordinates this simplifies to, following similar steps as in Surface pressure tendency

Let     and  , then this can also be summarised as

The are constants and can be precomputed. The surface pressure flux   has to be computed, so does the vertical sigma-weighted average from top to  , which is done when computing other vertical averages for the Surface pressure tendency.

Semi-implicit temperature equation

For the semi-implicit scheme we need to split the temperature equation into linear and non-linear terms, as the linear terms need to be evaluated at the previous time step. Decomposing temperature into    with the reference profile and its anomaly , the temperature equation becomes

Note that we do not change the adiabatic conversion term. While its linear component (the subscript for Virtual temperature as been raised) would need to be evaluated at the previous time step, we still evaluate this term at the current time step and move it within the semi-implicit corrections to the previous time step afterwards.

Semi-implicit time stepping

Conceptually, the semi-implicit time stepping in the Primitive equation model is the same as in the Shallow water model, but

• tendencies for divergence , logarithm of surface pressure but also temperature are computed semi-implicitly,

• the vertical layers are coupled, creating a linear equation system that is solved via matrix inversion.

The linear terms of the primitive equations follow a linearization around a state of rest without orography and a reference vertical temperature profile. The scheme described here largely follows Hoskins and Simmons [^HS75], which has also been used in Simmons and Burridge [^SB81].

As before, let   be the tendency we need for the Leapfrog time stepping. With the implicit time step  ,   we have

with being the explicitly-treated non-linear terms and the implicitly-treated linear terms, such that is a linear operator. We can therefore solve for by inverting ,

where we gathered the uncorrected right-hand side as

So for every linear term in we have two options corresponding to two sides of this equation

1. Evaluate it at the previous time step  

2. Or, evaluate it at the current time step as , but then move it back to the previous time step   by adding (in spectral space) the linear operator evaluated with the difference between the two time steps.

If there is a tendency that is easily evaluated in spectral space it is easier to follow 1. However, a term that is costly to evaluate in grid-point space should usually follow the latter. The reason is that the previous time step is generally not available in grid-point space (unless recalculated through a costly transform or stored with additional memory requirements) so it is easier to follow 2 where the is available in spectral space. For the adiabatic conversion term in the Temperature equation we follow 2 as one would otherwise need to split this term into a non-linear and linear term, evaluating it essentially twice in grid-point space.

So what is in the Primitive equation model?

is for the divergence, pressure and temperature equation the "uncorrected" tendency. Moving time step    we would be back with a fully explicit scheme. In the divergence equation the Geopotential is calculated from temperature at the previous time step   (denoted as superscript) and the "linear" Pressure gradient from the logarithm of surface pressure at the previous time step. One can think of these two calculations as linear operators, and . We will shortly discuss their properties. While we could combine them with the Laplace operator (which is also linear) we do not do this as do not depend on the degree and order of the spherical harmonics (their wavenumber) but on the vertical, but does not depend on the vertical, only on the wavenumber. All other terms are gathered in (subscript has been raised) and calculated as described in the respective section at the current time step .

For the pressure tendency, the subtraction with the thickness-weighted vertical average is the linear term that is treated implicitly. We call this operator . For the temperature tendency, we evaluate all terms explicitly at the current time step in but then move the linear term in the adiabatic conversion term with the operator back to the previous time step. For details see Semi-implicit temperature equation.

The operators are all linear, meaning that we can apply them in spectral space to each spherical harmonic independently – the vertical is coupled however. With being the number of vertical levels and the prognostic variables like temperature for a given degree and order being a column vector in the vertical,  , these operators have the following shapes

is an integration in the vertical hence it is an upper triangular matrix such that the first (an top-most)   element of the resulting vector depends on all vertical levels of the temperature mode , but the surface   only on the temperature mode at the surface. takes the surface value of the mode of the logarithm of surface pressure and multiplies it element-wise with the reference temperature profile and the dry gas constant. So the result is a column vector. is an   matrix as the adiabatic conversion term couples all layers. is a row vector as it represents the vertical averaging of the spherical harmonics of a divergence profile. So, is a scalar product for every giving a contribution of all vertical layers in divergence to the (single-layer!) logarithm of surface pressure tendency.

With the s defined we can now write the semi-implicit tendencies , , as (first equation in this section)

Solving for with the "combined" tendency

via

( is a matrix of size  ) yields

The other tendencies and are then obtained through insertion above. We may call the operator to be inverted which is of size   , hence for every degree of the spherical harmonics (which the Laplace operator depends on) a   matrix coupling the vertical levels. Furthermore, depends

• through on the time step ,

• through on the vertical level spacing

• through on the reference temperature profile

so for any changes of these the matrix inversion of has to be recomputed. Otherwise the algorithm for the semi-implicit scheme is as follows

0. Precompute the linear operators and with them the matrix inversion .

Then for every time step

1. Compute the uncorrected tendencies evaluated at the current time step for the explicit terms and the previous time step for the implicit terms.

2. Exception in SpeedyWeather.jl is the adiabatic conversion term, which is, using moved afterwards from the current to the previous time step  .

3. Compute the combined tendency from the uncorrected tendencies , , .

4. With the inverted operator get the corrected tendency for divergence,  .

5. Obtain the corrected tendencies for temperature and surface pressure from .

6. Apply Horizontal diffusion (which is only mentioned here as it further updates the tendencies).

7. Use , and in the Leapfrog time integration.

Horizontal diffusion

Horizontal diffusion in the primitive equations is applied to vorticity , divergence , temperature and humidity . In short, all variables that are advected. For the dry equations,   and no diffusion has to be applied.

The horizontal diffusion is applied implicitly in spectral space, as already described in Horizontal diffusion for the barotropic vorticity equation.

Algorithm

The following algorithm describes a time step of the PrimitiveWetModel, for the PrimitiveDryModel humidity can be set to zero and respective steps skipped.

0. Start with initial conditions of relative vorticity , divergence , temperature , humidity and the logarithm of surface pressure in spectral space. Variables are defined on all vertical levels, the logarithm of surface pressure only at the surface. Transform this model state to grid-point space, obtaining velocities is done as in the shallow water model

• 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

Additionally we

• Transform , , to in grid-point space

• Compute the (non-linearized) Virtual temperature in grid-point space.

Now loop over

1. Compute all tendencies of due to physical parameterizations in grid-point space.

2. Compute the gradient of the logarithm of surface pressure in spectral space and convert the two fields to grid-point space. Unscale the on the fly.

3. For every layer compute the pressure flux   in grid-point space.

4. For every layer compute a linearized Virtual temperature in spectral space.

5. For every layer compute a temperature anomaly (virtual and absolute) relative to a vertical reference profile in grid-point space.

6. Compute the Geopotential by integrating the virtual temperature vertically in spectral space from surface to top.

7. Integrate vertically to obtain in grid-point space and also in spectral space. Store on the fly also for every layer the partial integration from 1 to   (top to layer above). These will be used in the adiabatic term of the Temperature equation.

8. Compute the Surface pressure tendency with the vertical averages from the previous step. For the semi-implicit time stepping

9. For every layer compute the Vertical velocity.

10. For every layer add the linear contribution of the Pressure gradient to the geopotential in spectral space.

11. For every layer compute the Vertical advection for and add it to the respective tendency.

12. For every layer compute the tendency of due to Vorticity advection and the Pressure gradient and add to the respective existing tendency. Unscale , transform to spectral space, take curl and divergence to obtain tendencies for .

13. For every layer compute the adiabatic term and the horizontal advection in the Temperature equation in grid-point space, add to existing tendency and transform to spectral.

14. For every layer compute the horizontal advection of humidity in the Humidity equation in grid-point space, add to existing tendency and transform to spectral.

15. For every layer compute the kinetic energy  , transform to spectral and add to the Geopotential. For the semi-implicit time stepping also add the linear pressure gradient calculated from the previous time step. Now apply the Laplace operator and subtract from the divergence tendency.

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

17. Compute the horizontal diffusion for the advected variables

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

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

20. Possibly do some output

21. Repeat from 1.

Scaled primitive equations

TODO

References

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

[^GFDL2]: Geophysical Fluid Dynamics Laboratory, The Spectral Dynamical Core

[^Vallis]: GK Vallis, 2006. Atmopsheric and Ocean Fluid Dynamics, Cambridge University Press.

[^SB81]: Simmons and Burridge, 1981. An Energy and Angular-Momentum Conserving Vertical Finite-Difference Scheme and Hybrid Vertical Coordinates, Monthly Weather Review. DOI: 10.1175/1520-0493(1981)109<0758:AEAAMC>2.0.CO;2.

[^HS75]: Hoskins and Simmons, 1975. A multi-layer spectral model and the semi-implicit method, Quart. J. R. Met. Soc. DOI: 10.1002/qj.49710142918