pkg/shelfice provides a thermodynamic model for basal melting
underneath floating ice shelves.
CPP options enable or disable different aspects of the package
(Section 8.6.3.2). Run-time options, flags, filenames and
field-related dates/times are described in Section 8.6.3.3. A description of key subroutines is given
in Section 8.6.3.5. Available diagnostics output is listed in
Section 8.6.3.6.
As with all MITgcm packages, pkg/shelfice can be turned on or off at compile
time:
using the packages.conf file by adding shelfice to it,
or using genmake2 adding -enable=shelfice or disable=shelfice switches
pkg/shelfice does not require any additional packages, but it will only
work with conventional vertical \(z\)-coordinates (pressure
coordinates are not implemented). If you use it together with
vertical mixing schemes, be aware that non-local parameterizations
are turned off, e.g., such as pkg/kpp.
Parts of the pkg/shelfice code can be enabled or disabled at compile time
via CPP preprocessor flags. These options are set in SHELFICE_OPTIONS.h:
The data file specifying under-ice topography of ice shelves (SHELFICEtopoFile) is in meters; upwards is positive,
and as for the bathymetry files, negative values are required for topography below the sea-level.
The data file for the pressure load anomaly at the bottom of the ice shelves SHELFICEloadAnomalyFile is in pressure
units (Pa). This field is absolutely required to avoid large
excursions of the free surface during initial adjustment processes,
obtained by integrating an approximate density from the surface at
\(z=0\) down to the bottom of the last fully dry cell within the
ice shelf, see (8.30). Note however the file SHELFICEloadAnomalyFile must
not be \(p_{top}\), but
\(p_{top}-g\sum_{k'=1}^{n-1}\rho_{0}\Delta{z}_{k'}\), with
\(\rho_{0}=\)rhoConst, so that in the absence of a \(\rho^{*}\)
that is different from \(\rho_{0}\), the anomaly is zero.
Table 8.21 Run-time parameters and default values; all parameters are in namelist group SHELFICE_PARM01¶
In the light of isomorphic equations for pressure and height
coordinates, the ice shelf topography on top of the water column has a
similar role as (and in the language of Marshall et al. (2004) [MAC+04],
is isomorphic to) the orography and the pressure boundary conditions at
the bottom of the fluid for atmospheric and oceanic models in pressure
coordinates. The total pressure \(p_{tot}\) in the ocean can be
divided into the pressure at the top of the water column
\(p_{top}\), the hydrostatic pressure and the non-hydrostatic
pressure contribution \(p_{NH}\):
with the gravitational acceleration \(g\), the density
\(\rho\), the vertical coordinate \(z\) (positive upwards), and
the dynamic sea-surface height \(\eta\). For the open ocean,
\(p_{top}=p_{a}\) (atmospheric pressure) and \(h=0\). Underneath
an ice-shelf that is assumed to be floating in isostatic equilibrium,
\(p_{top}\) at the top of the water column is the atmospheric
pressure \(p_{a}\) plus the weight of the ice-shelf. It is this
weight of the ice-shelf that has to be provided as a boundary condition
at the top of the water column (in run-time parameter SHELFICEloadAnomalyFile). The weight is
conveniently computed by integrating a density profile \(\rho^*\),
that is constant in time and corresponds to the sea-water replaced by
ice, from \(z=0\) to a “reference” ice-shelf draft at \(z=-h\) (Beckmann et al. (1999)
[BHT99]), so that
Underneath the ice shelf, the “sea-surface height” \(\eta\) is the
deviation from the “reference” ice-shelf draft \(h\). During a model
integration, \(\eta\) adjusts so that the isostatic equilibrium is
maintained for sufficiently slow and large scale motion.
In MITgcm, the total pressure anomaly \(p'_{tot}\) which is used
for pressure gradient computations is defined by substracting a purely
depth dependent contribution \(-g\rho_{0}z\) with a constant
reference density \(\rho_{0}\) from \(p_{tot}\).
(8.27) becomes
(8.29)¶\[p_{tot} = p_{top} - g \rho_0 (z+h) + g \rho_0 \eta + \, \int_z^{\eta-h}{ g (\rho-\rho_0) \, dz} + \, p_{NH}\]
with \(p'_{tot} = p_{tot} + g\,\rho_0\,z\) and
\(p'_{top} = p_{top} -
g\,\rho_0\,h\). The non-hydrostatic pressure contribution \(p_{NH}\)
is neglected in the following.
In practice, the ice shelf contribution to \(p_{top}\) is computed
by integrating (8.28) from \(z=0\) to the bottom of the
last fully dry cell within the ice shelf:
where \(n\) is the vertical index of the first (at least partially)
“wet” cell and \(\Delta{z_{k'}}\) is the thickness of the
\(k'\)-th layer (counting downwards). The pressure anomaly for
evaluating the pressure gradient is computed in the center of the “wet”
cell \(k\) as
where \(H(k'-k)=1\) for \(k'<k\) and \(0\) otherwise.
Figure 8.11 Schematic of a vertical section of the grid at the base of an ice shelf. Grid lines are thin;
the thick line is the model’s representation of the ice shelf-water interface. Plus signs mark the position
of pressure points for pressure gradient computations. The letters A, B, and C mark specific grid cells for
reference. \(h_k\) is the fractional cell thickness so that \(h_k \Delta z_k\) is the actual cell thickness.¶
Setting SHELFICEboundaryLayer=.TRUE. introduces a simple boundary layer that reduces the potential
noise problem at the cost of increased vertical mixing. For this purpose
the water temperature at the \(k\)-th layer abutting ice shelf
topography for use in the heat flux parameterizations is computed as a
mean temperature \(\overline{\theta}_{k}\) over a boundary layer of
the same thickness as the layer thickness \(\Delta{z}_{k}\):
where \(h_{k}\in[0,1]\) is the fractional layer thickness of the
\(k\)-th layer (see Figure 8.11). The original contributions due to ice shelf-ocean
interaction \(g_{\theta}\) to the total tendency terms
\(G_{\theta}\) in the time-stepping equation
\(\theta^{n+1} = f(\theta^{n},\Delta{t},G_{\theta}^{n})\) are
for layers \(k\) and \(k+1\) (\(c_{p}\) is the heat
capacity). Averaging these terms over a layer thickness
\(\Delta{z_{k}}\) (e.g., extending from the ice shelf base down to
the dashed line in cell C) and applying the averaged tendency to cell A
(in layer \(k\)) and to the appropriate fraction of cells C (in
layer \(k+1\)) yields
(8.35) describes averaging over the part of the grid
cell \(k+1\) that is part of the boundary layer with tendency
\(g_{\theta,k}^*\) and the part with no tendency. Salinity is
treated in the same way. The momentum equations are not modified.
Freezing and melting form a boundary layer between ice shelf and ocean.
Phase transitions at the boundary between saline water and ice imply the
following fluxes across the boundary: the freshwater mass flux \(q\)
(\(<0\) for melting); the heat flux that consists of the diffusive
flux through the ice, the latent heat flux due to melting and freezing
and the heat that is carried by the mass flux; and the salinity that is
carried by the mass flux, if the ice has a non-zero salinity
\(S_I\). Further, the position of the interface between ice and
ocean changes because of \(q\), so that, say, in the case of melting
the volume of sea water increases. As a consequence salinity and
temperature are modified.
The turbulent exchange terms for tracers at the ice-ocean interface are
generally expressed as diffusive fluxes. Following Jenkins et al. (2001)
[JHH01], the boundary conditions for a tracer
take into account that this boundary is not a material surface. The
implied upward freshwater flux \(q\) (in mass units, negative for
melting) is included in the boundary conditions for the temperature and
salinity equation as an advective flux:
(8.36)¶\[{\rho}K\frac{\partial{X}}{\partial{z}}\biggl|_{b}
= (\rho\gamma_{X}-q) ( X_{b} - X )\]
where tracer \(X\) stands for either temperature \(T\) or
salinity \(S\). \(X_b\) is the tracer at the interface (taken to
be at freezing), \(X\) is the tracer at the first interior grid
point, \(\rho\) is the density of seawater, and \(\gamma_X\) is
the turbulent exchange coefficient (in units of an exchange velocity).
The left hand side of (8.36) is shorthand for the
(downward) flux of tracer \(X\) across the boundary. \(T_b\),
\(S_b\) and the freshwater flux \(q\) are obtained from solving
a system of three equations that is derived from the heat and freshwater
balance at the ice ocean interface.
In this so-called three-equation-model (e.g., Hellmer and Olbers (1989)
[HO89], Jenkins et al. (2001) [JHH01])
the heat balance at the ice-ocean interface is expressed as
where \(\rho\) is the density of sea-water,
\(c_{p} = 3974 \, \text{J kg}^{-1} \text{K}^{-1}\) is the specific heat
capacity of water and \(\gamma_T\) the turbulent exchange
coefficient of temperature. The value of \(\gamma_T\) is discussed
in Holland and Jenkins (1999) [HJ99]. \(L =
334000 \, \text{J kg}^{-1}\) is the latent heat of fusion.
\(\rho_{I} = 920 \, \text{kg m}^{-3}\), \(c_{p,I} =
2000 \, \text{J kg}^{-1} \text{K}^{-1}\), and \(T_{S}\) are the density,
heat capacity and the surface temperature of the ice shelf;
\(\kappa=1.54\times10^{-6} \, \text{m}^2 \text{s}^{-1}\) is the heat
diffusivity through the ice-shelf and \(h\) is the ice-shelf draft.
The second term on the right hand side describes the heat flux through
the ice shelf. A constant surface temperature \(T_S=-20^{\circ}\text{C}\) is
imposed. \(T\) is the temperature of the model cell adjacent to the
ice-water interface. The temperature at the interface \(T_{b}\) is
assumed to be the in-situ freezing point temperature of sea-water
\(T_{f}\) which is computed from a linear equation of state
with the salinity \(S_{b}\) and the pressure \(p_{b}\) (in dBar)
in the cell at the ice-water interface. From the salt budget, the salt
flux across the shelf ice-ocean interface is equal to the salt flux due
to melting and freezing:
(8.39)¶\[\rho \gamma_{S} (S - S_{b}) = - q\,(S_{b}-S_{I})\]
where \(\gamma_S = 5.05\times10^{-3}\gamma_T\) is the turbulent
salinity exchange coefficient, and \(S\) and \(S_{b}\) are
defined in analogy to temperature as the salinity of the model cell
adjacent to the ice-water interface and at the interface, respectively.
Note, that the salinity of the ice shelf is generally neglected
(\(S_{I}=0\)). (8.37) to
(8.39) can be solved for \(S_{b}\),
\(T_{b}\), and the freshwater flux \(q\) due to melting. These
values are substituted into expression (8.36) to obtain the
boundary conditions for the temperature and salinity equations of the
ocean model. This formulation tends to yield smaller melt rates than the
simpler formulation of the ISOMIP protocol because the freshwater flux
due to melting decreases the salinity which raises the freezing point
temperature and thus leads to less melting at the interface. For a
simpler thermodynamics model where \(S_b\) is not computed
explicitly, for example as in the ISOMIP protocol,
(8.36) cannot be applied directly. In this case
(8.39) can be used with (8.36) to obtain:
This formulation can be used for all cases for which
(8.39) is valid. Further, in this
formulation it is obvious that melting (\(q<0\)) leads to a
reduction of salinity.
The default value of SHELFICEconserve=.FALSE. removes the contribution \(q\, ( X_{b}-X )\) from
(8.36), making the boundary conditions for temperature
non-conservative.
There has been some confusion about the three-equations system, so we document the solution in the code here:
We use (8.38)\(T_{b} = a_{0} S_{b} + \epsilon_{4}\) to eliminate \(T_{b}\)
from (8.37) and find an expression for the freshwater flux \(q\):
where the abbrevations \(\epsilon_{1} = c_{p} \rho \gamma_{T}\),
\(\epsilon_{2} = \rho L \gamma_{S}\), \(\epsilon_{3} =
\frac{\rho_{I} c_{p,I} \kappa}{h}\), \(\epsilon_{4}=b_{0}p +
c_{0}\), \(\epsilon_{q} = \epsilon_{1}\,(\epsilon_{4} - T) +
\epsilon_{3}\,(\epsilon_{4} - T_{S})\) have been introduced. The
quadratic equation in \(S_{b}\) is solved and the smaller
non-negative root is used. In the MITgcm code, the ice shelf salinity
\(S_{I}\) is always zero and the quadratic equation simplifies to
A simpler formulation follows the ISOMIP protocol. The
freezing and melting in the boundary layer between ice shelf and ocean
is parameterized following Grosfeld et al. (1997) [GGD97]. In this
formulation (8.37) reduces to
The default exchange coefficents \(\gamma_{T/S}\) are constant and
set by the run-time parameters SHELFICEheatTransCoeff and
SHELFICEsaltTransCoeff (see
Table 8.21). If
SHELFICEuseGammaFrict=.TRUE., exchange coefficients
are computed from drag laws and friction velocities estimated from
ocean speeds following Holland and Jenkins (1999)
[HJ99]. This computation can be modified using run-time
parameters and user is referred to S/R
pkg/shelfice/shelfice_readparms.F for details.
The shelfice package and experiments demonstrating its strengths and
weaknesses are also described in Losch (2008) [Los08]. Unfortunately however,
the description of the thermodynamics in the
appendix of Losch (2008) is wrong.
C !CALLING SEQUENCE:
C ...
C |-FORWARD_STEP :: Step forward a time-step ( AT LAST !!! )
C ...
C | |-DO_OCEANIC_PHY :: Control oceanic physics and parameterization
C ...
C | | |-SHELFICE_THERMODYNAMICS :: main routine for thermodynamics
C with diagnostics
C ...
C | |-THERMODYNAMICS :: theta, salt + tracer equations driver.
C ...
C | | |-EXTERNAL_FORCING_T :: Problem specific forcing for temperature.
C | | |-SHELFICE_FORCING_T :: apply heat fluxes from ice shelf model
C ...
C | | |-EXTERNAL_FORCING_S :: Problem specific forcing for salinity.
C | | |-SHELFICE_FORCING_S :: apply fresh water fluxes from ice shelf model
C ...
C | |-DYNAMICS :: Momentum equations driver.
C ...
C | | |-MOM_FLUXFORM :: Flux form mom eqn. package ( see
C ...
C | | | |-SHELFICE_U_DRAG :: apply drag along ice shelf to u-equation
C with diagnostics
C ...
C | | |-MOM_VECINV :: Vector invariant form mom eqn. package ( see
C ...
C | | | |-SHELFICE_V_DRAG :: apply drag along ice shelf to v-equation
C with diagnostics
C ...
C o
