Warning, /doc/phys_pkgs/shelfice.rst is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 00f81e67 on 2024-05-24 21:00:12 UTC
b938a3c63b antn*0001 .. _sub_phys_pkg_shelfice:
0002
274298f12b Jeff*0003 SHELFICE Package
0004 ----------------
0005
0006 Authors: Martin Losch, Jean-Michel Campin
0007
0a8794f5ee Mart*0008 Introduction
274298f12b Jeff*0009 ~~~~~~~~~~~~
0010
0011 :filelink:`pkg/shelfice` provides a thermodynamic model for basal melting
0012 underneath floating ice shelves.
0013
0014 CPP options enable or disable different aspects of the package
0a8794f5ee Mart*0015 (:numref:`shelfice_config`). Run-time options, flags, filenames and
274298f12b Jeff*0016 field-related dates/times are described in :numref:`shelfice_runtime`. A description of key subroutines is given
0a8794f5ee Mart*0017 in :numref:`shelfice_subroutines`. Available diagnostics output is listed in
274298f12b Jeff*0018 :numref:`shelfice_diagnostics`.
0019
0020
0021 .. _shelfice_config:
0022
0023 SHELFICE configuration
0024 ~~~~~~~~~~~~~~~~~~~~~~
0025
0026 As with all MITgcm packages, :filelink:`pkg/shelfice` can be turned on or off at compile
0027 time:
0028
0029 - using the ``packages.conf`` file by adding ``shelfice`` to it,
0030
0031 - or using :filelink:`genmake2 <tools/genmake2>` adding ``-enable=shelfice`` or ``disable=shelfice`` switches
0032
0033 :filelink:`pkg/shelfice` does not require any additional packages, but it will only
0034 work with conventional vertical :math:`z`-coordinates (pressure
0035 coordinates are not implemented). If you use it together with
0036 vertical mixing schemes, be aware that non-local parameterizations
0037 are turned off, e.g., such as :filelink:`pkg/kpp`.
0038
0039 Parts of the :filelink:`pkg/shelfice` code can be enabled or disabled at compile time
0040 via CPP preprocessor flags. These options are set in :filelink:`SHELFICE_OPTIONS.h <pkg/shelfice/SHELFICE_OPTIONS.h>`:
0041
7b8b86ab99 Timo*0042 .. tabularcolumns:: |\Y{.32}|\Y{.1}|\Y{.570}|
0043
0044 .. table:: Compile-time parameters
0045 :name: tab_phys_pkg_shelfice_compileparms
0046
0047 +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
0048 | CPP Flag Name | Default | Description |
0049 +===============================================+=========+======================================================================================================================+
0050 | :varlink:`ALLOW_SHELFICE_DEBUG` | #undef | include code for enhanced diagnostics and debug output |
0051 +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
0052 | :varlink:`ALLOW_ISOMIP_TD` | #define | include code for for simplified ISOMIP thermodynamics |
0053 +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
0054 | :varlink:`SHI_ALLOW_GAMMAFRICT` | #define | allow friction velocity-dependent transfer coefficient following Holland and Jenkins (1999) :cite:`holland:99` |
0055 +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
274298f12b Jeff*0056
0057 .. _shelfice_runtime:
0058
0a8794f5ee Mart*0059 SHELFICE run-time parameters
274298f12b Jeff*0060 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0061
7b8b86ab99 Timo*0062 :filelink:`pkg/shelfice` is switched on/off at run time by setting :varlink:`useSHELFICE` to ``.TRUE.`` in file ``data.pkg``.
274298f12b Jeff*0063 Run-time parameters are set in file ``data.shelfice`` (read in :filelink:`pkg/shelfice/shelfice_readparms.F`),as listed below.
0064
0065 The data file specifying under-ice topography of ice shelves (:varlink:`SHELFICEtopoFile`) is in meters; upwards is positive,
0066 and as for the bathymetry files, negative values are required for topography below the sea-level.
0067 The data file for the pressure load anomaly at the bottom of the ice shelves :varlink:`SHELFICEloadAnomalyFile` is in pressure
0068 units (Pa). This field is absolutely required to avoid large
0a8794f5ee Mart*0069 excursions of the free surface during initial adjustment processes,
274298f12b Jeff*0070 obtained by integrating an approximate density from the surface at
0071 :math:`z=0` down to the bottom of the last fully dry cell within the
0072 ice shelf, see :eq:`surfacepressure`. Note however the file :varlink:`SHELFICEloadAnomalyFile` must
0073 not be :math:`p_{top}`, but
67dae37801 Eli *0074 :math:`p_{top}-g\sum_{k'=1}^{n-1}\rho_c \Delta{z}_{k'}`, with
0075 :math:`\rho_c =` :varlink:`rhoConst`, so that in the absence of a :math:`\rho^{*}`
0076 that is different from :math:`\rho_c`, the anomaly is zero.
274298f12b Jeff*0077
7b8b86ab99 Timo*0078 .. tabularcolumns:: |\Y{.275}|\Y{.28}|\Y{.455}|
274298f12b Jeff*0079
7b8b86ab99 Timo*0080 .. table:: Run-time parameters and default values; all parameters are in namelist group ``SHELFICE_PARM01``
0a8794f5ee Mart*0081 :name: tab_phys_pkg_shelfice_runtimeparms
0082 :class: longtable
0083
7b8b86ab99 Timo*0084 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0085 | Parameter | Default | Description |
0086 +========================================+============================================+=========================================================================================================+
0087 | :varlink:`useISOMIPTD` | FALSE | use simplified ISOMIP thermodynamics on/off flag |
0088 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0089 | :varlink:`SHELFICEconserve` | FALSE | use conservative form of temperature boundary conditions on/off flag |
0090 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0091 | :varlink:`SHELFICEboundaryLayer` | FALSE | use simple boundary layer mixing parameterization on/off flag |
0092 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0093 | :varlink:`SHI_withBL_realFWflux` | FALSE | with :varlink:`SHELFICEboundaryLayer`, allow to use real-FW flux |
0094 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0095 | :varlink:`SHI_withBL_uStarTopDz` | FALSE | with :varlink:`SHELFICEboundaryLayer`, compute uStar from uVel,vVel averaged over top Dz thickness |
0096 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0097 | :varlink:`SHELFICEloadAnomalyFile` | :kbd:`' '` | initial geopotential anomaly |
0098 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0099 | :varlink:`SHELFICEtopoFile` | :kbd:`' '` | filename for under-ice topography of ice shelves |
0100 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0101 | :varlink:`SHELFICEmassFile` | :kbd:`' '` | filename for mass of ice shelves |
0102 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0103 | :varlink:`SHELFICEMassDynTendFile` | :kbd:`' '` | filename for mass tendency of ice shelves |
0104 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0105 | :varlink:`SHELFICETransCoeffTFile` | :kbd:`' '` | filename for spatially varying transfer coefficients |
0106 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0107 | :varlink:`SHELFICElatentHeat` | 334.0E+03 | latent heat of fusion (J/kg) |
0108 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0109 | :varlink:`SHELFICEHeatCapacity_Cp` | 2000.0E+00 | specific heat capacity of ice (J/kg/K) |
0110 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0111 | :varlink:`rhoShelfIce` | 917.0E+00 | (constant) mean density of ice shelf (kg/m\ :sup:`3`) |
0112 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
00f81e6785 Ou W*0113 | :varlink:`SHELFICEsalinity` | 0.0E+00 | (constant) salinity of ice shelf |
0114 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
7b8b86ab99 Timo*0115 | :varlink:`SHELFICEheatTransCoeff` | 1.0E-04 | transfer coefficient (exchange velocity) for temperature (m/s) |
0116 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0117 | :varlink:`SHELFICEsaltTransCoeff` | :varlink:`SHELFICEsaltToHeatRatio` * | transfer coefficient (exchange velocity) for salinity (m/s) |
0118 | | :varlink:`SHELFICEheatTransCoeff` | |
0119 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0120 | :varlink:`SHELFICEsaltToHeatRatio` | 5.05E-03 | ratio of salinity to temperature transfer coefficients (non-dim.) |
0121 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0122 | :varlink:`SHELFICEkappa` | 1.54E-06 | temperature diffusion coefficient of the ice shelf (m\ :sup:`2`\ /s) |
0123 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0124 | :varlink:`SHELFICEthetaSurface` | -20.0E+00 | (constant) surface temperature above the ice shelf (:sup:`o`\ C) |
0125 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0126 | :varlink:`no_slip_shelfice` | :varlink:`no_slip_bottom` | slip along bottom of ice shelf on/off flag |
0127 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0128 | :varlink:`SHELFICEDragLinear` | :varlink:`bottomDragLinear` | linear drag coefficient at bottom ice shelf (m/s) |
0129 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0130 | :varlink:`SHELFICEDragQuadratic` | :varlink:`bottomDragQuadratic` | quadratic drag coefficient at bottom ice shelf (non-dim.) |
0131 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0132 | :varlink:`SHELFICEselectDragQuadr` | -1 | select form of quadratic drag coefficient (non-dim.) |
0133 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0134 | :varlink:`SHELFICEMassStepping` | FALSE | recalculate ice shelf mass at every time step |
0135 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0136 | :varlink:`SHELFICEDynMassOnly` | FALSE | if :varlink:`SHELFICEmassStepping` = TRUE, exclude freshwater flux contribution |
0137 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0138 | :varlink:`SHELFICEadvDiffHeatFlux` | FALSE | use advective-diffusive heat flux into ice shelf instead of default diffusive heat flux |
0139 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0140 | :varlink:`SHELFICEuseGammaFrict` | FALSE | use velocity dependent exchange coefficients (Holland and Jenkins 1999 :cite:`holland:99`) |
0141 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0142 | :varlink:`SHELFICE_oldCalcUStar` | FALSE | use old uStar averaging expression |
0143 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0144 | :varlink:`SHELFICEwriteState` | FALSE | write ice shelf state to file on/off flag |
0145 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0146 | :varlink:`SHELFICE_dumpFreq` | :varlink:`dumpFreq` | dump frequency (s) |
0147 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0148 | :varlink:`SHELFICE_dump_mnc` | :varlink:`snapshot_mnc` | write snapshot using MNC on/off flag |
0149 +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0a8794f5ee Mart*0150
0151 SHELFICE description
274298f12b Jeff*0152 ~~~~~~~~~~~~~~~~~~~~
0153
0154 In the light of isomorphic equations for pressure and height
0155 coordinates, the ice shelf topography on top of the water column has a
bd7056e9aa Jeff*0156 similar role as (and in the language of Marshall et al. (2004) :cite:`marshall:04`,
274298f12b Jeff*0157 is isomorphic to) the orography and the pressure boundary conditions at
0158 the bottom of the fluid for atmospheric and oceanic models in pressure
0bad585a21 Navi*0159 coordinates. The total pressure :math:`p_{\rm tot}` in the ocean can be
274298f12b Jeff*0160 divided into the pressure at the top of the water column
0bad585a21 Navi*0161 :math:`p_{\rm top}`, the hydrostatic pressure and the non-hydrostatic
0162 pressure contribution :math:`p_{\rm nh}`:
274298f12b Jeff*0163
0164 .. math::
0bad585a21 Navi*0165 p_{\rm tot} = p_{\rm top} + \int_z^{\eta-h} g\,\rho\,dz + p_{\rm nh}
274298f12b Jeff*0166 :label: pressureocean
0167
0168
0169 with the gravitational acceleration :math:`g`, the density
0170 :math:`\rho`, the vertical coordinate :math:`z` (positive upwards), and
0171 the dynamic sea-surface height :math:`\eta`. For the open ocean,
0bad585a21 Navi*0172 :math:`p_{\rm top}=p_{a}` (atmospheric pressure) and :math:`h=0`. Underneath
274298f12b Jeff*0173 an ice-shelf that is assumed to be floating in isostatic equilibrium,
0bad585a21 Navi*0174 :math:`p_{\rm top}` at the top of the water column is the atmospheric
274298f12b Jeff*0175 pressure :math:`p_{a}` plus the weight of the ice-shelf. It is this
0176 weight of the ice-shelf that has to be provided as a boundary condition
0177 at the top of the water column (in run-time parameter :varlink:`SHELFICEloadAnomalyFile`). The weight is
0178 conveniently computed by integrating a density profile :math:`\rho^*`,
0179 that is constant in time and corresponds to the sea-water replaced by
0180 ice, from :math:`z=0` to a “reference” ice-shelf draft at :math:`z=-h` (Beckmann et al. (1999)
0181 :cite:`beckmann:99`), so that
0182
0183 .. math::
0bad585a21 Navi*0184 p_{\rm top} = p_{a} + \int_{-h}^{0}g\,\rho^{*}\,dz
274298f12b Jeff*0185 :label: ptop
0186
0187 Underneath the ice shelf, the “sea-surface height” :math:`\eta` is the
0188 deviation from the “reference” ice-shelf draft :math:`h`. During a model
0189 integration, :math:`\eta` adjusts so that the isostatic equilibrium is
0190 maintained for sufficiently slow and large scale motion.
0191
0bad585a21 Navi*0192 In MITgcm, the total pressure anomaly :math:`p'_{\rm tot}` which is used
67dae37801 Eli *0193 for pressure gradient computations is defined by subtracting a purely
0194 depth dependent contribution :math:`-g\rho_c z` using constant
0bad585a21 Navi*0195 reference density :math:`\rho_c` from :math:`p_{\rm tot}`.
274298f12b Jeff*0196 :eq:`pressureocean` becomes
0197
0198 .. math::
0bad585a21 Navi*0199 p_{\rm tot} = p_{\rm top} - g \rho_c (z+h) + g \rho_c \eta + \, \int_z^{\eta-h}{ g (\rho-\rho_c) \, dz} + \, p_{\rm nh}
274298f12b Jeff*0200 :label: pressure
0201
0202 and after rearranging
0203
0204 .. math::
0bad585a21 Navi*0205 p'_{\rm tot} = p'_{\rm top} + g \rho_c \eta + \, \int_z^{\eta-h}{g (\rho-\rho_c) \, dz} + \, p_{\rm nh}
274298f12b Jeff*0206
0bad585a21 Navi*0207 with :math:`p'_{\rm tot} = p_{\rm tot} + g\,\rho_c\,z` and
0208 :math:`p'_{\rm top} = p_{\rm top} - g\,\rho_c\,h`.
0209 The non-hydrostatic pressure contribution :math:`p_{\rm nh}`
274298f12b Jeff*0210 is neglected in the following.
0211
0bad585a21 Navi*0212 In practice, the ice shelf contribution to :math:`p_{\rm top}` is computed
274298f12b Jeff*0213 by integrating :eq:`ptop` from :math:`z=0` to the bottom of the
0214 last fully dry cell within the ice shelf:
0215
0216 .. math::
0bad585a21 Navi*0217 p_{\rm top} = g\,\sum_{k'=1}^{n-1}\rho_{k'}^{*}\Delta{z_{k'}} + p_{a}
274298f12b Jeff*0218 :label: surfacepressure
0219
0220 where :math:`n` is the vertical index of the first (at least partially)
0221 “wet” cell and :math:`\Delta{z_{k'}}` is the thickness of the
0222 :math:`k'`-th layer (counting downwards). The pressure anomaly for
0223 evaluating the pressure gradient is computed in the center of the “wet”
0224 cell :math:`k` as
0225
0226 .. math::
0bad585a21 Navi*0227 p'_{k} = p'_{\rm top} + g\rho_{n}\eta +
67dae37801 Eli *0228 g\,\sum_{k'=n}^{k}\left((\rho_{k'}-\rho_c)\Delta{z_{k'}}
274298f12b Jeff*0229 \frac{1+H(k'-k)}{2}\right)
0230 :label: discretizedpressure
0231
0232 where :math:`H(k'-k)=1` for :math:`k'<k` and :math:`0` otherwise.
0233
bd7056e9aa Jeff*0234 .. figure:: figs/gridschematic.*
0235 :width: 80%
0236 :align: center
0237 :alt: schematic of vertical section of grid
0238 :name: shelfice_grid
0239
0a8794f5ee Mart*0240 Schematic of a vertical section of the grid at the base of an ice shelf. Grid lines are thin;
0241 the thick line is the model’s representation of the ice shelf-water interface. Plus signs mark the position
0242 of pressure points for pressure gradient computations. The letters A, B, and C mark specific grid cells for
0243 reference. :math:`h_k` is the fractional cell thickness so that :math:`h_k \Delta z_k` is the actual cell thickness.
bd7056e9aa Jeff*0244
0245
274298f12b Jeff*0246 Setting :varlink:`SHELFICEboundaryLayer` ``=.TRUE.`` introduces a simple boundary layer that reduces the potential
0247 noise problem at the cost of increased vertical mixing. For this purpose
0248 the water temperature at the :math:`k`-th layer abutting ice shelf
0249 topography for use in the heat flux parameterizations is computed as a
0250 mean temperature :math:`\overline{\theta}_{k}` over a boundary layer of
0251 the same thickness as the layer thickness :math:`\Delta{z}_{k}`:
0252
0253 .. math::
0254 \overline{\theta}_{k} = \theta_{k} h_{k} + \theta_{k+1} (1-h_{k})
0255 :label: thetabl
0256
0257 where :math:`h_{k}\in[0,1]` is the fractional layer thickness of the
bd7056e9aa Jeff*0258 :math:`k`-th layer (see :numref:`shelfice_grid`). The original contributions due to ice shelf-ocean
274298f12b Jeff*0259 interaction :math:`g_{\theta}` to the total tendency terms
0260 :math:`G_{\theta}` in the time-stepping equation
0261 :math:`\theta^{n+1} = f(\theta^{n},\Delta{t},G_{\theta}^{n})` are
0262
0263 .. math::
67dae37801 Eli *0264 g_{\theta,k} = \frac{Q}{\rho_c c_{p} h_{k} \Delta{z}_{k}}
274298f12b Jeff*0265 \text{ and } g_{\theta,k+1} = 0
0266 :label: orgtendency
0267
0268 for layers :math:`k` and :math:`k+1` (:math:`c_{p}` is the heat
0269 capacity). Averaging these terms over a layer thickness
0270 :math:`\Delta{z_{k}}` (e.g., extending from the ice shelf base down to
0271 the dashed line in cell C) and applying the averaged tendency to cell A
0272 (in layer :math:`k`) and to the appropriate fraction of cells C (in
0273 layer :math:`k+1`) yields
0274
0275 .. math::
67dae37801 Eli *0276 g_{\theta,k}^* = \frac{Q}{\rho_c c_{p} \Delta{z}_{k}}
274298f12b Jeff*0277 :label: tendencyk
0278
0279 .. math::
0280 g_{\theta,k+1}^*
67dae37801 Eli *0281 = \frac{Q}{\rho_c c_{p} \Delta{z}_{k}}
274298f12b Jeff*0282 \frac{ \Delta{z}_{k} ( 1- h_{k} )}{\Delta{z}_{k+1}}
0283 :label: tendencykp1
0284
0285 :eq:`tendencykp1` describes averaging over the part of the grid
0286 cell :math:`k+1` that is part of the boundary layer with tendency
0287 :math:`g_{\theta,k}^*` and the part with no tendency. Salinity is
0288 treated in the same way. The momentum equations are not modified.
0289
0a8794f5ee Mart*0290 Three-equations thermodynamics
274298f12b Jeff*0291 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0292
67dae37801 Eli *0293 Freezing and melting form a boundary layer between the ice shelf and the
0294 ocean that is represented in the model by an infinitesimal layer used to
0295 calculate the exchanges between the ocean and the ice. Melting and
0296 freezing at the boundary between saline water and ice imply a freshwater
0297 mass flux :math:`q` (:math:`<0` for melting); the relevant heat fluxes
0298 into and out of the boundary layer therefore include a diffusive flux
0299 through the ice, the latent heat flux due to melting and freezing, and
0300 the advective heat that is carried by the freshwater mass flux. There is
0301 a salt flux carried by the mass flux if the ice has a non-zero salinity
274298f12b Jeff*0302 :math:`S_I`. Further, the position of the interface between ice and
0303 ocean changes because of :math:`q`, so that, say, in the case of melting
67dae37801 Eli *0304 the volume of sea water increases. As a consequence of these fluxes,
0305 salinity and temperature are affected.
0306
0307 The turbulent tracer exchanges between the infinitesimal boundary layer and
0308 the ocean are expressed as diffusive fluxes. Following Jenkins et
0309 al. (2001) :cite:`jenkins:01`, the boundary conditions for a tracer take
0310 into account that this boundary may not move with the ice ocean interface
0311 (for example, in a linear free surface model). The implied upward freshwater
0312 flux :math:`q` is therefore included in the boundary conditions for the
0313 temperature and salinity equation as an advective flux.
0314
0315 The boundary conditions for tracer :math:`X=S,T` (tracer :math:`X` stands
0316 for either in-situ temperature :math:`T` or salinity :math:`S`, located at
0317 the first interior ocean grid point) in the ocean are expressed as the sum
0318 of advective and diffusive fluxes
274298f12b Jeff*0319
0320 .. math::
67dae37801 Eli *0321 F_X = (\rho_c \, \gamma_{X} -q ) ( X_{b} - X )
274298f12b Jeff*0322 :label: jenkinsbc
0323
67dae37801 Eli *0324 where the diffusive flux has been parameterized as a turbulent exchange
0325 :math:`\rho_c \, \gamma_{X}( X_{b} - X )` following Holland and
0326 Jenkins (1999) :cite:`holland:99` or Jenkins et al. (2001)
0327 :cite:`jenkins:01`. :math:`X_b` indicates the tracer in the boundary layer,
0328 :math:`\rho_c` the density of seawater (parameter :varlink:`rhoConst`), and
0329 :math:`\gamma_X` is the turbulent exchange (or transfer) coefficient
0330 (parameters :varlink:`SHELFICEheatTransCoeff` and
0331 :varlink:`SHELFICEsaltTransCoeff`), in units of an exchange
0332 velocity. In-situ temperature, computed locally from tracer potential
0333 temperature, is required here to accurately compute the in-situ freezing
0334 point of seawater in order to determine ice melt.
0335
0336 The tracer budget for the infinitesimal boundary layer takes the general
0337 form:
0338
0339 .. math::
0340 {\rho_I} K_{I,X} \frac{\partial{X_I}}{\partial{z}}\biggl|_{b}
0341 = \rho_c \, \gamma_{X} ( X_{b} - X ) - q ( X_{b} - X_{I} )
0342 :label: jenkinsgenbudget
0343
0344 where the LHS represents diffusive flux from the ice evaluated at the
0345 interface between the infinitesimal boundary layer and the ice, and the RHS
0346 represents the turbulent and advective exchanges between the infinitesimal
0347 layer and the ocean and the advective exchange between the boundary layer
0348 and the ice (:math:`qX_{I}`, this flux will be zero if the ice contains no
0349 tracer: :math:`X_I=0`). The temperature in the boundary layer (:math:`T_b`)
0350 is taken to be at the freezing point, which is a function of pressure and
0351 salinity, :math:`\rho_I` is ice density (:varlink:`rhoShelfIce`), and
0352 :math:`K_{I,X}` the appropriate ice diffusivity.
0353
0354 For any material tracer such as salinity, the LHS in :eq:`jenkinsgenbudget`
0355 vanishes (no material diffusion into the ice), while for temperature, the
0356 term :math:`q\,( T_{b}-T_{I} )` vanishes, because both the boundary layer
0357 and the ice are at the freezing point. Instead, the latent heat of freezing
0358 is included as an additional term to take into account the conversion of
0359 ice to water:
274298f12b Jeff*0360
0361 .. math::
67dae37801 Eli *0362 {\rho_I} \, c_{p,I} \, \kappa_{I,T}
0363 \frac{\partial{T_I}}{\partial{z}}\biggl|_{b}
0364 = c_{p} \, \rho_c \, \gamma_{T} ( T_{b} - T )+ L q.
0365 :label: jenkinsheatbudget
0366
0367 where :math:`\kappa_{I,T}` is the thermal ice diffusivity
0368 (:varlink:`SHELFICEkappa`), :math:`L` is the latent heat of fusion
0369 (:varlink:`SHELFICElatentHeat`), :math:`c_{p}` is the specific heat
0370 capacity of water (:varlink:`HeatCapacity_Cp`), and :math:`c_{p,I}` the
0371 heat capacity of the ice shelf (:varlink:`SHELFICEHeatCapacity_Cp`). A
0372 reasonable choice for :math:`\gamma_T` (:varlink:`SHELFICEheatTransCoeff`),
0373 the turbulent exchange coefficient of temperature, is discussed in Holland
0374 and Jenkins (1999) :cite:`holland:99` (see
0375 :numref:`shelfice_exchange_coeff`). The temperature at the interface
0376 :math:`T_{b}` is assumed to be the in-situ freezing point temperature of
0377 sea-water :math:`T_{f}`, which is computed from a linear equation of state:
274298f12b Jeff*0378
0379 .. math::
67dae37801 Eli *0380 T_{f} = (0.0901 - 0.0575\ S_{b})
0381 - 7.61 \times 10^{-4} p_{b}.
0a8794f5ee Mart*0382 :label: hellmerfreeze
274298f12b Jeff*0383
67dae37801 Eli *0384 where :math:`T_f` is given in :sup:`o`\ C and :math:`p_{b}` is in dBar. In
0385 :eq:`jenkinsheatbudget`, the diffusive heat flux at the ice-ocean interface
0386 can be appproximated by assuming a linear temperature profile in the ice
0387 and approximating the vertical derivative of temperature in the ice as the
0388 difference between the ice surface and ice bottom temperatures divided by
0389 the ice thickness, so that the left-hand-side of :eq:`jenkinsheatbudget`
0390 becomes
274298f12b Jeff*0391
0392 .. math::
67dae37801 Eli *0393 {\rho_I} \, c_{p,I} \, \kappa_{I,T}
0394 \frac{\partial{T_I}}{\partial{z}}\biggl|_{b}
0395 \approx \rho_{I} \, c_{p,I} \, \kappa_{I,T} \frac{(T_{S} - T_{b})}{h}
0396 :label: dTdzdiffus
0397
0398 where :math:`T_{S}` the (surface) temperature of the ice shelf
0399 (:varlink:`SHELFICEthetaSurface`) and :math:`h` is the ice-shelf
0400 draft. Alternatively, assuming that the ice is "advected" vertically as
0401 implied by the meltflux :math:`q`, the diffusive flux can be approximated
0402 as :math:`\min(q,0)\,c_{p,I} (T_{S} - T_{b})` (runtime flag
0403 :varlink:`SHELFICEadvDiffHeatFlux`; see Holland and Jenkins, 1999
0404 :cite:`holland:99` for details).
0405
0406 From the salt budget, the salt flux across the shelf ice-ocean interface is
0407 equal to the salt flux due to melting and freezing:
0408
0409 .. math::
0410 \rho_c \, \gamma_{S} (S - S_{b}) = - q\,(S_{b}-S_{I})
274298f12b Jeff*0411 :label: hellmersaltbalance
0412
67dae37801 Eli *0413 where :math:`\gamma_S =` :varlink:`SHELFICEsaltToHeatRatio` :math:`*
0414 \gamma_T` is the turbulent salinity exchange coefficient. Note, it is
0415 assumed that :math:`\kappa_{I,S} =0`; moreover, the salinity of the ice
0416 shelf is generally neglected (:math:`S_{I}=0`).
0417
c717173819 Matt*0418 The budget equations for temperature :eq:`jenkinsheatbudget` (with
0419 :eq:`dTdzdiffus`) and salinity
67dae37801 Eli *0420 :eq:`hellmersaltbalance`, together with the freezing point temperature of
0421 sea-water :eq:`hellmerfreeze`, form the so-called three-equation-model
0422 (e.g., Hellmer and Olbers (1989) :cite:`hellmer:89`, Jenkins et al. (2001)
0423 :cite:`jenkins:01`). These equations are solved to obtain :math:`S_b, T_b,
0424 q`, which are then substituted into :eq:`jenkinsbc` to obtain the boundary
0425 conditions for the temperature and salinity equations of the ocean
0426 model. Note that with :math:`S_{I}=0` and :eq:`hellmersaltbalance`, the
0427 boundary flux for salinity becomes :math:`F_S = q\,S`, which is the flux
0428 that is necessary to account for the dilution of salinity in the case of
0429 melting.
0430
0431 The three-equation-model tends to yield smaller melt rates than the simpler
0432 formulation of the ISOMIP protocol (:numref:`shelfice_isomip`) because the
0433 freshwater flux (due to melting) decreases the salinity, which in turn
0434 raises the freezing point temperature and thus leads to less melting at the
0435 interface. For a simpler thermodynamics model where :math:`S_b` is not
0436 computed explicitly, for example as in the ISOMIP protocol, :eq:`jenkinsbc`
0437 cannot be applied directly. In this case :eq:`hellmersaltbalance` can be
0438 used with :eq:`jenkinsbc` to obtain:
0439
0440 .. math:: F_X = q\,(S-S_I)
274298f12b Jeff*0441
0442 This formulation can be used for all cases for which
67dae37801 Eli *0443 :eq:`hellmersaltbalance` is valid. Further, in this formulation it is
0444 obvious that melting (:math:`q<0`) leads to a reduction of salinity.
274298f12b Jeff*0445
67dae37801 Eli *0446 The default value of :varlink:`SHELFICEconserve` ``=.FALSE.`` removes the
0447 contribution :math:`q\, ( X_{b}-X )` from :eq:`jenkinsbc`, making the
0448 boundary conditions non-conservative.
274298f12b Jeff*0449
0a8794f5ee Mart*0450 Solving the three-equations system
0451 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0452
67dae37801 Eli *0453 There has been some confusion about the three-equations system, so we
0454 document the solution in the code here: We use :eq:`hellmerfreeze`
0455 :math:`T_{b} = a_{0} S_{b} + \epsilon_{4}` to eliminate :math:`T_{b}`
0456 from :eq:`jenkinsheatbudget` with :eq:`dTdzdiffus` and find an
0457 expression for the freshwater flux :math:`q`:
0a8794f5ee Mart*0458
0459 .. math::
0460 \begin{aligned}
0461 -Lq &= \epsilon_{1} (T - a_{0} S_{b} - \epsilon_{4})
0462 + \epsilon_{3} (T_{S} - a_{0} S_{b} - \epsilon_{4}) \\
0463 \Leftrightarrow Lq &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}
0464 + \epsilon_{q}
0465 \end{aligned}
0466 :label: solvedmeltrate
0467
0468 to be substituted into :eq:`hellmersaltbalance`:
0469
0470 .. math::
0471 \begin{aligned}
0472 \epsilon_{2}\,(S - S_{b}) &= - Lq\,(S_{b}-S_{I})
0473 = - (a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}
0474 + \epsilon_{q})\,(S_{b}-S_{I}) \\
0475 \Leftrightarrow 0 &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}^{2}
0476 + \{ \epsilon_{q} - \epsilon_{2}
0477 - a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{I} \}\,S_{b}
00f81e6785 Ou W*0478 + \epsilon_{2}\,S - \epsilon_{q}\,S_{I} \\
0479 \Leftrightarrow 0 &= A\,S_{b}^{2} + B\,S_{b} + C \\
0480 \Rightarrow S_{b} &= \frac{-B \pm \sqrt{ B^{2} - 4AC }}{2A}
0481 \end{aligned}
0482
0483 with the abbrevations
0484
0485 .. math::
0486 \begin{aligned}
0487 \epsilon_{1} &= c_{p} \, \rho_c \, \gamma_{T}, \quad
0488 \epsilon_{2} = \rho_c L \, \gamma_{S}, \quad
0489 \epsilon_{3} = \frac{\rho_{I} \, c_{p,I} \, \kappa}{h}, \quad
0490 \epsilon_{4} = b_{0}p + c_{0}, \\
0491 \epsilon_{q} &= \epsilon_{1}\,(\epsilon_{4} - T)
0492 + \epsilon_{3}\,(\epsilon_{4} - T_{S}), \\
0493 A &= a_{0}\,(\epsilon_{1} + \epsilon_{3}), \quad
0494 B = \epsilon_{q} - \epsilon_{2} - a_{0}\,(\epsilon_{1} +
0495 \epsilon_{3})\,S_{I}, \quad
0496 C = \epsilon_{2}\,S -\epsilon_{q}\,S_{I}.
0a8794f5ee Mart*0497 \end{aligned}
0498
00f81e6785 Ou W*0499 The smaller non-negative root of the quadratic equation in :math:`S_{b}` is
0500 used. By default, the ice shelf salinity :math:`S_{I}` is zero and the
0501 quadratic equation simplifies to
0a8794f5ee Mart*0502
0503 .. math::
0504 \begin{aligned}
0505 0 &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}^{2}
0506 + (\epsilon_{q} - \epsilon_{2}) \,S_{b} + \epsilon_{2}\,S \\
0507 S_{b} &= \frac{\epsilon_{2} - \epsilon_{q}\mp
0508 \sqrt{(\epsilon_{q} - \epsilon_{2})^2
b8f9311cc6 Matt*0509 - 4\, a_{0}\,(\epsilon_{1} + \epsilon_{3})\,\epsilon_{2}\,S}}
0a8794f5ee Mart*0510 {2\,a_{0}\,(\epsilon_{1} + \epsilon_{3})}
0511 \end{aligned}
0512
00f81e6785 Ou W*0513 With :math:`S_b`, the boundary layer temperature :math:`T_b` and the melt rate
0514 :math:`q` are known through :eq:`hellmerfreeze` and :eq:`solvedmeltrate`.
0a8794f5ee Mart*0515
67dae37801 Eli *0516 .. _shelfice_isomip:
0517
0a8794f5ee Mart*0518 ISOMIP thermodynamics
274298f12b Jeff*0519 ^^^^^^^^^^^^^^^^^^^^^
0520
0521 A simpler formulation follows the ISOMIP protocol. The
0522 freezing and melting in the boundary layer between ice shelf and ocean
0523 is parameterized following Grosfeld et al. (1997) :cite:`grosfeld:97`. In this
67dae37801 Eli *0524 formulation :eq:`jenkinsheatbudget` reduces to
274298f12b Jeff*0525
0526 .. math::
67dae37801 Eli *0527 c_{p} \, \rho_c \, \gamma_T (T - T_{b}) = -L q
274298f12b Jeff*0528 :label: isomipheatbalance
0529
0530 and the fresh water flux :math:`q` is computed from
0531
0532 .. math::
67dae37801 Eli *0533 q = - \frac{c_{p} \, \rho_c \, \gamma_T (T - T_{b})}{L}
274298f12b Jeff*0534 :label: isomipfwflx
0535
67dae37801 Eli *0536 In order to use this formulation, set runtime parameter
0537 :varlink:`useISOMIPTD` ``=.TRUE.`` in ``data.shelfice``.
0538
0539 .. _shelfice_exchange_coeff:
274298f12b Jeff*0540
0a8794f5ee Mart*0541 Exchange coefficients
0542 ^^^^^^^^^^^^^^^^^^^^^
0543
0544 The default exchange coefficents :math:`\gamma_{T/S}` are constant and
7b8b86ab99 Timo*0545 set by the run-time parameters :varlink:`SHELFICEheatTransCoeff` and
0a8794f5ee Mart*0546 :varlink:`SHELFICEsaltTransCoeff` (see
0547 :numref:`tab_phys_pkg_shelfice_runtimeparms`). If
0548 :varlink:`SHELFICEuseGammaFrict` ``=.TRUE.``, exchange coefficients
0549 are computed from drag laws and friction velocities estimated from
0550 ocean speeds following Holland and Jenkins (1999)
67dae37801 Eli *0551 :cite:`holland:99`. This computation can be modified using runtime
0552 parameters and the user is referred to
0a8794f5ee Mart*0553 :filelink:`pkg/shelfice/shelfice_readparms.F` for details.
0554
274298f12b Jeff*0555 Remark
0556 ^^^^^^
0557
0558 The shelfice package and experiments demonstrating its strengths and
67dae37801 Eli *0559 weaknesses are also described in Losch (2008)
0560 :cite:`losch:08`. Unfortunately, the description of the
0561 thermodynamics in the appendix of Losch (2008) is wrong.
274298f12b Jeff*0562
0563 .. _shelfice_subroutines:
0564
0565 Key subroutines
0566 ~~~~~~~~~~~~~~~
0567
0568 The main routine is :filelink:`shelfice_thermodynamics.F <pkg/shelfice/shelfice_thermodynamics.F>`
0569 but note that :filelink:`/pkg/shelfice` routines are also called when solving the momentum equations.
0570
0571 ::
0572
0573 C !CALLING SEQUENCE:
0574 C ...
0575 C |-FORWARD_STEP :: Step forward a time-step ( AT LAST !!! )
0576 C ...
0577 C | |-DO_OCEANIC_PHY :: Control oceanic physics and parameterization
0578 C ...
0579 C | | |-SHELFICE_THERMODYNAMICS :: main routine for thermodynamics
0580 C with diagnostics
0581 C ...
0582 C | |-THERMODYNAMICS :: theta, salt + tracer equations driver.
0583 C ...
0584 C | | |-EXTERNAL_FORCING_T :: Problem specific forcing for temperature.
0585 C | | |-SHELFICE_FORCING_T :: apply heat fluxes from ice shelf model
0586 C ...
0587 C | | |-EXTERNAL_FORCING_S :: Problem specific forcing for salinity.
0588 C | | |-SHELFICE_FORCING_S :: apply fresh water fluxes from ice shelf model
0589 C ...
0590 C | |-DYNAMICS :: Momentum equations driver.
0591 C ...
0592 C | | |-MOM_FLUXFORM :: Flux form mom eqn. package ( see
0593 C ...
0594 C | | | |-SHELFICE_U_DRAG :: apply drag along ice shelf to u-equation
0595 C with diagnostics
0596 C ...
0597 C | | |-MOM_VECINV :: Vector invariant form mom eqn. package ( see
0598 C ...
0599 C | | | |-SHELFICE_V_DRAG :: apply drag along ice shelf to v-equation
0600 C with diagnostics
0601 C ...
0602 C o
0603
0604
0605
0606 .. _shelfice_diagnostics:
0607
0608 SHELFICE diagnostics
0609 ~~~~~~~~~~~~~~~~~~~~
0610
0611 Diagnostics output is available via the diagnostics package (see
0612 :numref:`outp_pack`). Available output fields are summarized as follows:
0613
0614
0615 ::
0616
0617 ---------+----+----+----------------+-----------------
0618 <-Name->|Levs|grid|<-- Units -->|<- Tile (max=80c)
0619 ---------+----+----+----------------+-----------------
0620 SHIfwFlx| 1 |SM |kg/m^2/s |Ice shelf fresh water flux (positive upward)
0621 SHIhtFlx| 1 |SM |W/m^2 |Ice shelf heat flux (positive upward)
0622 SHIUDrag| 30 |UU |m/s^2 |U momentum tendency from ice shelf drag
0623 SHIVDrag| 30 |VV |m/s^2 |V momentum tendency from ice shelf drag
0624 SHIForcT| 1 |SM |W/m^2 |Ice shelf forcing for theta, >0 increases theta
0625 SHIForcS| 1 |SM |g/m^2/s |Ice shelf forcing for salt, >0 increases salt
0626
0627 Experiments and tutorials that use shelfice
0628 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0629
bd7056e9aa Jeff*0630 See the verification experiment :filelink:`isomip <verification/isomip>` for example usage of :filelink:`pkg/shelfice`.