## Dr. Gavin A. Schmidt

### Using tracer moments within the GISS GCM

**Version 1.2 Dec 1999**

**Gavin Schmidt and Gary Russell**

## Introduction

The GISS GCM is unique in using prognostic first and second order moments as well as the mean for all the tracers (including water vapour and potential enthalpy) within each grid box. These were introduced to improve the advection of tracers, and have not as yet truly been used to their maximum potential. In theory, the quadratic (linear) moments provide sub-grid scale horizontal resolution three (twice) times more detailed than the coarser resolution of the U and V fields. That is to say that using the linear moments is (almost) equivalent to four boxes within each grid cell, while using the quadratic moments gives 9 boxes per grid cell. (Strictly speaking, the moments are not exactly equivalent to the higher resolution, information is lost going from 9 horizontal boxes to the 6 horizontal moments. However, this is a minor point). This takes 4x5 resolution down to 1.3x1.7 resolution for the tracers. In the vertical it can also double or triple the resolution.

This paper is intended to be a simple introduction to using the
moments more effectively in the GCM (i.e. not just in the dynamics).
We will use a generic concentration tracer `TR(I,J,L)` to describe
how the moments should be changed based on common processes in the
GCM. We should note that although water vapour `Q` and potential
enthalpy `T` are saved as concentrations, the tracers (such as
Radon, O-18 and the like) are generally saved as masses. This is
simply a unit issue and does not effect any of the following
advice. In the GCM, for quantities in mass units, the moments are
generally denoted with an `M` in the variable name (such as
`QXM`).

**Figure 1**. The quadratic and linear profiles within a gird box.

For each tracer there is the mean value `TR0`; the three first
order (linear) moments, `TRX, TRY,` and `TRZ`; and the 6
second order (quadratic) moments `TRXX, TRXY, TRYY, TRYZ,
TRZZ,` and `TRZX` (some older code uses the variable name
`TRXZ` instead). They are scaled so that they are generally the
same order of magnitude as the mean. A good way to visualize them in
the linear case (for a one dimensional tracer) is to think of the
value of the tracer at three points in the box (Figure 1). The left
hand most value is `TR-TRX`, the mid point is `TR` and
the right hand most value, `TR+TRX`. As a function of the
distance *x* from the mid point (scaled so that the box is of
length *2x* with midpoint *x*=0), the tracer has the
value

`TR`(*x*)= `TR0` + *x* `TRX`,
-1 < *x* < 1

Including the second order moments is also straightforward,

`TR`(*x*)= `TR0` + *x* `TRX` + (3/2
*x*^{2}-1/2)`TRXX`, -1 < *x* < 1

Note that the mid point value is now `TR0 - TRXX/2`. These
equations can of course be generalized in the three directions
(making sure to include the cross terms for the quadratic case).

`TR`(*x,y,z*)= `TR0` + *x* `TRX` +
*y* `TRY` + *z* `TRZ` + *xy*
`TRXY` + *yz* `TRYZ` + *zx* `TRZX`
+ (3/2 *x*^{2}-1/2) `TRXX`
+ (3/2 *y*^{2}-1/2) `TRYY`
+ (3/2 *z*^{2}-1/2) `TRZZ`

with -1 < *x,y,z *< 1. For mass tracers, each equation should be
divided by the air mass to get concentration.

There are two stages to using the moments more effectively. Firstly, by making sure that all operations on the mean fields are carried over consistently to the moments, and secondly, actually using the sub-grid scale resolution to perform more detailed physics within each grid box. The following sections outline the details of how this should be done.

## Rules for manipulating moments consistently

In this section we list the most common operations to the tracers and how the gradients should be manipulated consistently. The overriding principle is that operators should be applied to the whole of the profile (equation 3). Since it is linear, the operators can be applied to each of the moments independently. The secondary principle is that the profile should approximate the least squares quadratic fit to the tracer distribution, hence any change that is not uniform within a grid box can be incorporated by adjusting the profile accordingly.

- Complete mixing of two boxes:
- All gradients should also be mixed, except those in the direction
of the mixing which should be set to zero. For instance, during
dry convective mixing set
`TRZ=TRZZ=TRYZ=TRZX=0`in each box and each other moment to its mass-weighted average over both boxes. - Removing a fraction of tracer uniformly over a grid box:
- All moments should be multiplied by
*(1-a)*, where*a*is the fraction of the tracer that is removed. For instance, this should be done for radioactive decay, or condensation of water vapour. - Adding tracer uniformly over a grid box:
- This is somewhat more
subjective, but for stability the gradients should not be changed at
all. Otherwise, gradients can become unreasonably large and
destabilize the model. The only time this possibly should cause a
change for the gradients is when the process that is increasing the
tracer is a function of the tracer concentration itself (such as an
auto-catalytic reaction?). This asymmetry between adding and removing
tracer uniformly is slightly diffusive. In the Michael Prather's
chemistry models, the changes in moments of each species depend on the
moments of the other species involved in its production.
- Defining a vertically moving plume:
- All horizontal moments
should be preserved for the plume. All vertical moments (including
`TRYZ`and`TRZX`) and should be set to zero within the plume. - Removing a vertically moving plume:
- If the plume originates
uniformly from the box, all moments should be reduced by
*(1-a)*, where*a*is the fraction of the tracer that is removed. If the plume originates somewhere specific then the vertical moments should be re-calculated according to the least square fit (see`/HOSTS/Isis/u/cmglr/DOC/PLUMES.060777`for exact details on how this is done for the linear scheme). - Adding a vertically moving plume:
- If it is inserted uniformly,
all moments (in mass units) within the plume should be added to the
moments (again, in mass units) of the receiving grid box. Again, if
the plume is inserted in a specific point the vertical moments should
be re-calculated according to the least square fit (see
`/HOSTS/Isis/u/cmglr/DOC/PLUMES.060777`). - Diffusing tracer within a grid box:
- For diffusion in one
dimension (for instance vertically), the same operator that applies to
the mean should also apply to the moments in the unrelated
directions (i.e
`QX, QY, QXX, QYY, QXY`). The boundary conditions on these gradients should be no-flux at the surface (unless of course some sub-grid scale processes have been calculated using moments as well). For the moments involving`Z`, (i.e.`QZ, QZZ, QZX,`and`QYZ`) the equations for modifying them come from looking at weighted integrals of the exact diffusion equation within the box and depend on the values of the diffusion and the fluxes at the top and bottom of the box. See Appendix A for the relevant equations. - Using a separately calculated sub-gridscale profile:
- The PBL code and the radiation schemes (?) use vertical profiles
that are calculated separately from the moments. These should be
used to update the vertical moments so that the rest of the code
can use them. All that is needed is the least squares fit of
equation 2 in the vertical, to the internally calculated
profile. Possibly, the vertical profile should be used as to
adjust the initial values for the internal profiles (since other
processes may have changed the vertical moments since the last
time the internal profiles were calculated). See Appendix B for
how the equations for this.
- Adding a point source (usually at the surface):
- When tracer sources occur at the surface without any associated
mass flux, a delta function is the best approximation of the change. A
least squares fit to the new profile would then give the new
gradients. However, (depending on the amount of tracer present in the
box already relative to the addition), such a scheme would violate the
positive definite criteria. Hence the changes to the gradients in this
case must be modified. We choose to modify them based on two criteria:
firstly, that at no place does the implied tracer concentration
decrease, and secondly, that the changes are monotonic to the opposite
edge. With these constraints, for a point source of amount
`DTR`at the bottom surface (`z=-1`):

`TR0`=`TR0`+`DTR`;`TRZ`=`TRZ`- 3/2`DTR`;`TRZZ`=`TRZZ`+ 1/2`DTR`

`Q`), then the sub-gridscale profile also should be positive definite. However, this is not absolutely guaranteed by the dynamics. Applying the following limits (taken from Prather, 1986) will ensure this is the case, however, it does not ensure that the profile is smooth across grid-boxes. If these limits need to be applied, it is possible that the profile has already become very noisy and the numerics already compromised.

`TRX `=min(1.5 `TR0`,max(-1.5 `TR0`, `TRX`))

`TRXX`=min(2 `TR0`-| `TRX`|/3,max(| `TRX`|- `TR0`, `TRXX`))

`TRXY`=min( `TR0`,max(- `TR0`, `TRXY`))

`TRZX`=min( `TR0`,max(- `TR0`, `TRXZ`))

In order to check the fully resolved profile, the program
`ra:/u/exec/qusplt` will take as input lines of the form
`TR0 TRX TRXX` and plot the resulting profile. It will work
with either piped input, or read from a file. The output is displayed
on the screen and saved in x.ps (using `aplotX`). See file
`ra:/u/exec/qusplt.doc` for more info.

## Incorporating non-uniform sub-gridscale physics

There are two ways to incorporate processes that are not uniform over
the grid box. Firstly, we can breakdown the grid box into component
sub-boxes, or we can directly manipulate the tracer moments. It is
more straightforward to split the box into sub-boxes, but for column
physics and the like, this leads to a large increase in the time
required for the routines. Changing the gradients directly (as is done
in `DYNAM`) is more efficient, but also more complicated.

A note of caution: Although the tracer concentration profile can be
extended to the edges of any grid box, the errors in using the profile
increase with *x*, and the linear or quadratic profile is not
guaranteed to be continuous across the edge to the next box. Hence it
is more conservative to use values of *x* less than unit magnitude.
The C-grid model code generally uses *x*=+/-2/sqrt(12)=0.57735
since those are the points at which the linear least squares fit
crosses the quadratic profile. These points also have nice properties
for evaluating integrals of the tracer over the box. However, for
other cases it may be more straightforward to evenly split a box in
two (*x*=+/-1/2) or three (*x*=-2/3,0,2/3).

If we decide to split the box into sub-boxes, define the horizontal
sub-boxes using `IQ=1,3` and `JQ=1,3` and set `T(IQ,JQ)=TR(2*(IQ-2)/3, 2*(JQ-2)/3)) ` using equation 3 in two dimensions.
The physics can now be applied for each `(IQ,JQ)` pair.
Recalculating the moments is straightforward although there may be
round off errors to beware of (see
`/HOSTS/Isis/u/gavin/gissgcm/doc/sigfig.doc` for more
details). In the vertical a similar process can be used. For example,
testing for stability within a grid box can be done twice (between
*z*=-2/3 and *z*=0 and again between *z*=0 and
*z*=2/3). This could allow for more low level convective
cloud for instance.

Examples of manipulating the gradients directly (i.e. the Prather chemical schemes) depend heavily on the process involved and no examples will be given here.

## Conclusions

Even if much of the physics is calculated at the original resolution,
there is no excuse for not manipulating the gradients in a consistent
manner. Eventually, everywhere the tracers are modified some thought
must be given to the effect of the moments. Examples of how this works
can be found in the condensation and moist convection routines (
`CB263.S`), in `SURFCE` ( `SB354BM15`),
`DRYCNV` ( `PB354AM15`), and in `FILTER` (
`MB354AM15`).

It is possible that the noise introduced by the inconsistency in treating the gradients is beneficial to the model, increasing diffusion or triggering more convection. However, it does not seem a good idea to rely on these mechanisms. If more diffusion is wanted, then it should be explicitly entered in a controlled, and hopefully more physical, manner.

## Appendix A: Modifying the moments during diffusion

We will consider vertical diffusion with a non-dimensional *z*
(-1 < *z* < 1), and with *K* in units of per second. We
define *K'=(K _{1}+K_{-1})/2*,

*dK=(K*and similarly for the fluxes

_{1}-K_{-1})*F = K d*(however the gradient is calculated). Then if

`Q`/dz*dt*is the timestep, the equations for the moments are,

`QZX`_{+}=`QZX`_{-} - 3 dt K' `QZX`

`QYZ`_{+}=`QYZ`_{-} - 3 dt K' `QYZ`

`QZ`_{+}=`QZ`_{-} + dt (3F' - 3 K' `QZ` -
1.5 dK `QZZ`)

`QZZ`_{+}=`QZZ`_{-}+dt (2.5 dF - 2.5 dK
`QZ` - 15 K' `QZZ`)

Note that if there are asymmetries in the fluxes or diffusion coefficient, the first and second moments are entangled. If we assume an implicit calculation, the new moments are

`QZX`_{+}= `QZX`_{-} / (1 + 3 dt K'),
(and similarly for `QYZ`) and

`QZ`_{+} = [(1+15 dt K') ( `QZ`_{-}+3 dt F') -
1.5 dt dK ( `QZZ`_{-}+2.5 dt dF)]/ D

`QZZ`_{+}= [(1+3 dt K') ( `QZZ`_{-}+2.5 dt
dF) - 2.5 dt dK ( `QZ`_{-}+3 dt F')]/ D

where *D=(1 + 3 dt K')(1 + 15 dt K')- 3.75 dt ^{2}
dK^{2} *, (

*D > 0*if

*dt >= 0*and

*K >= 0*). This simplifies considerably in the case of constant diffusivity (

*dK = 0*).

In dimensional terms, every *K* should be replaced with
4 *K*/*dz*^{2} , every *F* with 2*F*/ *dz
* where *dz* is the height of the box in whatever unit is
appropriate (either m (height coordinates) for *K* in
m^{2}/s, or kg/m^{2} (mass coordinates) for *K*
in units kg^{2}/s m^{4}).

(Thanks to Max Kelley for helping correct this section).

## Appendix B: Quadratic least squares fit of an internal profile

This should be useful for using changes in internal profiles (such as
used in the PBL code) into changes to the vertical moments. For a
uniformly spaced internal vertical profile *tr(z _{i})*, define
departures from the mean

*tr'(z*. If the

_{i})*z*coordinate is scaled such that -1 <

*z*<1 over the grid box, the quadratic least squares fit is given by,

_{i}
`TRZ` = (T_{1}S_{22} - T_{2}
S_{12})/(S_{11} S_{22} -
S_{12}^{2})

`TRZZ`=(T_{1} S_{12} - T_{2} S_{11})/(S_{11} S_{22} - S_{12}^{2})

where

T_{1=}sum_{i} z_{i} tr'(z_{i}),
T_{2}=1/2 sum_{i} (1 - 3 z_{i}^{2})
tr'(z_{i})

S_{11}=sum_{i} z_{i}^{2},
S_{12}=1/2 sum_{i} z_{i} (1 - 3
z_{i}^{2}),

S_{22}=1/4 sum_{i} (1 - 3 z_{i}^{2})

If the *z _{i}* are not uniform then all the sums will
include extra weighting terms.

#### Simple case for half-box profiles

If the internal vertical profile is valid from the mid point down (i.e. over only half the box), and an even weighting between the upper and lower half of the boxes is used, the equations are very simple:
`TRZ ^{+}` = 3/2 (

`tz`+ 1/2

`TR0`+ 1/3

^{-}`TRZ`+ 1/8

^{-}`TRZZ`)

^{-}`TRZZ`= 5/2 (

^{+}`tzz`+ 1/8

`TRZ`+ 1/5

^{-}`TRZZ`)

^{-}
where `tz,tzz` are the integrals of `z*t(z)` and 1/2
`(3z ^{2}-1)*t(zz)` over the interval [-1,0].