# On the vertical coordinate in PISM, and a critical change of variable¶

In PISM all fields in the ice, including enthalpy, age, velocity, and so on, evolve within
an ice fluid domain of *changing geometry*. See Fig. 47. In
particular, the upper and lower surfaces of the ice fluid move with respect to the geoid.

Note

FIXME: This figure should show the floating case too.

The \((x,y,z)\) coordinates in Fig. 47 are supposed to be from an orthogonal coordinate system with \(z\) in the direction anti-parallel to gravity, so this is a flat-earth approximation. In practice, the data inputs to PISM are in some particular projection, of course.

We make a change of the independent variable \(z\) which simplifies how PISM deals with the changing geometry of the ice, especially in the cases of a non-flat or moving bed. We replace the vertical coordinate relative to the geoid with the vertical coordinate relative to the base of the ice. Let

where \(H = h - b\) is the ice thickness and \(\rho_i, \rho_o\) are densities of ice and ocean, respectively.

Now we make the change of variables

throughout the PISM code. This replaces \(z=b\) by \(s=0\) as the equation of the base surface of the ice. The ice fluid domain in the new coordinates only has a free upper surface as shown in Fig. 48.

Note

FIXME: This figure should show the floating case too, and bedrock.”

In PISM the computational domain (region)

is divided into a three-dimensional grid. See `IceGrid`

.

The change of variable \(z\to s\) used here *is not* the [170] change of variable
\(\tilde s=(z-b)/H\) . That change causes the conservation of energy equation to
become singular at the boundaries of the ice sheet. Specifically, the Jenssen change
replaces the vertical conduction term by a manifestly-singular term at ice sheet margins
where \(H\to 0\), because

A singular coefficient of this type can be assumed to affect the stability of all time-stepping schemes. The current change \(s=z-b\) has no such singularizing effect though the change does result in added advection terms in the conservation of energy equation, which we now address. See BOMBPROOF, PISM’s numerical scheme for conservation of energy for more general considerations about the conservation of energy equation.

The new coordinates \((x,y,s)\) are not orthogonal.

Recall that if \(f=f(x,y,z,t)\) is a function written in the old variables and if \(\tilde f(x,y,s,t)=f(x,y,z(x,y,s,t),t)\) is the “same” function written in the new variables, equivalently \(f(x,y,z,t)=\tilde f(x,y,s(x,y,z,t),t)\) , then

Similarly,

On the other hand,

The following table records some important changes to formulae related to conservation of energy:

Note \(E\) is the ice enthalpy and \(T\) is the ice temperature (which is a
function of the enthalpy; see `EnthalpyConverter`

), \(P\) is the ice pressure
(assumed hydrostatic), \(\mathbf{U}\) is the depth-dependent horizontal velocity, and
\(Q\) is the strain-heating term.

Now the vertical velocity is computed by
`StressBalance::compute_vertical_velocity(...)`

. In the old coordinates
\((x,y,z,t)\) it has this formula:

Here \(S\) is the basal melt rate, positive when ice is being melted. We have used the basal kinematical equation and integrated the incompressibility statement

In the new coordinates we have

(Note that the term \(\mathbf{U}(s) \cdot \nabla b\) evaluates the horizontal velocity at level \(s\) and not at the base.)

Let

This quantity is the vertical velocity of the ice *relative to the location on the bed
immediately below it*. In particular, \(\tilde w=0\) for a slab sliding down a
non-moving inclined plane at constant horizontal velocity, if there is no basal melt rate.
Also, \(\tilde w(s=0)\) is nonzero only if there is basal melting or freeze-on, i.e.
when \(S\ne 0\). Within PISM, \(\tilde w\) is written with name \(wvel_rel\) into an
input file. Comparing the last two equations, we see how
`StressBalance::compute_vertical_velocity(...)`

computes \(\tilde w\) :

The conservation of energy equation is now, in the new coordinate \(s\) and newly-defined relative vertical velocity,

Thus it looks just like the conservation of energy equation in the original vertical
velocity \(z\). This is the form of the equation solved by `EnthalpyModel`

using
`enthSystemCtx::solve()`

.

Under option `-o_size big`

, all of these vertical velocity fields are available as
fields in the output NetCDF file. The vertical velocity relative to the geoid, as a
three-dimensional field, is written as the diagnostic variable `wvel`

. This is the
“actual” vertical velocity \(w = \tilde w + \diff{b}{t} + \mathbf{U}(s)\cdot\nabla b\)
. Its surface value is written as `wvelsurf`

, and its basal value as `wvelbase`

. The
relative vertical velocity \(\tilde w\) is written to the NetCDF output file as
`wvel_rel`

.

Previous | Up | Next |