Authors: Martin Losch, Jean-Michel Campin
:filelink:`pkg/shelfice` provides a thermodynamic model for basal melting underneath floating ice shelves.
CPP options enable or disable different aspects of the package (:numref:`shelfice_config`). Run-time options, flags, filenames and field-related dates/times are described in :numref:`shelfice_runtime`. A description of key subroutines is given in :numref:`shelfice_subroutines`. Available diagnostics output is listed in :numref:`shelfice_diagnostics`.
As with all MITgcm packages, :filelink:`pkg/shelfice` can be turned on or off at compile time:
- using the
packages.conffile by addingshelficeto it, - or using :filelink:`genmake2 <tools/genmake2>` adding
-enable=shelficeordisable=shelficeswitches
:filelink:`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 :filelink:`pkg/kpp`.
Parts of the :filelink:`pkg/shelfice` code can be enabled or disabled at compile time via CPP preprocessor flags. These options are set in :filelink:`SHELFICE_OPTIONS.h <pkg/shelfice/SHELFICE_OPTIONS.h>`:
.. tabularcolumns:: |\Y{.32}|\Y{.1}|\Y{.570}|
Compile-time parameters
| CPP Flag Name | Default | Description |
|---|---|---|
| :varlink:`ALLOW_SHELFICE_DEBUG` | #undef | include code for enhanced diagnostics and debug output |
| :varlink:`ALLOW_ISOMIP_TD` | #define | include code for for simplified ISOMIP thermodynamics |
| :varlink:`SHI_ALLOW_GAMMAFRICT` | #define | allow friction velocity-dependent transfer coefficient following Holland and Jenkins (1999) :cite:`holland:99` |
:filelink:`pkg/shelfice` is switched on/off at run time by setting :varlink:`useSHELFICE` to .TRUE. in file data.pkg.
Run-time parameters are set in file data.shelfice (read in :filelink:`pkg/shelfice/shelfice_readparms.F`),as listed below.
The data file specifying under-ice topography of ice shelves (:varlink:`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 :varlink:`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 :eq:`surfacepressure`. Note however the file :varlink:`SHELFICEloadAnomalyFile` must not be p_{top}, but p_{top}-g\sum_{k'=1}^{n-1}\rho_{0}\Delta{z}_{k'}, with \rho_{0}= :varlink:`rhoConst`, so that in the absence of a \rho^{*} that is different from \rho_{0}, the anomaly is zero.
.. tabularcolumns:: |\Y{.275}|\Y{.28}|\Y{.455}|
Run-time parameters and default values; all parameters are in namelist group SHELFICE_PARM01| Parameter | Default | Description |
|---|---|---|
| :varlink:`useISOMIPTD` | FALSE | use simplified ISOMIP thermodynamics on/off flag |
| :varlink:`SHELFICEconserve` | FALSE | use conservative form of temperature boundary conditions on/off flag |
| :varlink:`SHELFICEboundaryLayer` | FALSE | use simple boundary layer mixing parameterization on/off flag |
| :varlink:`SHI_withBL_realFWflux` | FALSE | with :varlink:`SHELFICEboundaryLayer`, allow to use real-FW flux |
| :varlink:`SHI_withBL_uStarTopDz` | FALSE | with :varlink:`SHELFICEboundaryLayer`, compute uStar from uVel,vVel averaged over top Dz thickness |
| :varlink:`SHELFICEloadAnomalyFile` | ' ' | initial geopotential anomaly |
| :varlink:`SHELFICEtopoFile` | ' ' | filename for under-ice topography of ice shelves |
| :varlink:`SHELFICEmassFile` | ' ' | filename for mass of ice shelves |
| :varlink:`SHELFICEMassDynTendFile` | ' ' | filename for mass tendency of ice shelves |
| :varlink:`SHELFICETransCoeffTFile` | ' ' | filename for spatially varying transfer coefficients |
| :varlink:`SHELFICElatentHeat` | 334.0E+03 | latent heat of fusion (J/kg) |
| :varlink:`SHELFICEHeatCapacity_Cp` | 2000.0E+00 | specific heat capacity of ice (J/kg/K) |
| :varlink:`rhoShelfIce` | 917.0E+00 | (constant) mean density of ice shelf (kg/m3) |
| :varlink:`SHELFICEheatTransCoeff` | 1.0E-04 | transfer coefficient (exchange velocity) for temperature (m/s) |
| :varlink:`SHELFICEsaltTransCoeff` | :varlink:`SHELFICEsaltToHeatRatio` * :varlink:`SHELFICEheatTransCoeff` | transfer coefficient (exchange velocity) for salinity (m/s) |
| :varlink:`SHELFICEsaltToHeatRatio` | 5.05E-03 | ratio of salinity to temperature transfer coefficients (non-dim.) |
| :varlink:`SHELFICEkappa` | 1.54E-06 | temperature diffusion coefficient of the ice shelf (m2/s) |
| :varlink:`SHELFICEthetaSurface` | -20.0E+00 | (constant) surface temperature above the ice shelf (oC) |
| :varlink:`no_slip_shelfice` | :varlink:`no_slip_bottom` | slip along bottom of ice shelf on/off flag |
| :varlink:`SHELFICEDragLinear` | :varlink:`bottomDragLinear` | linear drag coefficient at bottom ice shelf (m/s) |
| :varlink:`SHELFICEDragQuadratic` | :varlink:`bottomDragQuadratic` | quadratic drag coefficient at bottom ice shelf (non-dim.) |
| :varlink:`SHELFICEselectDragQuadr` | -1 | select form of quadratic drag coefficient (non-dim.) |
| :varlink:`SHELFICEMassStepping` | FALSE | recalculate ice shelf mass at every time step |
| :varlink:`SHELFICEDynMassOnly` | FALSE | if :varlink:`SHELFICEmassStepping` = TRUE, exclude freshwater flux contribution |
| :varlink:`SHELFICEadvDiffHeatFlux` | FALSE | use advective-diffusive heat flux into ice shelf instead of default diffusive heat flux |
| :varlink:`SHELFICEuseGammaFrict` | FALSE | use velocity dependent exchange coefficients (Holland and Jenkins 1999 :cite:`holland:99`) |
| :varlink:`SHELFICE_oldCalcUStar` | FALSE | use old uStar averaging expression |
| :varlink:`SHELFICEwriteState` | FALSE | write ice shelf state to file on/off flag |
| :varlink:`SHELFICE_dumpFreq` | :varlink:`dumpFreq` | dump frequency (s) |
| :varlink:`SHELFICE_dump_mnc` | :varlink:`snapshot_mnc` | write snapshot using MNC on/off flag |
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) :cite:`marshall: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}:
p_{tot} = p_{top} + \int_z^{\eta-h} g\,\rho\,dz + 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 :varlink:`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) :cite:`beckmann:99`), so that
p_{top} = p_{a} + \int_{-h}^{0}g\,\rho^{*}\,dz
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}. :eq:`pressureocean` becomes
p_{tot} = p_{top} - g \rho_0 (z+h) + g \rho_0 \eta + \, \int_z^{\eta-h}{ g (\rho-\rho_0) \, dz} + \, p_{NH}
and after rearranging
p'_{tot} = p'_{top} + 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 :eq:`ptop` from z=0 to the bottom of the last fully dry cell within the ice shelf:
p_{top} = g\,\sum_{k'=1}^{n-1}\rho_{k'}^{*}\Delta{z_{k'}} + p_{a}
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
p'_{k} = p'_{top} + g\rho_{n}\eta +
g\,\sum_{k'=n}^{k}\left((\rho_{k'}-\rho_{0})\Delta{z_{k'}}
\frac{1+H(k'-k)}{2}\right)
where H(k'-k)=1 for k'<k and 0 otherwise.
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 :varlink:`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}:
\overline{\theta}_{k} = \theta_{k} h_{k} + \theta_{k+1} (1-h_{k})
where h_{k}\in[0,1] is the fractional layer thickness of the k-th layer (see :numref:`shelfice_grid`). 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
g_{\theta,k} = \frac{Q}{\rho_{0} c_{p} h_{k} \Delta{z}_{k}}
\text{ and } g_{\theta,k+1} = 0
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
g_{\theta,k}^* = \frac{Q}{\rho_{0} c_{p} \Delta{z}_{k}}
g_{\theta,k+1}^*
= \frac{Q}{\rho_{0} c_{p} \Delta{z}_{k}}
\frac{ \Delta{z}_{k} ( 1- h_{k} )}{\Delta{z}_{k+1}}
:eq:`tendencykp1` 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) :cite:`jenkins:01`, 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:
{\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 :eq:`jenkinsbc` 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) :cite:`hellmer:89`, Jenkins et al. (2001) :cite:`jenkins:01`) the heat balance at the ice-ocean interface is expressed as
c_{p} \rho \gamma_T (T - T_{b})
+\rho_{I} c_{p,I} \kappa \frac{(T_{S} - T_{b})}{h} = -Lq
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) :cite:`holland:99`. 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
T_{f} = (0.0901 - 0.0575\ S_{b})^{\circ}
- 7.61 \times 10^{-4}\frac{\text{K}}{\text{dBar}}\ p_{b}
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:
\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). :eq:`hellmerheatbalance` to :eq:`hellmersaltbalance` can be solved for S_{b}, T_{b}, and the freshwater flux q due to melting. These values are substituted into expression :eq:`jenkinsbc` 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, :eq:`jenkinsbc` cannot be applied directly. In this case :eq:`hellmersaltbalance` can be used with :eq:`jenkinsbc` to obtain:
\rho{K}\frac{\partial{S}}{\partial{z}}\biggl|_{b} = q\,(S-S_I)
This formulation can be used for all cases for which :eq:`hellmersaltbalance` is valid. Further, in this formulation it is obvious that melting (q<0) leads to a reduction of salinity.
The default value of :varlink:`SHELFICEconserve` =.FALSE. removes the contribution q\, ( X_{b}-X ) from
:eq:`jenkinsbc`, 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 :eq:`hellmerfreeze` T_{b} = a_{0} S_{b} + \epsilon_{4} to eliminate T_{b} from :eq:`hellmerheatbalance` and find an expression for the freshwater flux q:
\begin{aligned}
-Lq &= \epsilon_{1} (T - a_{0} S_{b} - \epsilon_{4})
+ \epsilon_{3} (T_{S} - a_{0} S_{b} - \epsilon_{4}) \\
\Leftrightarrow Lq &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}
+ \epsilon_{q}
\end{aligned}
to be substituted into :eq:`hellmersaltbalance`:
\begin{aligned}
\epsilon_{2}\,(S - S_{b}) &= - Lq\,(S_{b}-S_{I})
= - (a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}
+ \epsilon_{q})\,(S_{b}-S_{I}) \\
\Leftrightarrow 0 &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}^{2}
+ \{ \epsilon_{q} - \epsilon_{2}
- a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{I} \}\,S_{b}
+ \epsilon_{2}\,S - \epsilon_{q}\,S_{I}
\end{aligned}
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
\begin{aligned}
0 &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}^{2}
+ (\epsilon_{q} - \epsilon_{2}) \,S_{b} + \epsilon_{2}\,S \\
S_{b} &= \frac{\epsilon_{2} - \epsilon_{q}\mp
\sqrt{(\epsilon_{q} - \epsilon_{2})^2
- 4\, a_{0}\,(\epsilon_{1} + \epsilon_{3})\,\epsilon_{2}}}
{2\,a_{0}\,(\epsilon_{1} + \epsilon_{3})}
\end{aligned}
With S_b, the boundary layer temperature T_b and the melt rate q are known through :eq:`hellmerfreeze` and :eq:`solvedmeltrate`.
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) :cite:`grosfeld:97`. In this formulation :eq:`hellmerheatbalance` reduces to
c_{p} \rho \gamma_T (T - T_{b}) = -Lq
and the fresh water flux q is computed from
q = - \frac{c_{p} \rho \gamma_T (T - T_{b})}{L}
In order to use this formulation, set run-time parameter :varlink:`useISOMIPTD` =.TRUE. in
data.shelfice.
The default exchange coefficents \gamma_{T/S} are constant and
set by the run-time parameters :varlink:`SHELFICEheatTransCoeff` and
:varlink:`SHELFICEsaltTransCoeff` (see
:numref:`tab_phys_pkg_shelfice_runtimeparms`). If
:varlink:`SHELFICEuseGammaFrict` =.TRUE., exchange coefficients
are computed from drag laws and friction velocities estimated from
ocean speeds following Holland and Jenkins (1999)
:cite:`holland:99`. This computation can be modified using run-time
parameters and user is referred to S/R
:filelink:`pkg/shelfice/shelfice_readparms.F` for details.
The shelfice package and experiments demonstrating its strengths and weaknesses are also described in Losch (2008) :cite:`losch:08`. Unfortunately however, the description of the thermodynamics in the appendix of Losch (2008) is wrong.
The main routine is :filelink:`shelfice_thermodynamics.F <pkg/shelfice/shelfice_thermodynamics.F>` but note that :filelink:`/pkg/shelfice` routines are also called when solving the momentum equations.
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
Diagnostics output is available via the diagnostics package (see :numref:`outp_pack`). Available output fields are summarized as follows:
---------+----+----+----------------+----------------- <-Name->|Levs|grid|<-- 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
See the verification experiment :filelink:`isomip <verification/isomip>` for example usage of :filelink:`pkg/shelfice`.