Home Contact Us Site Map  
 
       
    next up previous contents
Next: 6.6.4 STREAMICE Package Up: 6.6 Sea Ice Packages Previous: 6.6.2 SEAICE Package   Contents

Subsections


6.6.3 SHELFICE Package

Authors: Martin Losch, Jean-Michel Campin


6.6.3.1 Introduction

Package ``shelfice'' provides a thermodynamic model for basal melting underneath floating ice shelves.

CPP options enable or disable different aspects of the package (Section 6.6.3.2). Run-Time options, flags, filenames and field-related dates/times are set in data.shelfice (Section 6.6.3.3). A description of key subroutines is given in Section 6.6.3.5. Input fields, units and sign conventions are summarized in Section 6.6.3.3, and available diagnostics output is listed in Section 6.6.3.6.

6.6.3.2 SHELFICE configuration, compiling & running


6.6.3.2.1 Compile-time options

 

As with all MITgcm packages, 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
  • required packages and CPP options:
    SHELFICE does not require any additional packages, but it will only work with conventional vertical $ z$ -coordinates (pressure coordinates are not implemented, yet). If you use it together with vertical mixing schemes, be aware, that non-local parameterizations have been turned off, e.g. for KPP (6.4.2).
(see Section 3.4).

Parts of the SHELFICE code can be enabled or disabled at compile time via CPP preprocessor flags. These options are set SHELFICE_OPTIONS.h. Table 6.6.3.2 summarizes these options.


Table 6.21: Available CPP-flags to be set in SHELFICE_OPTIONS.h
CPP option Description
ALLOW_SHELFICE_DEBUG Include code for enhanged diagnosis
ALLOW_ISOMIP_TD Include code for simplifed ISOMIP thermodynamics



6.6.3.3 Run-time parameters

Run-time parameters are set in files data.pkg (read in packages_readparms.F), and data.shelfice (read in shelfice_readparms.F).

6.6.3.3.1 Enabling the package

 
A package is switched on/off at run-time by setting (e.g. for SHELFICE) useSHELFICE = .TRUE. in data.pkg.

6.6.3.3.2 General flags and parameters

 
Table 6.22 lists all run-time parameters.

Table 6.22: Run-time parameters and default values
Name Default value Description Reference
useISOMIPTD F use simplified ISOMIP thermodynamics instead of default
SHELFICEconserve F use conservative form of temperature boundary conditions
SHELFICEboundaryLayer F use simple boundary layer mixing parameterization
SHELFICEloadAnomalyFile UNSET inital geopotential anomaly
SHELFICEtopoFile UNSET under-ice topography of ice shelves
SHELFICElatentHeat 334.0E+03 latent heat of fusion ($ L$ )
SHELFICEHeatCapacity_Cp 2000.0E+00 specific heat capacity of ice ($ c_{p,I}$ )
rhoShelfIce 917.0E+00 (constant) mean density of ice shelf ($ \rho_{I}$ )
SHELFICEheatTransCoeff 1.0E-04 transfer coefficient (exchange velocity) for temperature ($ \gamma_T$ )
SHELFICEsaltTransCoeff 5.05E-03 $ \times $  SHELFICEheatTransCoeff transfer coefficient (exchange velocity) for salinity ($ \gamma_S$ )
SHELFICEkappa 1.54E-06 temperature diffusion coefficient of the ice shelf ($ \kappa$ )
SHELFICEthetaSurface -20.0E+00 (constant) surface temperature above the ice shelf ($ T_{S}$ )
no_slip_shelfice no_slip_bottom (data) flag for slip along bottom of ice shelf
SHELFICEDragLinear bottomDragLinear (data) linear drag coefficient at bottom ice shelf
SHELFICEDragQuadratic bottomDragQuadratic (data) quadratic drag coefficient at bottom ice shelf
SHELFICEwriteState F write ice shelf state to file
SHELFICE_dumpFreq dumpFreq (data) dump frequency
SHELFICE_taveFreq taveFreq (data) time-averaging frequency
SHELFICE_tave_mnc timeave_mnc (data.mnc) write snap-shot using MNC
SHELFICE_dump_mnc snapshot_mnc (data.mnc) write TimeAverage using MNC


6.6.3.3.3 Input fields and units

SHEFLICEtopoFile:
under-ice topography of ice shelves in meters; upwards is positive, that as for the bathymetry files, negative values are required for topography below the sea-level;
SHEFLICEloadAnomalyFile:
pressure load anomaly at the bottom of the ice shelves 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. (6.75); however, the file SHEFLICEloadAnomalyFile 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 absenses of a $ \rho^{*}$ that is different from $ \rho_{0}$ , the anomaly is zero.


6.6.3.4 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] 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}$ :

$\displaystyle p_{tot} = p_{top} + \int_z^{\eta-h} g\,\rho\,dz + p_{NH},$ (6.71)

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], so that

$\displaystyle p_{top} = p_{a} + \int_{-h}^{0}g\,\rho^{*}\,dz.$ (6.72)

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 the 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. (6.71) becomes

$\displaystyle p_{tot} =$ $\displaystyle \,p_{top} - g\,\rho_0\,(z+h)$ $\displaystyle + g\,\rho_0\,\eta + \int_z^{\eta-h} g\,(\rho-\rho_0)\,dz + p_{NH},$ (6.73)

and after rearranging


$\displaystyle p'_{tot} =$ $\displaystyle \,p'_{top}$ $\displaystyle + g\,\rho_0\,\eta + \int_z^{\eta-h} g\,(\rho-\rho_0)\,dz + p_{NH},$ (6.74)

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. (6.72) from $ z=0$ to the bottom of the last fully dry cell within the ice shelf:

$\displaystyle p_{top} = g\,\sum_{k'=1}^{n-1}\rho_{k'}^{*}\Delta{z_{k'}} + p_{a}$ (6.75)

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

$\displaystyle 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)$ (6.76)

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}$ :

$\displaystyle \overline{\theta}_{k} = \theta_{k} h_{k} + \theta_{k+1} (1-h_{k})$ (6.77)

where $ h_{k}\in(0,1]$ is the fractional layer thickness of the $ k$ -th layer. 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

$\displaystyle g_{\theta,k} = \frac{Q}{\rho_{0} c_{p} h_{k} \Delta{z}_{k}}$    and $\displaystyle g_{\theta,k+1} = 0$ (6.78)

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


Eq. (6.80) 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.


6.6.3.4.1 Three-Equations-Thermodynamics

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], 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:

$\displaystyle {\rho}K\frac{\partial{X}}{\partial{z}}\biggl\vert _{b} = (\rho\gamma_{X}-q) ( X_{b} - X )$ (6.81)

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. (6.81) 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; Jenkins et al., 2001] the heat balance at the ice-ocean interface is expressed as

$\displaystyle c_{p} \rho \gamma_T (T - T_{b}) +\rho_{I} c_{p,I} \kappa \frac{(T_{S} - T_{b})}{h} = -Lq$ (6.82)

where $ \rho $ is the density of sea-water, $ c_{p} = 3974$ J kg$^-1$ 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]. $ L = 334000$ J kg$^-1$ is the latent heat of fusion. $ \rho_{I} = 920$ kg m$^-3$ , $ c_{p,I} = 2000$ J kg$^-1$ K$^-1$ , and $ T_{S}$ are the density, heat capacity and the surface temperature of the ice shelf; $ \kappa=1.54\times10^{-6}$ m$^2$ 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}$ 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

$\displaystyle T_{f} = (0.0901 - 0.0575\ S_{b})^{\circ} - 7.61 \times 10^{-4}\frac{\text{K}}{\text{dBar}}\ p_{b}$ (6.83)

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:

$\displaystyle \rho \gamma_{S} (S - S_{b}) = - q\,(S_{b}-S_{I}),$ (6.84)

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$ ). Equations (6.82) to (6.84) can be solved for $ S_{b}$ , $ T_{b}$ , and the freshwater flux $ q$ due to melting. These values are substituted into expression (6.81) 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, equation (6.81) cannot be applied directly. In this case equation (6.84) can be used with Eq. (6.81) to obtain:

$\displaystyle \rho{K}\frac{\partial{S}}{\partial{z}}\biggl\vert _{b} = q\,(S-S_I).$ (6.85)

This formulation can be used for all cases for which equation (6.84) 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 Eq. (6.81), making the boundary conditions for temperature non-conservative.


6.6.3.4.2 ISOMIP-Thermodynamics

A simpler formulation follows the ISOMIP protocol (http://efdl.cims.nyu.edu/project_oisi/isomip/overview.html). The freezing and melting in the boundary layer between ice shelf and ocean is parameterized following Grosfeld et al. [1997]. In this formulation Eq. (6.82) reduces to

$\displaystyle c_{p} \rho \gamma_T (T - T_{b}) = -Lq$ (6.86)

and the fresh water flux $ q$ is computed from

$\displaystyle q = - \frac{c_{p} \rho \gamma_T (T - T_{b})}{L}.$ (6.87)

In order to use this formulation, set run-time parameter useISOMIPTD=.true. in data.shelfice.

6.6.3.4.3 Remark

The shelfice package and experiments demonstrating its strenghts and weaknesses are also described in Losch [2008]. However, note that unfortunately the description of the thermodynamics in the appendix of Losch [2008] is wrong.


6.6.3.5 Key subroutines

Top-level routine: shelfice_model.F 

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


6.6.3.6 SHELFICE diagnostics

Diagnostics output is available via the diagnostics package (see Section 7.1). Available output fields are summarized in Table 6.6.3.6.


Table 6.23: Available diagnostics of the shelfice-package
\begin{table}\centering {\footnotesize \begin{verbatim}---------+----+----+---... ...\vert Ice shelf forcing for salt, >0 increases salt\end{verbatim} }\end{table}



6.6.3.7 Experiments and tutorials that use shelfice

next up previous contents
Next: 6.6.4 STREAMICE Package Up: 6.6 Sea Ice Packages Previous: 6.6.2 SEAICE Package   Contents
mitgcm-support@mitgcm.org
Copyright © 2006 Massachusetts Institute of Technology Last update 2018-01-23