8.6.3. SHELFICE Package
Authors: Martin Losch, JeanMichel Campin
8.6.3.1. Introduction
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). Runtime options, flags, filenames and fieldrelated 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.
8.6.3.2. SHELFICE configuration
As with all MITgcm packages, pkg/shelfice can be turned on or off at compile time:
using the
packages.conf
file by addingshelfice
to it,or using genmake2 adding
enable=shelfice
ordisable=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 nonlocal 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:
CPP Flag Name 
Default 
Description 

#undef 
include code for enhanced diagnostics and debug output 

#define 
include code for for simplified ISOMIP thermodynamics 

#define 
allow friction velocitydependent transfer coefficient following Holland and Jenkins (1999) [HJ99] 
8.6.3.3. SHELFICE runtime parameters
pkg/shelfice is switched on/off at run time by setting useSHELFICE to .TRUE.
in file data.pkg
.
Runtime parameters are set in file data.shelfice
(read in pkg/shelfice/shelfice_readparms.F),as listed below.
The data file specifying underice 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 sealevel. 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.31). Note however the file SHELFICEloadAnomalyFile must not be \(p_{top}\), but \(p_{top}g\sum_{k'=1}^{n1}\rho_c \Delta{z}_{k'}\), with \(\rho_c =\) rhoConst, so that in the absence of a \(\rho^{*}\) that is different from \(\rho_c\), the anomaly is zero.
Parameter 
Default 
Description 

FALSE 
use simplified ISOMIP thermodynamics on/off flag 

FALSE 
use conservative form of temperature boundary conditions on/off flag 

FALSE 
use simple boundary layer mixing parameterization on/off flag 

FALSE 
with SHELFICEboundaryLayer, allow to use realFW flux 

FALSE 
with SHELFICEboundaryLayer, compute uStar from uVel,vVel averaged over top Dz thickness 

' ' 
initial geopotential anomaly 

' ' 
filename for underice topography of ice shelves 

' ' 
filename for mass of ice shelves 

' ' 
filename for mass tendency of ice shelves 

' ' 
filename for spatially varying transfer coefficients 

334.0E+03 
latent heat of fusion (J/kg) 

2000.0E+00 
specific heat capacity of ice (J/kg/K) 

917.0E+00 
(constant) mean density of ice shelf (kg/m^{3}) 

0.0E+00 
(constant) salinity of ice shelf 

1.0E04 
transfer coefficient (exchange velocity) for temperature (m/s) 

transfer coefficient (exchange velocity) for salinity (m/s) 

5.05E03 
ratio of salinity to temperature transfer coefficients (nondim.) 

1.54E06 
temperature diffusion coefficient of the ice shelf (m^{2}/s) 

20.0E+00 
(constant) surface temperature above the ice shelf (^{o}C) 

slip along bottom of ice shelf on/off flag 

linear drag coefficient at bottom ice shelf (m/s) 

quadratic drag coefficient at bottom ice shelf (nondim.) 

1 
select form of quadratic drag coefficient (nondim.) 

FALSE 
recalculate ice shelf mass at every time step 

FALSE 
if SHELFICEmassStepping = TRUE, exclude freshwater flux contribution 

FALSE 
use advectivediffusive heat flux into ice shelf instead of default diffusive heat flux 

FALSE 
use velocity dependent exchange coefficients (Holland and Jenkins 1999 [HJ99]) 

FALSE 
use old uStar averaging expression 

FALSE 
write ice shelf state to file on/off flag 

dump frequency (s) 

write snapshot using MNC on/off flag 
8.6.3.4. SHELFICE description
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_{\rm tot}\) in the ocean can be divided into the pressure at the top of the water column \(p_{\rm top}\), the hydrostatic pressure and the nonhydrostatic pressure contribution \(p_{\rm nh}\):
with the gravitational acceleration \(g\), the density \(\rho\), the vertical coordinate \(z\) (positive upwards), and the dynamic seasurface height \(\eta\). For the open ocean, \(p_{\rm top}=p_{a}\) (atmospheric pressure) and \(h=0\). Underneath an iceshelf that is assumed to be floating in isostatic equilibrium, \(p_{\rm top}\) at the top of the water column is the atmospheric pressure \(p_{a}\) plus the weight of the iceshelf. It is this weight of the iceshelf that has to be provided as a boundary condition at the top of the water column (in runtime parameter SHELFICEloadAnomalyFile). The weight is conveniently computed by integrating a density profile \(\rho^*\), that is constant in time and corresponds to the seawater replaced by ice, from \(z=0\) to a “reference” iceshelf draft at \(z=h\) (Beckmann et al. (1999) [BHT99]), so that
Underneath the ice shelf, the “seasurface height” \(\eta\) is the deviation from the “reference” iceshelf 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'_{\rm tot}\) which is used for pressure gradient computations is defined by subtracting a purely depth dependent contribution \(g\rho_c z\) using constant reference density \(\rho_c\) from \(p_{\rm tot}\). (8.28) becomes
and after rearranging
with \(p'_{\rm tot} = p_{\rm tot} + g\,\rho_c\,z\) and \(p'_{\rm top} = p_{\rm top}  g\,\rho_c\,h\). The nonhydrostatic pressure contribution \(p_{\rm nh}\) is neglected in the following.
In practice, the ice shelf contribution to \(p_{\rm top}\) is computed by integrating (8.29) 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.
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 shelfocean interaction \(g_{\theta}\) to the total tendency terms \(G_{\theta}\) in the timestepping 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.36) 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.
8.6.3.4.1. Threeequations thermodynamics
Freezing and melting form a boundary layer between the ice shelf and the ocean that is represented in the model by an infinitesimal layer used to calculate the exchanges between the ocean and the ice. Melting and freezing at the boundary between saline water and ice imply a freshwater mass flux \(q\) (\(<0\) for melting); the relevant heat fluxes into and out of the boundary layer therefore include a diffusive flux through the ice, the latent heat flux due to melting and freezing, and the advective heat that is carried by the freshwater mass flux. There is a salt flux carried by the mass flux if the ice has a nonzero 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 of these fluxes, salinity and temperature are affected.
The turbulent tracer exchanges between the infinitesimal boundary layer and the ocean are expressed as diffusive fluxes. Following Jenkins et al. (2001) [JHH01], the boundary conditions for a tracer take into account that this boundary may not move with the ice ocean interface (for example, in a linear free surface model). The implied upward freshwater flux \(q\) is therefore included in the boundary conditions for the temperature and salinity equation as an advective flux.
The boundary conditions for tracer \(X=S,T\) (tracer \(X\) stands for either insitu temperature \(T\) or salinity \(S\), located at the first interior ocean grid point) in the ocean are expressed as the sum of advective and diffusive fluxes
where the diffusive flux has been parameterized as a turbulent exchange \(\rho_c \, \gamma_{X}( X_{b}  X )\) following Holland and Jenkins (1999) [HJ99] or Jenkins et al. (2001) [JHH01]. \(X_b\) indicates the tracer in the boundary layer, \(\rho_c\) the density of seawater (parameter rhoConst), and \(\gamma_X\) is the turbulent exchange (or transfer) coefficient (parameters SHELFICEheatTransCoeff and SHELFICEsaltTransCoeff), in units of an exchange velocity. Insitu temperature, computed locally from tracer potential temperature, is required here to accurately compute the insitu freezing point of seawater in order to determine ice melt.
The tracer budget for the infinitesimal boundary layer takes the general form:
where the LHS represents diffusive flux from the ice evaluated at the interface between the infinitesimal boundary layer and the ice, and the RHS represents the turbulent and advective exchanges between the infinitesimal layer and the ocean and the advective exchange between the boundary layer and the ice (\(qX_{I}\), this flux will be zero if the ice contains no tracer: \(X_I=0\)). The temperature in the boundary layer (\(T_b\)) is taken to be at the freezing point, which is a function of pressure and salinity, \(\rho_I\) is ice density (rhoShelfIce), and \(K_{I,X}\) the appropriate ice diffusivity.
For any material tracer such as salinity, the LHS in (8.38) vanishes (no material diffusion into the ice), while for temperature, the term \(q\,( T_{b}T_{I} )\) vanishes, because both the boundary layer and the ice are at the freezing point. Instead, the latent heat of freezing is included as an additional term to take into account the conversion of ice to water:
where \(\kappa_{I,T}\) is the thermal ice diffusivity (SHELFICEkappa), \(L\) is the latent heat of fusion (SHELFICElatentHeat), \(c_{p}\) is the specific heat capacity of water (HeatCapacity_Cp), and \(c_{p,I}\) the heat capacity of the ice shelf (SHELFICEHeatCapacity_Cp). A reasonable choice for \(\gamma_T\) (SHELFICEheatTransCoeff), the turbulent exchange coefficient of temperature, is discussed in Holland and Jenkins (1999) [HJ99] (see Section 8.6.3.4.4). The temperature at the interface \(T_{b}\) is assumed to be the insitu freezing point temperature of seawater \(T_{f}\), which is computed from a linear equation of state:
where \(T_f\) is given in ^{o}C and \(p_{b}\) is in dBar. In (8.39), the diffusive heat flux at the iceocean interface can be appproximated by assuming a linear temperature profile in the ice and approximating the vertical derivative of temperature in the ice as the difference between the ice surface and ice bottom temperatures divided by the ice thickness, so that the lefthandside of (8.39) becomes
where \(T_{S}\) the (surface) temperature of the ice shelf (SHELFICEthetaSurface) and \(h\) is the iceshelf draft. Alternatively, assuming that the ice is “advected” vertically as implied by the meltflux \(q\), the diffusive flux can be approximated as \(\min(q,0)\,c_{p,I} (T_{S}  T_{b})\) (runtime flag SHELFICEadvDiffHeatFlux; see Holland and Jenkins, 1999 [HJ99] for details).
From the salt budget, the salt flux across the shelf iceocean interface is equal to the salt flux due to melting and freezing:
where \(\gamma_S =\) SHELFICEsaltToHeatRatio \(* \gamma_T\) is the turbulent salinity exchange coefficient. Note, it is assumed that \(\kappa_{I,S} =0\); moreover, the salinity of the ice shelf is generally neglected (\(S_{I}=0\)).
The budget equations for temperature (8.39) (with (8.41)) and salinity (8.42), together with the freezing point temperature of seawater (8.40), form the socalled threeequationmodel (e.g., Hellmer and Olbers (1989) [HO89], Jenkins et al. (2001) [JHH01]). These equations are solved to obtain \(S_b, T_b, q\), which are then substituted into (8.37) to obtain the boundary conditions for the temperature and salinity equations of the ocean model. Note that with \(S_{I}=0\) and (8.42), the boundary flux for salinity becomes \(F_S = q\,S\), which is the flux that is necessary to account for the dilution of salinity in the case of melting.
The threeequationmodel tends to yield smaller melt rates than the simpler formulation of the ISOMIP protocol (Section 8.6.3.4.3) because the freshwater flux (due to melting) decreases the salinity, which in turn 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.37) cannot be applied directly. In this case (8.42) can be used with (8.37) to obtain:
This formulation can be used for all cases for which (8.42) 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.37), making the
boundary conditions nonconservative.
8.6.3.4.2. Solving the threeequations system
There has been some confusion about the threeequations system, so we document the solution in the code here: We use (8.40) \(T_{b} = a_{0} S_{b} + \epsilon_{4}\) to eliminate \(T_{b}\) from (8.39) with (8.41) and find an expression for the freshwater flux \(q\):
to be substituted into (8.42):
with the abbrevations
The smaller nonnegative root of the quadratic equation in \(S_{b}\) is used. By default, the ice shelf salinity \(S_{I}\) is zero and the quadratic equation simplifies to
With \(S_b\), the boundary layer temperature \(T_b\) and the melt rate \(q\) are known through (8.40) and (8.43).
8.6.3.4.3. ISOMIP thermodynamics
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.39) reduces to
and the fresh water flux \(q\) is computed from
In order to use this formulation, set runtime parameter
useISOMIPTD =.TRUE.
in data.shelfice
.
8.6.3.4.4. Exchange coefficients
The default exchange coefficents \(\gamma_{T/S}\) are constant and
set by the runtime parameters SHELFICEheatTransCoeff and
SHELFICEsaltTransCoeff (see
Table 8.23). 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 runtime
parameters and the user is referred to
pkg/shelfice/shelfice_readparms.F for details.
8.6.3.4.5. Remark
The shelfice package and experiments demonstrating its strengths and weaknesses are also described in Losch (2008) [Los08]. Unfortunately, the description of the thermodynamics in the appendix of Losch (2008) is wrong.
8.6.3.5. Key subroutines
The main routine is shelfice_thermodynamics.F but note that /pkg/shelfice routines are also called when solving the momentum equations.
C !CALLING SEQUENCE:
C ...
C FORWARD_STEP :: Step forward a timestep ( 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 uequation
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 vequation
C with diagnostics
C ...
C o
8.6.3.6. SHELFICE diagnostics
Diagnostics output is available via the diagnostics package (see Section 9). Available output fields are summarized as follows:
++++
<Name>Levsgrid< Units >< Tile (max=80c)
++++
SHIfwFlx 1 SM kg/m^2/s Ice shelf fresh water flux (positive upward)
SHIhtFlx 1 SM W/m^2 Ice shelf heat flux (positive upward)
SHIUDrag 30 UU m/s^2 U momentum tendency from ice shelf drag
SHIVDrag 30 VV m/s^2 V momentum tendency from ice shelf drag
SHIForcT 1 SM W/m^2 Ice shelf forcing for theta, >0 increases theta
SHIForcS 1 SM g/m^2/s Ice shelf forcing for salt, >0 increases salt
8.6.3.7. Experiments and tutorials that use shelfice
See the verification experiment isomip for example usage of pkg/shelfice.