8.6.2. SEAICE Package¶
Authors: Martin Losch, Dimitris Menemenlis, An Nguyen, JeanMichel Campin, Patrick Heimbach, Chris Hill, Jinlun Zhang, and Damien Ringeisen
8.6.2.1. Introduction¶
Package seaice provides a dynamic and thermodynamic interactive sea ice model.
CPP options enable or disable different aspects of the package
(Section 8.6.2.2). Runtime options, flags, filenames and
fieldrelated dates/times are set in data.seaice
(Section 8.6.2.3). A description of key subroutines is
given in Section 8.6.2.8. Available diagnostics
output is listed in Section 8.6.2.9.
8.6.2.2. SEAICE configuration and compiling¶
8.6.2.2.1. Compiletime options¶
As with all MITgcm packages, SEAICE can be turned on or off at compile time (see Section 3.5)
using the
packages.conf
file by addingseaice
to itor using genmake2 adding
enable=seaice
ordisable=seaice
switchesrequired packages and CPP options: seaice requires the external forcing package pkg/exf to be enabled; no additional CPP options are required.
Parts of the seaice code can be enabled or disabled at compile time via CPP preprocessor flags. These options are set in SEAICE_OPTIONS.h. Table 8.18 summarizes the most important ones. For more options see SEAICE_OPTIONS.h.
CPP option 
Default 
Description” 

#undef 
enhance STDOUT for debugging 

#define 
sea ice dynamics code 

#define 
LSR solver on Cgrid (rather than original Bgrid) 

#define 
enable use of EVP rheology solver 

#define 
enable use of JFNK rheology solver 

#define 
enable use of Krylov rheology solver 

#undef 
enable use of the truncated ellipse method (TEM) and coulombic yield curve 

#undef 
enable use of MohrCoulomb yield curve with shear flow rule 

#undef 
enable use of MohrCoulomb yield curve with elliptical plastic potential 

#undef 
enable use of teardrop and parabolic Lens yield curves with normal flow rules 

#undef 
use a coloring method for LSR solver 

#undef 
enable solve approximate sea ice momentum equation and bypass solving for sea ice internal stress 

#define 
use pkg/exfcomputed fluxes as starting point 

#define 
use differentiable regularization for viscosities 

#undef 
use differentiable regularization for \(1/\Delta\) 

#undef 
enable grounding parameterization for improved fastice in shallow seas 

#undef 
run with dynamical sea Ice Thickness Distribution (ITD) 

#undef 
enable sea ice with variable salinity 

#undef 
enable to limit seaice load (siceLoad) on the sea surface 

#undef 
enable sea ice tracer package 

#undef 
Bgrid only for backward compatiblity: turn on icestress on ocean 

#undef 
Bgrid only for backward compatiblity: use ETAN for tilt computations rather than geostrophic velocities 

#undef 
use of adjointable but more simplified sea ice thermodynamics model in seaice_growth_adx.F instead of seaice_growth.F 
8.6.2.3. Runtime parameters¶
Runtime parameters (see Table 8.19) are set in
data.seaice
(read in pkg/seaice/seaice_readparms.F).
8.6.2.3.1. Enabling the package¶
seaice package is switched on/off at runtime by
setting useSEAICE = .TRUE.,
in data.pkg
.
8.6.2.3.2. General flags and parameters¶
Table 8.19 lists most runtime parameters.
Name 
Default value 
Description 

FALSE 
write sea ice state to file 

TRUE 
use dynamics 

FALSE 
use the JFNKsolver 

FALSE 
use truncated ellipse method or coulombic yield curve 

FALSE 
use the MohrCoulomb yield curve with shear flow rule 

FALSE 
use the MohrCoulomb yield curve with elliptical plastic potential 

FALSE 
use the teardrop yield curve with normal flow rule 

FALSE 
use the parabolic Lens yield curve with normal flow rule 

FALSE 
use strength implicit coupling in LSR/JFNK 

TRUE 
use metric terms in dynamics 

TRUE 
use EVP pickups 

FALSE 
solve approximate momentum equation, bypassing rheology 

TRUE 
use flux form for 2nd central difference advection scheme 

FALSE 
enable restoring to climatology under ice 

TRUE 
update ocean surface stress accounting for sea ice cover 

TRUE 
scale atmosphere and oceansurface stress on ice by concentration (AREA) 

TRUE 
in computing seaiceMass, add snow contribution 

FALSE 
turn on iceocean stress coupling following 

TRUE 
flag to turn off zerolayerthermodynamics for testing 

TRUE 
flag to turn off advection of scalar variable HEFF 

TRUE 
flag to turn off advection of scalar variable AREA 

TRUE 
flag to turn off advection of scalar variable HSNOW 

TRUE 
flag to turn off advection of scalar variable HSALT 

77 
set advection scheme for seaice scalar state variables 

TRUE 
use floodfreeze algorithm 

1.0 
over/undershoot factor for seaice advective term in forward/adjoint (SEAICE_USE_GROWTH_ADX only) 

FALSE 
use noslip boundary conditions instead of freeslip 

dTtracerLev (1) 
time step for seaice thermodynamics (s) 

dTtracerLev (1) 
time step for seaice dynamics (s) 

0.0 
EVP subcycling time step (s); values \(>\) 0 turn on EVP 

TRUE 
use modified EVP* instead of EVP, following [LKT+12] 

TRUE 
“revisited form” variation on EVP*, following [BFLM13] 

unset 
number of modified EVP* iterations 

unset 
EVP* parameter (nondim.), to replace 2*SEAICE_evpTauRelax/SEAICE_deltaTevp 

unset 
EVP* parameter (nondim.), to replace SEAICE_deltaTdyn/SEAICE_deltaTevp 

unset 
largest stabilized frequency for adaptive EVP (nondim.) 

4.0 
aEVP multiple of stability factor (nondim.), see [KDL16] \(\alpha * \beta = c^\ast * \gamma\) 

5.0 
aEVP lower limit of alpha and beta (nondim.), see [KDL16] 

0.33333333 
EVP parameter \(E_0\) (nondim.), sets relaxation timescale SEAICE_evpTauRelax = SEAICE_elasticParm * SEAICE_deltaTdyn 

relaxation time scale \(T\) for EVP waves (s) 

OLx  2 
overlap for LSRsolver or preconditioner, \(x\)dimension 

OLy  2 
overlap for LSRsolver or preconditioner, \(y\)dimension 

2/10 
maximum number of nonlinear (outer loop) iterations (LSR/JFNK) 

1500/10 
maximum number of linear iterations (LSR/JFNK) 

(off) 
start line search after “lsIter” Newton iterations 

4 
maximum number of line search steps 

0.5 
line search step size parameter 

1.0E05 
nonlinear tolerance parameter for JFNK solver 

0.10 
minimum tolerance parameter for linear JFNK solver 

0.99 
maximum tolerance parameter for linear JFNK solver 

unset 
tolerance parameter for FGMRES residual 

1.0E06 
step size for the FDgradient in s/r seaice_jacvec 

dumpFreq 
dump frequency (s) 

TRUE 
write snapshot using /pkg/mdsio 

FALSE 
write snapshot using /pkg/mnc 

0.0 
initial sea ice thickness averaged over grid cell (m) 

1.0E03 
airice drag coefficient (nondim.) 

1.0E03 
airocean drag coefficient (nondim.) 

5.5E03 
waterice drag coefficient (nondim.) 

0.75 
winter sea ice albedo 

0.66 
summer sea ice albedo 

0.84 
dry snow albedo 

0.70 
wet snow albedo 

0.10 
water albedo (not used if #define SEAICE_EXTERNAL_FLUXES) 

2.75E+04 
sea ice strength constant \(P^{\ast}\) (N/m^{2}) 

20.0 
sea ice strength constant \(C^{\ast}\) (nondim.) 

2.0 
VP rheology ellipse aspect ratio \(e\) 

sea ice plastic potential ellipse aspect ratio \(e_G\) 

1.0 
slope of the MohrCoulomb yield curve 

1.0 
use replacement pressure (0.01.0) 

0.0 
tensile factor for the yield curve 

1.3 (or pkg/exf value) 
density of air (kg/m^{3}) 

1004.0 (or pkg/exf value) 
specific heat of air (J/kg/K) 

2.5E+06 (or pkg/exf value) 
latent heat of evaporation (J/kg) 

3.34E+05 (or pkg/exf value) 
latent heat of fusion (J/kg) 

1.75E03 
iceocean transfer coefficient for latent and sensible heat (nondim.) 

FALSE 
use Maykut polynomial to compute saturation vapor pressure 

2.16560E+00 
sea ice conductivity (W m^{1} K^{1}) 

3.10000E01 
snow conductivity (W m^{1} K^{1}) 

0.970018 (or pkg/exf value) 
longwave ocean surface emissivity (nondim.) 

0.15 
cutoff snow thickness to use snow albedo (m) 

0.30 
ice penetration shortwave radiation factor (nondim.) 

0.0 
salinity newly formed ice (as fraction of ocean surface salinity) 

1.0 (or computed from other parms) 
frazil to sea ice conversion rate, as fraction (relative to the local freezing point of sea ice water) 

1.0 
scaling factor for ice area in computing total ocean stress (nondim.) 

unset 
filename for initial sea ice eff. thickness field HEFF (m) 

unset 
filename for initial fraction sea ice cover AREA (nondim.) 

unset 
filename for initial eff. snow thickness field HSNOW (m) 

unset 
filename for initial eff. sea ice salinity field HSALT (g/m^{2}) 

1.0E05 
sets accuracy of LSR solver 

0.0 
parameter used in advect.F 

0.5 
lead closing parameter \(h_0\) (m); demarcation thickness between thick and thin ice which determines partition between vertical and lateral ice growth 

50.0 
minimum air temperature (^{o}C) 

60.0 
minimum downward longwave (W/m^{2}) 

50.0 
minimum ice temperature (^{o}C) 

10 
number of iterations for ice surface temperature solution 

1.0E10 
a “small number” used in various routines 

1.0E5 
minimum concentration to regularize ice thickness 

0.05 
minimum ice thickness (m) for regularization 

1 
number of ice categories for thermodynamics 

TRUE 
use same fixed pdf for snow as for multithicknesscategory ice 
The following dynamical ice thickness distribution and ridging parameters in Table 8.20 are only active with #define SEAICE_ITD. All parameters are nondimensional unless indicated.
Name 
Default value 
Description 

TRUE 
use [Hib79] ice strength; do not use [Rot75] with #define SEAICE_ITD 

TRUE 
use simple ridging a la [Hib79] 

17.0 
scaling parameter of [Rot75] ice strength parameterization 

0 
use partition function of [TRMC75] 

0 
use redistribution function of [Hib80] 

10 
maximum number of ridging sweeps 

0.5 
fraction of shear to be used for ridging 

0.15 
max. ice conc. that participates in ridging [TRMC75] 

25.0 

0.05 
similar to SEAICEgStar for [LHMJ07] participation function 

3.0 
similar to SEAICEhStar for [LHMJ07] ridging function 

1.0 
regularization parameter for rafting 

0.5 
fraction of snow that remains on ridged ice 

TRUE 
use linear remapping scheme of [Lip01] 

unset 
nITD+1array of ice thickness category limits (m) 

3.0, 15.0, 3.0 
when Hlimit is not set, then these parameters determine Hlimit from a simple function following [Lip01] 
8.6.2.4. Description¶
The MITgcm sea ice model is based on a variant of the viscousplastic (VP) dynamicthermodynamic sea ice model (Zhang and Hibler 1997 [ZH97]) first introduced in Hibler (1979) and Hibler (1980) [Hib79, Hib80]. In order to adapt this model to the requirements of coupled iceocean state estimation, many important aspects of the original code have been modified and improved, see Losch et al. (2010) [LMC+10]:
the code has been rewritten for an Arakawa Cgrid, both B and Cgrid variants are available; the Cgrid code allows for noslip and freeslip lateral boundary conditions;
three different solution methods for solving the nonlinear momentum equations have been adopted: LSOR (Zhang and Hibler 1997 [ZH97]), EVP (Hunke and Dukowicz 1997 [HD97]), JFNK (Lemieux et al. 2010 [LTSedlavcek+10], Losch et al. 2014 [LFLV14]);
iceocean stress can be formulated as in Hibler and Bryan (1987) [HB87] or as in Campin et al. (2008) [CMF08];
ice variables are advected by sophisticated, conservative advection schemes with flux limiting;
growth and melt parameterizations have been refined and extended in order to allow for more stable automatic differentiation of the code.
The sea ice model is tightly coupled to the ocean compontent of the MITgcm. Heat, fresh water fluxes and surface stresses are computed from the atmospheric state and, by default, modified by the ice model at every time step.
The ice dynamics models that are most widely used for largescale climate studies are the viscousplastic (VP) model (Hilber 1979 [Hib79]), the cavitating fluid (CF) model (Flato and Hibler 1992 [FWDH92]), and the elasticviscousplastic (EVP) model (Hunke and Dukowicz 1997 [HD97]). Compared to the VP model, the CF model does not allow ice shear in calculating ice motion, stress, and deformation. EVP models approximate VP by adding an elastic term to the equations for easier adaptation to parallel computers. Because of its higher accuracy in plastic solution and relatively simpler formulation, compared to the EVP model, we decided to use the VP model as the default dynamic component of our ice model. To do this we extended the line successive over relaxation (LSOR) method of Zhang and Hibler (1997) [ZH97] for use in a parallel configuration. An EVP model and a freedrift implementation can be selected with runtime flags.
pkg/seaice includes the original socalled zerolayer thermodynamics with a snow cover as in the appendix of Semtner (1976) [Sem76]. Two versions of this zerolayer thermodynamic code exist, with a more developed version seaice_growth.F and a simplified version seaice_growth_adx.F based on Fenty (2013) [FH13] that excludes physics such as ITD, treatment for sublimation, and frazil ice but provides a stable sea ice adjointable with physical sensitivity. When the seaice_growth_adx code is enabled (by defining SEAICE_USE_GROWTH_ADX in SEAICE_OPTIONS.h), the regularization parameter SINegFac is set to zero in adjoint mode to disable the potential propagation of unphysical terms associated with sea ice dynamics.
8.6.2.4.1. Compatibility with icethermodynamics package pkg/thsice¶
The zerolayer thermodynamic model assumes that ice does not store heat and, therefore, tends to exaggerate the seasonal variability in ice thickness. This exaggeration can be significantly reduced by using Winton’s (Winton 2000 [Win00]) threelayer thermodynamic model that permits heat storage in ice.
The Winton (2000) seaice thermodynamics have been ported to MITgcm; they
currently reside under pkg/thsice, described in
Section 8.6.1. It is fully compatible with the packages
seaice and exf. When turned on
together with seaice, the zerolayer thermodynamics
are replaced by the Winton thermodynamics. In order to use package
seaice with the thermodynamics of
pkg/thsice, compile both packages and turn both package on in
data.pkg
; see an example in
verification/global_ocean.cs32x15/input.icedyn. Note, that once
thsice is turned on, the variables and diagnostics
associated to the default thermodynamics are meaningless, and the diagnostics
of thsice must be used instead.
8.6.2.4.2. Surface forcing¶
The sea ice model requires the following input fields: 10 m winds, 2 m air temperature and specific humidity, downward longwave and shortwave radiations, precipitation, evaporation, and river and glacier runoff. The sea ice model also requires surface temperature from the ocean model and the top level horizontal velocity. Output fields are surface wind stress, evaporation minus precipitation minus runoff, net surface heat flux, and net shortwave flux. The seaice model is global: in icefree regions bulk formulae (by default computed in package exf) are used to estimate oceanic forcing from the atmospheric fields.
8.6.2.5. Dynamics¶
The momentum equation of the seaice model is
where \(m=m_{i}+m_{s}\) is the ice and snow mass per unit area; \(\mathbf{u}=u\hat{\mathbf{i}}+v\hat{\mathbf{j}}\) is the ice velocity vector; \(\hat{\mathbf{i}}\), \(\hat{\mathbf{j}}\), and \(\hat{\mathbf{k}}\) are unit vectors in the \(x\), \(y\), and \(z\) directions, respectively; \(f\) is the Coriolis parameter; \(\mathbf{\tau}_\mathrm{air}\) and \(\mathbf{\tau}_\mathrm{ocean}\) are the windice and oceanice stresses, respectively; \(g\) is the gravity accelation; \(\nabla\phi(0)\) is the gradient (or tilt) of the sea surface height; \(\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}\) is the sea surface height potential in response to ocean dynamics (\(g\eta\)), to atmospheric pressure loading (\(p_{a}/\rho_{0}\), where \(\rho_{0}\) is a reference density) and a term due to snow and ice loading ; and \(\mathbf{F}= \nabla \cdot\sigma\) is the divergence of the internal ice stress tensor \(\sigma_{ij}\). Advection of seaice momentum is neglected. The wind and iceocean stress terms are given by
where \(\mathbf{U}_\mathrm{air/ocean}\) are the surface winds of the atmosphere and surface currents of the ocean, respectively; \(C_\mathrm{air/ocean}\) are air and ocean drag coefficients; \(\rho_\mathrm{air/ocean}\) are reference densities; and \(R_\mathrm{air/ocean}\) are rotation matrices that act on the wind/current vectors.
8.6.2.5.1. ViscousPlastic (VP) Rheology¶
For an isotropic system the stress tensor \(\sigma_{ij}\) (\(i,j=1,2\)) can be related to the ice strain rate and strength by a nonlinear viscousplastic (VP) constitutive law:
The ice strain rate is given by
The maximum ice pressure \(P_{\max}\) (variable PRESS0 in the code), a measure of ice strength, depends on both thickness \(h\) and compactness (concentration) \(c\):
with the constants \(P^{\ast}\) (runtime parameter SEAICE_strength) and \(C^{\ast}\) (runtime parameter SEAICE_cStar). By default, \(P\) (variable PRESS in the code) is the replacement pressure
(8.6)¶\[P = (1k_t)\,P_{\max} \left( (1  f_{r}) + f_{r} \frac{\Delta}{\Delta_{\rm reg}} \right)\]
where \(f_{r}\) is runtime parameter SEAICEpressReplFac (default = 1.0), and \(\Delta_{\rm reg}\) is a regularized form of \(\Delta = \left[ \left(\dot{\epsilon}_{11}+\dot{\epsilon}_{22}\right)^2 + e^{2}\left( \left(\dot{\epsilon}_{11}\dot{\epsilon}_{22} \right)^2 + \dot{\epsilon}_{12}^2 \right) \right]^{\frac{1}{2}}\), for example \(\Delta_{\rm reg} = \max(\Delta,\Delta_{\min})\).
The tensile strength factor \(k_t\) (runtime parameter SEAICE_tensilFac) determines the ice tensile strength \(T = k_t\cdot P_{\max}\), as defined by König Beatty and Holland (2010) [KBH10]. SEAICE_tensilFac is zero by default.
Different VP rheologies can be used to model sea ice dynamics. The different rheologies are characterized by different definitions of the bulk and shear viscosities \(\zeta\) and \(\eta\) in (8.4). The following Table 8.21 is a summary of the available choices with recommended (sensible) parameter values. All the rheologies presented here depend on the ice strength \(P\) (8.6).
Name 
CPP flags 
Runtime flags (recommended value) 

None (default) 


None 











Note: With the exception of the default rheology and the TEM (with SEAICEmcMU : \(\mu=1.0\)), these rheologies are not implemented in EVP (Section 8.6.2.5.3).
Elliptical yield curve with normal flow rule¶
The default rheology in the sea ice module of the MITgcm implements the widely used elliptical yield curve with a normal flow rule [Hib79]. For this yield curve, the nonlinear bulk and shear viscosities \(\zeta\) and \(\eta\) are functions of ice strain rate invariants and ice strength such that the principal components of the stress lie on an elliptical yield curve with the ratio of major to minor axis \(e = 2.0\) (runtime parameter SEAICE_eccen); they are given by:
with the abbreviation
\[\Delta = \left[ \left(\dot{\epsilon}_{11}+\dot{\epsilon}_{22}\right)^2 + e^{2}\left( \left(\dot{\epsilon}_{11}\dot{\epsilon}_{22} \right)^2 + \dot{\epsilon}_{12}^2 \right) \right]^{\frac{1}{2}}\]
The bulk viscosities are bounded above by imposing both a minimum \(\Delta_{\min}\) (for numerical reasons, runtime parameter SEAICE_deltaMin is set to a default value of \(10^{10}\,\text{s}^{1}\), the value of SEAICE_EPS) and a maximum \(\zeta_{\max} = P_{\max}/(2\Delta^\ast)\), where \(\Delta^\ast=(2\times10^4/5\times10^{12})\,\text{s}^{1}\) \(= 2\times10^{9}\,\text{s}^{1}\). Obviously, this corresponds to regularizing \(\Delta\) with the typical value of SEAICE_deltaMin \(= 2\times10^{9}\). Clearly, some of this regularization is redundant. (There is also the option of bounding \(\zeta\) from below by setting runtime parameter SEAICE_zetaMin \(>0\), but this is generally not recommended). For stress tensor computation the replacement pressure \(P = 2\,\Delta\zeta\) is used so that the stress state always lies on the elliptic yield curve by definition.
Defining the CPPflag SEAICE_ZETA_SMOOTHREG in SEAICE_OPTIONS.h before compiling replaces the method for bounding \(\zeta\) by a smooth (differentiable) expression:
where \(\Delta_{\min}=10^{20}\,\text{s}^{1}\) should be chosen to avoid divisions by zero.
In this default formulation the yield curve does not allow isotropic tensile stress, that is, sea ice can be “pulled apart” without any effort. Setting the parameter \(k_t\) (SEAICE_tensilFac) to a small value larger than zero, extends the yield curve into a region where the divergence of the stress \(\sigma_{11}+\sigma_{22} > 0\) to allow some tensile stress.
Besides this commonly used default rheology, a number of a alternative rheologies are implemented. Some of these are experiemental and should be used with caution.
Elliptical yield curve with nonnormal flow rule¶
Defining the runtime parameter SEAICE_eccfr with a value different from SEAICE_eccen allows one to use an elliptical yield curve with a nonnormal flow rule as described in Ringeisen et al. (2020) [RTL20]. In this case the viscosities are functions of \(e_F\) (SEAICE_eccen) and \(e_G\) (SEAICE_eccfr):
with the abbreviation
Note that if \(e_G=e_F=e\), these formulae reduce to the normal flow rule.
Truncated ellipse method (TEM) for elliptical yield curve¶
In the socalled truncated ellipse method, the shear viscosity \(\eta\) is capped to suppress any tensile stress:
To enable this method, set #define
SEAICE_ALLOW_TEM in
SEAICE_OPTIONS.h and turn it on with
SEAICEuseTEM =.TRUE.,
in data.seaice
. This parameter
combination implies the default of SEAICEmcMU \(= 1.0\).
Instead of an ellipse that is truncated by constant slope coulombic limbs, this
yield curve can also be seen as a MohrCoulomb yield curve with elliptical flow
rule that is truncated for high \(P\) by an ellipse. As a consequence, the
MohrCoulomb slope SEAICEmcMU can be set in data.seaice
to
values \(\ne 1.0\). This defines a coulombic yield curve similar to the
ones shown in Hibler and Schulson (2000) [HS00] and Ringeisen et
al. (2019) [RLTNH19].
For this rheology, it is recommended to use a nonzero tensile strength, so set
SEAICE_tensilFac \(=k_{t}>0\) in data.seaice
, e.g., \(=
0.05\) or 5%.
MohrCoulomb yield curve with elliptical plastic potential¶
To use a MohrCoulomb rheology, set #define
SEAICE_ALLOW_MCE in
SEAICE_OPTIONS.h and
SEAICEuseMCE = .TRUE.,
in data.seaice
. This MohrCoulomb
yield curve uses an elliptical plastic potential to define the flow rule. The
slope of the MohrCoulomb yield curve is defined by SEAICEmcMU in
data.seaice
, and the plastic potential ellipse aspect ratio is set by
SEAICE_eccfr in data.seaice
. For details of this rheology, see
https://doi.org/10.26092/elib/380, Chapter 2.
For this rheology, it is recommended to use a nonzero tensile strength, so set
SEAICE_tensilFac \(>0\) in data.seaice
, e.g., \(= 0.05\)
or 5%.
MohrCoulomb yield curve with shear flow rule¶
To use the specifc MohrCoulomb rheology as defined first by Ip et al. (1991)
[IHF91], set #define
SEAICE_ALLOW_MCS in
SEAICE_OPTIONS.h and
SEAICEuseMCS = .TRUE.,
in data.seaice
. The slope of the
MohrCoulomb yield curve is defined by SEAICEmcMU in
data.seaice
. For details of this rheology, including the tensile strength,
see https://doi.org/10.26092/elib/380, Chapter 2.
For this rheology, it is recommended to use a nonzero tensile strength, so set
SEAICE_tensilFac \(>0\) in data.seaice
, e.g., \(= 0.05\)
or 5%.
WARNING: This rheology is known to be unstable. Use with caution!
Teardrop yield curve with normal flow rule¶
The teardrop rheology was first described in Zhang and Rothrock (2005) [ZR05]. Here we implement a slightly modified version (See https://doi.org/10.26092/elib/380, Chapter 2).
To use this rheology, set #define
SEAICE_ALLOW_TEARDROP in
SEAICE_OPTIONS.h and
SEAICEuseTD = .TRUE.,
in data.seaice
. The size of the yield
curve can be modified by changing the tensile strength, using
SEAICE_tensFac in data.seaice
.
For this rheology, it is recommended to use a nonzero tensile strength, so set
SEAICE_tensilFac \(>0\) in data.seaice
, e.g., \(=
0.025\) or 2.5%.
Parabolic lens yield curve with normal flow rule¶
The parabolic lens rheology was first described in Zhang and Rothrock (2005) [ZR05]. Here we implement a slightly modified version (See https://doi.org/10.26092/elib/380, Chapter 2).
To use this rheology, set #define
SEAICE_ALLOW_TEARDROP in
SEAICE_OPTIONS.h and
SEAICEusePL = .TRUE.,
in data.seaice
. The size of the yield
curve can be modified by changing the tensile strength, using
SEAICE_tensFac in data.seaice
.
For this rheology, it is recommended to use a nonzero tensile strength, so set
SEAICE_tensilFac \(>0\) in data.seaice
, e.g., \(=
0.025\) or 2.5%.
8.6.2.5.2. LSR and JFNK solver¶
In matrix notation, the discretized momentum equations can be written as
The solution vector \(\mathbf{x}\) consists of the two velocity components \(u\) and \(v\) that contain the velocity variables at all grid points and at one time level. The standard (and default) method for solving Eq. (8.10) in the sea ice component of MITgcm is an iterative Picard solver: in the \(k\)th iteration a linearized form \(\mathbf{A}(\mathbf{x}^{k1})\,\mathbf{x}^{k} = \mathbf{b}(\mathbf{x}^{k1})\) is solved (in the case of MITgcm it is a Line Successive (over) Relaxation (LSR) algorithm). Picard solvers converge slowly, but in practice the iteration is generally terminated after only a few nonlinear steps and the calculation continues with the next time level. This method is the default method in MITgcm. The number of nonlinear iteration steps or pseudotime steps can be controlled by the runtime parameter SEAICEnonLinIterMax. This parameter’s default is 2, but using a number of at least 10 is recommended for better solutions that are converged at least in an energy norm sense (Zhang and Hibler 1997) [ZH97].
In order to overcome the poor convergence of the Picard solver, Lemieux et al. (2010) [LTSedlavcek+10] introduced a Jacobianfree NewtonKrylov solver for the sea ice momentum equations. This solver is also implemented in MITgcm (see Losch et al. 2014 [LFLV14]). The Newton method transforms minimizing the residual \(\mathbf{F}(\mathbf{x}) = \mathbf{A}(\mathbf{x})\,\mathbf{x}  \mathbf{b}(\mathbf{x})\) to finding the roots of a multivariate Taylor expansion of the residual \(\mathbf{F}\) around the previous (\(k1\)) estimate \(\mathbf{x}^{k1}\):
with the Jacobian \(\mathbf{J}\equiv\mathbf{F}'\). The root \(\mathbf{F}(\mathbf{x}^{k1}+\delta\mathbf{x}^{k})=0\) is found by solving
for \(\delta\mathbf{x}^{k}\). The next (\(k\)th) estimate is given by \(\mathbf{x}^{k}=\mathbf{x}^{k1}+(1\gamma_{\mathrm{LS}})^{l} \,\delta\mathbf{x}^{k}\).
By default \(l=0\), but in order to avoid overshoots, the step size factor
\((1\gamma_{\mathrm{LS}})^{l}\) with \(\gamma_{\mathrm{LS}}<1\) can be
iteratively reduced in a line search with \(l=0,1,2,\ldots\) until
\(\\mathbf{F}(\mathbf{x}^k)\ < \\mathbf{F}(\mathbf{x}^{k1})\\), where
\(\\cdot\=\int\cdot\,dx^2\) is the \(L_2\)norm. The line search
starts after SEAICE_JFNK_lsIter nonlinear Newton iterations (off by
default) to allow for full Newton steps at the beginning of the iteration. If
the line search is turned on by setting SEAICE_JFNK_lsIter to a
nonnegative value in data.seaice
, by default, the line search with
\(\gamma_\mathrm{LS}=\frac{1}{2}\) (runtime parameter
SEAICE_JFNK_lsGamma) is stopped after \(L_{\max}=4\) (runtime
parameter SEAICE_JFNK_lsLmax) steps.
Forming the Jacobian \(\mathbf{J}\) explicitly is often avoided as “too error prone and time consuming”. Instead, Krylov methods only require the action of \(\mathbf{J}\) on an arbitrary vector \(\mathbf{w}\) and hence allow a matrix free algorithm for solving (8.12). The action of \(\mathbf{J}\) can be approximated by a firstorder Taylor series expansion:
or computed exactly with the help of automatic differentiation (AD) tools. SEAICE_JFNKepsilon sets the step size \(\epsilon\).
We use the Flexible Generalized Minimum RESidual (FMGRES) method with righthand side preconditioning to solve (8.12) iteratively starting from a first guess of \(\delta\mathbf{x}^{k}_{0} = 0\). For the preconditioning matrix \(\mathbf{P}\) we choose a simplified form of the system matrix \(\mathbf{A}(\mathbf{x}^{k1})\) where \(\mathbf{x}^{k1}\) is the estimate of the previous Newton step \(k1\). The transformed equation (8.12) becomes
The Krylov method iteratively improves the approximate solution to (8.14) in subspace (\(\mathbf{r}_0\), \(\mathbf{J}\mathbf{P}^{1}\mathbf{r}_0\), \((\mathbf{J}\mathbf{P}^{1})^2\mathbf{r}_0\), \(\dots\), \((\mathbf{J}\mathbf{P}^{1})^m\mathbf{r}_0\)) with increasing \(m\); \(\mathbf{r}_0 = \mathbf{F}(\mathbf{x}^{k1}) \mathbf{J}(\mathbf{x}^{k1})\,\delta\mathbf{x}^{k}_{0}\) is the initial residual of (8.12); \(\mathbf{r}_0=\mathbf{F}(\mathbf{x}^{k1})\) with the first guess \(\delta\mathbf{x}^{k}_{0}=0\). We allow a Krylov subspace of dimension \(m=50\) and we do allow restarts for more than 50 Krylov iterations. The preconditioning operation involves applying \(\mathbf{P}^{1}\) to the basis vectors \(\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_m\) of the Krylov subspace. This operation is approximated by solving the linear system \(\mathbf{P}\,\mathbf{w}=\mathbf{v}_i\). Because \(\mathbf{P} \approx \mathbf{A}(\mathbf{x}^{k1})\), we can use the LSR algorithm already implemented in the Picard solver. Each preconditioning operation uses a fixed number of 10 LSR iterations avoiding any termination criterion. More details and results can be found in Losch et al. (2014) [LFLV14]).
To use the JFNK solver set SEAICEuseJFNK = .TRUE.,
in the
namelist file data.seaice
; #define
SEAICE_ALLOW_JFNK in
SEAICE_OPTIONS.h and we recommend
using a smooth regularization of \(\zeta\) by #define
SEAICE_ZETA_SMOOTHREG (see above) for better convergence. The
nonlinear Newton iteration is terminated when the \(L_2\)norm of the
residual is reduced by \(\gamma_{\mathrm{nl}}\) (runtime parameter
SEAICEnonLinTol = 1.E4,
will already lead to expensive
simulations) with respect to the initial norm:
\(\\mathbf{F}(\mathbf{x}^k)\ <
\gamma_{\mathrm{nl}}\\mathbf{F}(\mathbf{x}^0)\\). Within a nonlinear
iteration, the linear FGMRES solver is terminated when the residual is smaller
than \(\gamma_k\\mathbf{F}(\mathbf{x}^{k1})\\) where \(\gamma_k\) is
determined by
so that the linear tolerance parameter \(\gamma_k\) decreases with the nonlinear Newton step as the nonlinear solution is approached. This inexact Newton method is generally more robust and computationally more efficient than exact methods. Typical parameter choices are \(\gamma_0 =\) JFNKgamma_lin_max \(= 0.99\), \(\gamma_{\min} =\) JFNKgamma_lin_min \(= 0.1\), and \(r =\) JFNKres_tFac \(\times\\mathbf{F}(\mathbf{x}^{0})\\) with JFNKres_tFac \(= 0.5\). We recommend a maximum number of nonlinear iterations SEAICEnewtonIterMax \(= 100\) and a maximum number of Krylov iterations SEAICEkrylovIterMax \(= 50\), because the Krylov subspace has a fixed dimension of 50 (but restarts are allowed for SEAICEkrylovIterMax \(> 50\)).
Setting SEAICEuseStrImpCpl to .TRUE.
turns on “strength implicit
coupling” (see Hutchings et al. 2004 [HJL04]) in the LSR solver
and in the LSR preconditioner for the JFNK solver. In this mode, the different
contributions of the stress divergence terms are reordered so as to increase
the diagonal dominance of the system matrix. Unfortunately, the convergence
rate of the LSR solver is increased only slightly, while the JFNK convergence
appears to be unaffected.
8.6.2.5.3. ElasticViscousPlastic (EVP) Dynamics¶
Hunke and Dukowicz (1997) [HD97] introduced an elastic contribution to the strain rate in order to regularize (8.4) in such a way that the resulting elasticviscousplastic (EVP) and VP models are identical at steady state,
The EVP model uses an explicit time stepping scheme with a short timestep. According to the recommendation in Hunke and Dukowicz (1997) [HD97], the EVPmodel should be stepped forward in time 120 times (SEAICE_deltaTevp = SEAICE_deltaTdyn /120) within the physical ocean model time step (although this parameter is under debate), to allow for elastic waves to disappear. Because the scheme does not require a matrix inversion it is fast in spite of the small internal timestep and simple to implement on parallel computers. For completeness, we repeat the equations for the components of the stress tensor \(\sigma_{1} = \sigma_{11}+\sigma_{22}\), \(\sigma_{2}= \sigma_{11}\sigma_{22}\), and \(\sigma_{12}\). Introducing the divergence \(D_D = \dot{\epsilon}_{11}+\dot{\epsilon}_{22}\), and the horizontal tension and shearing strain rates, \(D_T = \dot{\epsilon}_{11}\dot{\epsilon}_{22}\) and \(D_S = 2\dot{\epsilon}_{12}\), respectively, and using the above abbreviations, the equations (8.16) can be written as:
Here, the elastic parameter \(E\) is redefined in terms of a damping timescale \(T\) for elastic waves
\(T=E_{0}\Delta{t}\) with the tunable parameter \(E_0<1\) and the external (long) timestep \(\Delta{t}\). \(E_{0} = \frac{1}{3}\) is the default value in the code and close to what Hunke and Dukowicz (1997) [HD97] recommend.
We do not recommend to use the EVP solver in its original form. Instead, use
mEVP or aEVP instead (see Section 8.6.2.5.4). If you
really need to use the original EVP solver, make sure that both #define
SEAICE_CGRID and #define
SEAICE_ALLOW_EVP are set in
SEAICE_OPTIONS.h (both are defined by
default). By default, the runtime parameters SEAICEuseEVPstar and
SEAICEuseEVPrev are set to .TRUE.
, which already improves the
behavoir of EVP, but for the original EVP they should be set to .FALSE.
. The
solver is turned on by setting the subcycling time step
SEAICE_deltaTevp to a value larger than zero. The choice of this
time step is under debate. Hunke and Dukowicz (1997) [HD97] recommend
order 120 time steps for the EVP solver within one model time step
\(\Delta{t}\) (deltaTmom). One can also choose order 120 time
steps within the forcing time scale, but then we recommend adjusting the
damping time scale \(T\) accordingly, by setting either
SEAICE_elasticParm (\(E_{0}\)), so that \(E_{0}\Delta{t}=\)
forcing time scale, or directly SEAICE_evpTauRelax (\(T\)) to
the forcing time scale. (NOTE: with the improved EVP variants of the next
section, the above recommendations are obsolete. Use mEVP or aEVP instead.)
8.6.2.5.4. More stable variants of ElasticViscousPlastic Dynamics: EVP*, mEVP, and aEVP¶
The genuine EVP scheme appears to give noisy solutions (see Hunke 2001, Lemieux et al. 2012, Bouillon et a1. 2013 [BFLM13, Hun01, LKT+12]). This has led to a modified EVP or EVP* (Lemieux et al. 2012, Bouillon et a1. 2013, Kimmritz et al. 2015 [BFLM13, KDL15, LKT+12]); here, we refer to these variants by modified EVP (mEVP) and adaptive EVP (aEVP). The main idea is to modify the “natural” timediscretization of the momentum equations:
where \(n\) is the previous time step index, and \(p\) is the previous
subcycling index. The extra “intertial” term
\(m\,(\mathbf{u}^{p+1}\mathbf{u}^{n})/\Delta{t})\) allows the definition
of a residual \(\mathbf{u}^{p+1}\mathbf{u}^{p}\) that, as
\(\mathbf{u}^{p+1} \rightarrow \mathbf{u}^{n+1}\), converges to
\(0\). In this way EVP can be reinterpreted as a pure iterative solver
where the subcycling has no association with timerelation (through
\(\Delta{t}_{\mathrm{EVP}}\)). With the setting of
SEAICEuseEVPstar to .TRUE.
(default), this form of EVP is used.
Using the terminology of Kimmritz et al. 2015 [KDL15], the evolution
equations of stress \(\sigma_{ij}\) and momentum \(\mathbf{u}\) can be
written as:
\(\mathbf{R}\) contains all terms in the momentum equations except for the rheology terms and the time derivative; \(\alpha\) and \(\beta\) are free parameters (SEAICE_evpAlpha, SEAICE_evpBeta) that replace the time stepping parameters SEAICE_deltaTevp (\(\Delta{t}_{\mathrm{EVP}}\)), SEAICE_elasticParm (\(E_{0}\)), or SEAICE_evpTauRelax (\(T\)). \(\alpha\) and \(\beta\) determine the speed of convergence and the stability. Usually, it makes sense to use \(\alpha = \beta\), and SEAICEnEVPstarSteps \(\gg (\alpha,\,\beta)\) (Kimmritz et al. 2015 [KDL15]). Currently, there is no termination criterion and the number of mEVP iterations is fixed to SEAICEnEVPstarSteps.
In order to use mEVP in MITgcm, compile with both #define
SEAICE_CGRID and #define
SEAICE_ALLOW_EVP in
SEAICE_OPTIONS.h (default) and make
sure that SEAICEuseEVPstar = .TRUE.,
(default) in data.seaice
.
By default SEAICEuseEVPrev is set to .TRUE.
and the
actual form of equations (8.21) and (8.22) is used
with fewer implicit terms and the factor of \(e^{2}\) dropped in the stress
equations (8.18) and (8.19). Although
this modifies the original EVP equations, it turns out to improve convergence
(Bouillon et al. 2013 [BFLM13]).
The aEVP scheme is an enhanced variant of mEVP (Kimmritz et al. 2016 [KDL16]), where the value of \(\alpha\) is set dynamically based on the stability criterion
with the grid cell area \(A_c\) and the ice and snow mass \(m\). This choice sacrifices speed of convergence for stability with the result that aEVP converges quickly to VP where \(\alpha\) can be small and more slowly in areas where the equations are stiff. In practice, aEVP leads to an overall better convergence than mEVP (Kimmritz et al. 2016 [KDL16]). To use aEVP in MITgcm set SEAICEaEVPcoeff \(= \tilde{c}\) (see (8.23); default is unset); this also sets the default values of SEAICEaEVPcStar (\(c=4\)) and SEAICEaEVPalphaMin (\(\alpha_{\min}=5\)). Good convergence has been obtained with these values (Kimmritz et al. 2016 [KDL16]):
SEAICEaEVPcoeff = 0.5,
SEAICEnEVPstarSteps = 500,
# The following two parameters are required by mEVP and aEVP,
# but they are TRUE by default:
SEAICEuseEVPstar = .TRUE.,
SEAICEuseEVPrev = .TRUE.,
Because of the Cgrid staggering of velocities and stresses, mEVP may not converge as successfully as in Kimmritz et al. (2015) [KDL15], see also Kimmritz et al. (2016) [KDL16]. Convergence at very high resolution (order 5 km) has not yet been studied.
8.6.2.5.5. IceOcean stress¶
Moving sea ice exerts a stress on the ocean which is the opposite of the stress
\(\mathbf{\tau}_\mathrm{ocean}\) in (8.3). This stress is
applied directly to the surface layer of the ocean model. An alternative ocean
stress formulation is given by Hibler and Bryan (1987)
[HB87]. Rather than applying \(\mathbf{\tau}_\mathrm{ocean}\)
directly, the stress is derived from integrating over the ice thickness to the
bottom of the oceanic surface layer. In the resulting equation for the
combined oceanice momentum, the interfacial stress cancels and the total
stress appears as the sum of windstress and divergence of internal ice
stresses: \(\delta(z) (\mathbf{\tau}_\mathrm{air} + \mathbf{F})/\rho_0\),
see also Eq. (2) of Hibler and Bryan (1987) [HB87]. The disadvantage
of this formulation is that now the velocity in the surface layer of the ocean
that is used to advect tracers, is really an average over the ocean surface
velocity and the ice velocity leading to an inconsistency as the ice
temperature and salinity are different from the oceanic variables. To turn on
the stress formulation of Hibler and Bryan (1987) [HB87], set
useHB87StressCoupling =.TRUE.,
, in data.seaice
.
8.6.2.5.6. Finitevolume discretization of the stress tensor divergence¶
On an Arakawa C grid, ice thickness and concentration and thus ice strength \(P\) and bulk and shear viscosities \(\zeta\) and \(\eta\) are naturally defined a Cpoints in the center of the grid cell. Discretization requires only averaging of \(\zeta\) and \(\eta\) to vorticity or Zpoints (or \(\zeta\)points, but here we use Z in order avoid confusion with the bulk viscosity) at the bottom left corner of the cell to give \(\overline{\zeta}^{Z}\) and \(\overline{\eta}^{Z}\). In the following, the superscripts indicate location at Z or C points, distance across the cell (F), along the cell edge (G), between \(u\)points (U), \(v\)points (V), and Cpoints (C). The control volumes of the \(u\) and \(v\)equations in the grid cell at indices \((i,j)\) are \(A_{i,j}^{w}\) and \(A_{i,j}^{s}\), respectively. With these definitions (which follow the model code documentation except that \(\zeta\)points have been renamed to Zpoints), the strain rates are discretized as:
so that the diagonal terms of the strain rate tensor are naturally defined at Cpoints and the symmetric offdiagonal term at Zpoints. Noslip boundary conditions (\(u_{i,j1}+u_{i,j}=0\) and \(v_{i1,j}+v_{i,j}=0\) across boundaries) are implemented via “ghostpoints”; for free slip boundary conditions \((\epsilon_{12})^Z=0\) on boundaries.
For a spherical polar grid, the coefficients of the metric terms are \(k_{1}=0\) and \(k_{2}=\tan\phi/a\), with the spherical radius \(a\) and the latitude \(\phi\); \(\Delta{x}_1 = \Delta{x} = a\cos\phi \Delta\lambda\), and \(\Delta{x}_2 = \Delta{y}=a\Delta\phi\). For a general orthogonal curvilinear grid, \(k_{1}\) and \(k_{2}\) can be approximated by finite differences of the cell widths:
The stress tensor is given by the constitutive viscousplastic relation \(\sigma_{\alpha\beta} = 2\eta\dot{\epsilon}_{\alpha\beta} + [(\zeta\eta)\dot{\epsilon}_{\gamma\gamma}  P/2 ]\delta_{\alpha\beta}\) . The stress tensor divergence \((\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}\), is discretized in finite volumes . This conveniently avoids dealing with further metric terms, as these are “hidden” in the differential cell widths. For the \(u\)equation (\(\alpha=1\)) we have:
with
Similarly, we have for the \(v\)equation (\(\alpha=2\)):
with
Again, noslip boundary conditions are realized via ghost points and \(u_{i,j1}+u_{i,j}=0\) and \(v_{i1,j}+v_{i,j}=0\) across boundaries. For freeslip boundary conditions the lateral stress is set to zeros. In analogy to \((\epsilon_{12})^Z=0\) on boundaries, we set \(\sigma_{21}^{Z}=0\), or equivalently \(\eta_{i,j}^{Z}=0\), on boundaries.
8.6.2.6. Thermodynamics¶
NOTE: THIS SECTION IS STILL NOT COMPLETE
In its original formulation the sea ice model uses simple 0layer thermodynamics following the appendix of Semtner (1976) [Sem76]. This formulation neglects storage of heat, that is, the heat capacity of ice is zero, and all internal heat sources so that the heat equation reduces to a constant conductive heat flux. This constant upward conductive heat flux together with a constant ice conductivity implies a linear temperature profile. The boundary conditions for the heat equations are: at the bottom of the ice \(T_{\rm bottom} = T_{\rm fr}\) (freezing point temperature of sea water), and at the surface: \(Q_{\rm top} = \frac{\partial{T}}{\partial{z}} = (K/h)(T_{0}T_{\rm fr})\), where \(K\) is the ice conductivity, \(h\) the ice thickness, and \(T_{0}T_{\rm fr}\) the difference between the ice surface temperature and the water temperature at the bottom of the ice (at the freezing point). The surface heat flux \(Q_{\rm top}\) is computed in a similar way to that of Parkinson and Washington (1979) [PW79] and Manabe et al. (1979) [MBS79]. The resulting equation for surface temperature is
where \(\epsilon\) is the emissivity of the surface (snow or ice), \(Q_{\rm S/LW\downarrow}\) the downwelling shortwave and longwave radiation to be prescribed, and \(Q_{\rm LW\uparrow}=\epsilon\sigma_B T_{0}^4\) the emitted long wave radiation with the StefanBoltzmann constant \(\sigma_B\). With explicit expressions in \(T_0\) for the turbulent fluxes of latent and sensible heat
(8.24) can be solved for \(T_0\) with an iterative
RalphsonNewton method, which usually converges very quickly in less that 10
iterations. In these equations, \(\rho_\mathrm{air}\) is the air density
(parameter SEAICE_rhoAir), \(C_E\) is the iceocean transfer
coefficient for sensible and latent heat (parameter SEAICE_dalton),
\(\Lambda_v\) and \(\Lambda_f\) are the latent heat of vaporization and
fusion, respectively (parameters SEAICE_lhEvap and
SEAICE_lhFusion), and \(c_p\) is the specific heat of air
(parameter SEAICE_cpAir). For the latent heat \(Q_{\rm LH}\) a
choice can be made between the old polynomial expression for saturation
humidity \(q_\mathrm{sat}(T_0)\) (by setting
useMaykutSatVapPoly to .TRUE.
) and the default exponential
relation approximation that is more accurate at low temperatures.
In the zerolayer model of Semtner (1976) [Sem76], the conductive
heat flux depends strongly on the ice thickness \(h\). However, the ice
thickness in the model represents a mean over a potentially very heterogeneous
thickness distribution. In order to parameterize a subgrid scale distribution
for heat flux computations, the mean ice thickness \(h\) is split into
\(N\) thickness categories \(H_{n}\) that are equally distributed
between \(2h\) and a minimum imposed ice thickness of \(5\,\text{cm}\)
by \(H_n= \frac{2n1}{7}\,h\) for \(n\in[1,N]\). The heat fluxes
computed for each thickness category are areaaveraged to give the total heat
flux (see Hibler 1984 [Hib84]). To use this thickness category
parameterization set SEAICE_multDim to the number of desired
categories in data.seaice
(7 is a good guess, for anything larger than 7
modify SEAICE_SIZE.h). Note that this
requires different restart files and switching this flag on in the middle of an
integration is not advised. As an alternative to the flat distribution, the
runtime parameter SEAICE_PDF (1Darray of lenght nITD)
can be used to prescribe an arbitrary distribution of ice thicknesses, for
example derived from observed distributions (CastroMorales et al. 2014
[CMKL+14]). In order to include the ice thickness distribution
also for snow, set SEAICE_useMultDimSnow to .TRUE.
(this is the
default); only then, the parameterization of always having a fraction of thin
ice is efficient and generally thicker ice is produced (see CastroMorales et
al. 2014 [CMKL+14]).
The atmospheric heat flux is balanced by an oceanic heat flux from below. The oceanic flux is proportional to \(\rho\,c_{p}\left(T_{w}T_{fr}\right)\) where \(\rho\) and \(c_{p}\) are the density and heat capacity of sea water and \(T_{\rm fr}\) is the local freezing point temperature that is a function of salinity. This flux is not assumed to instantaneously melt or create ice, but a time scale of three days (runtime parameter SEAICE_gamma_t) is used to relax \(T_{w}\) to the freezing point. The parameterization of lateral and vertical growth of sea ice follows that of Hibler (1979) and Hibler (1980) [Hib79, Hib80]; the socalled lead closing parameter \(h_{0}\) (runtime parameter HO) has a default value of 0.5 meters.
On top of the ice there is a layer of snow that modifies the heat flux and the albedo (Zhang et al. 1998 [ZWDHSR98]). Snow modifies the effective conductivity according to
where \(K_s\) is the conductivity of snow and \(h_s\) the snow
thickness. If enough snow accumulates so that its weight submerges the ice and
the snow is flooded, a simple mass conserving parameterization of snowice
formation (a floodfreeze algorithm following Archimedes’ principle) turns snow
into ice until the ice surface is back at \(z=0\) (see Leppäranta 1983
[Lepparanta83]). The floodfreeze algorithm is turned on with runtime
parameter SEAICEuseFlooding set to .TRUE.
.
8.6.2.6.1. Advection of thermodynamic variables¶
Effective ice thickness (ice volume per unit area, \(c h\)), concentration \(c\) and effective snow thickness (\(c h_s\)) are advected by ice velocities:
where \(\Gamma_X\) are the thermodynamic source terms and \(D_{X}\) the diffusive terms for quantities \(X= c h, c, c h_s\). From the various advection schemes that are available in MITgcm, we recommend fluxlimited schemes to preserve sharp gradients and edges that are typical of sea ice distributions and to rule out unphysical over and undershoots (negative thickness or concentration). These schemes conserve volume and horizontal area and are unconditionally stable, so that we can set \(D_{X}=0\). Runtime flags: SEAICEadvScheme (default=77, is a 2ndorder flux limited scheme), DIFF1 = \(D_{X}/\Delta{x}\) (default=0).
The MITgcm sea ice model provides the option to use the thermodynamics model of
Winton (2000) [Win00], which in turn is based on the 3layer model of
Semtner (1976) [Sem76] which treats brine content by means of
enthalpy conservation; the corresponding package thsice is described in section Section 8.6.1. This
scheme requires additional state variables, namely the enthalpy of the two ice
layers (instead of effective ice salinity), to be advected by ice
velocities. The internal sea ice temperature is inferred from ice enthalpy. To
avoid unphysical (negative) values for ice thickness and concentration, a
positive 2ndorder advection scheme with a SuperBee flux limiter (Roe 1985
[Roe85]) should be used to advect all seaicerelated quantities of the
Winton (2000) [Win00] thermodynamic model (runtime flag
thSIceAdvScheme \(= 77\) and thSIce_diffK \(=
D_{X} = 0\) in data.ice
, defaults are 0). Because of the nonlinearity of the
advection scheme, care must be taken in advecting these quantities: when simply
using ice velocity to advect enthalpy, the total energy (i.e., the volume
integral of enthalpy) is not conserved. Alternatively, one can advect the
energy content (i.e., product of icevolume and enthalpy) but then false
enthalpy extrema can occur, which then leads to unrealistic ice temperature. In
the currently implemented solution, the seaice mass flux is used to advect the
enthalpy in order to ensure conservation of enthalpy and to prevent false
enthalpy extrema.
8.6.2.6.2. Dynamical Ice Thickness Distribution (ITD)¶
The ice thickness distribution model used by MITgcm follows the implementation in the Los Alamos sea ice model CICE (https://github.com/CICEConsortium/CICE). There are two parts to it that are closely connected: the participation and ridging functions that determine which thickness classes take part in ridging and which thickness classes receive ice during ridging based on Thorndike et al. (1975) [TRMC75], and the ice strength parameterization by Rothrock (1975) [Rot75] which uses this information. The following description is slightly modified from Ungermann et al. (2017) [UTML17]. Verification experiment seaice_itd uses the ITD model.
Distribution, participation and redistribution functions in ridging¶
When SEAICE_ITD is defined in SEAICE_OPTIONS.h, the ice thickness is described by the ice thickness distribution \(g(h,\mathbf{x},t)\) for the subgridscale (see Thorndike et al. 1975 [TRMC75]), a probability density function for thickness \(h\) following the evolution equation
Here \(f=\frac{\mathrm{d} h}{\mathrm{d} t}\) is the thermodynamic growth rate and \(\Psi\) a function describing the mechanical redistribution of sea ice during ridging or lead opening.
The mechanical redistribution function \(\Psi\) generates open water in
divergent motion and creates ridged ice during convergent motion. The ridging
process depends on total strain rate and on the ratio between shear (runtime
parameter SEAICEshearParm) and divergent strain. In the single
category model, ridge formation is treated implicitly by limiting the ice
concentration to a maximum of one (see Hibler 1979 [Hib79]), so that
further volume increase in convergent motion leads to thicker ice. (This is
also the default for ITD models; to change from the default, set runtime
parameter SEAICEsimpleRidging =.FALSE.,
in data.seaice
). For
the ITD model, the ridging mode in convergence
gives the effective change for the ice volume with thickness between \(h\) and \(h+\textrm{d} h\) as the normalized difference between the ice \(n(h)\) generated by ridging and the ice \(a(h)\) participating in ridging.
The participation function \(a(h) = b(h)g(h)\) can be computed either following Thorndike et al. (1975) [TRMC75] (runtime parameter SEAICEpartFunc =0) or Lipscomb et al. (2007) [LHMJ07] (SEAICEpartFunc =1), and similarly the ridging function \(n(h)\) can be computed following Hilber (1980) [Hib80] (runtime parameter SEAICEredistFunc =0) or Lipscomb et al. (2007) [LHMJ07] (SEAICEredistFunc =1). As an example, we show here the functions that Lipscomb et al. (2007) [LHMJ07] suggested to avoid noise in the solutions. These functions are smooth and avoid nondifferentiable discontinuities, but so far we did not find any noise issues as in Lipscomb et al. (2007) [LHMJ07].
With SEAICEpartFunc = 1,
in data.seaice
, the participation
function with the relative amount of ice of thickness \(h\) weighted by an
exponential function
where \(G(h)=\int_0^h g(h) \textrm{d} h\) is the cumulative thickness distribution function, \(b_0\) a normalization factor, and \(a^*\) (SEAICEaStar) the exponential constant that determines which relative amount of thicker and thinner ice take part in ridging.
With SEAICEredistFunc = 1,
in data.seaice
, the ice generated by
ridging is calculated as
where the density function \(\gamma(h_1,h)\) of resulting thickness \(h\) for ridged ice with an original thickness of \(h_1\) is taken as
for \(h \geq h_{\min}\), with \(\gamma(h_1,h)=0\) for \(h < h_{\min}\). In this parameterization, the normalization factor \(k=\frac{h_{\min} + \lambda}{h_1}\), the efolding scale \(\lambda = \mu h_1^{1/2}\) and the minimum ridge thickness \(h_{\min}=\min(2h_1,h_1 + h_{\textrm{raft}})\) all depend on the original thickness \(h_1\). The maximal ice thickness allowed to raft \(h_{\textrm{raft}}\) is constant (SEAICEmaxRaft, default =1 m) and \(\mu\) (SEAICEmuRidging) is a tunable parameter.
In the numerical model these equations are discretized into a set of \(n\)
(nITD defined in SEAICE_SIZE.h) thickness categories employing the delta function
scheme of Bitz et al. (2001) [BHWE01]. For each thickness category in
an ITD configuration, the volume conservation equation (8.25) is
evaluated using the heat flux with the categoryspecific values for ice and
snow thickness, so there are no conceptual differences in the thermodynamics
between the single category and ITD configurations. The only difference is
that only in the thinnest category the creation of new ice of thickness
\(H_0\) (runtime parameter HO) is possible, all other
categories are limited to basal growth. The conservation of ice area is
replaced by the evolution equation of the ITD (8.26) that is discretized
in thickness space with \(n+1\) category limits given by runtime parameter
Hlimit. If Hlimit is not set in data.seaice
, a
simple recursive formula following Lipscomb (2001) [Lip01] is used
to compute Hlimit:
with \(H_\mathrm{limit}(0)=0\) m and \(H_\mathrm{limit}(n)=999.9\) m. The three constants are the runtime parameters Hlimit_c1, Hlimit_c2, and Hlimit_c3. The total ice concentration and volume can then be calculated by summing up the values for each category.
Ice strength parameterization¶
In the default approach of equation (8.5), the ice strength is parameterized following Hibler (1979) [Hib79] and \(P\) depends only on average ice concentration and thickness per grid cell and the constant ice strength parameters \(P^{\ast}\) (SEAICE_strength) and \(C^{\ast}\) (SEAICE_cStar). With an ice thickness distribution, it is possible to use a different parameterization following Rothrock (1975) [Rot75]
by considering the production of potential energy and the frictional energy
loss in ridging. The physical constant \(C_p = \rho_i (\rho_w  \rho_i)
\hat{g} / (2 \rho_w)\) is a combination of the gravitational acceleration
\(\hat{g}\) and the densities \(\rho_i\), \(\rho_w\) of ice and
water, and \(C_f\) (SEAICE_cf) is a scaling factor relating the
amount of work against gravity necessary for ridging to the amount of work
against friction. To calculate the integral, this parameterization needs
information about the ITD in each grid cell, while the default
parameterization (8.5) can be used for both ITD and single
thickness category models. In contrast to (8.5), which is based
on the plausible assumption that thick and compact ice is stronger than thin
and loose drifting ice, this parameterization (8.27) clearly
contains the more physical assumptions about energy conservation. For that
reason alone this parameterization is often considered to be more physically
realistic than (8.5), but in practice, the success is not so
clear (Ungermann et al. 2007 [UTML17]). Ergo, the default is to
use (8.5); set useHibler79IceStrength =.FALSE.,
in
data.seaice
to change this behavior.
8.6.2.7. Known issues and workarounds¶
An often encountered problem in long simulations with sea ice models is (local) perpetually increasing sea ice (plus snow) height; this is problematic when using a nonlinear free surface and useRealFreshWaterFlux set to
.TRUE.
, because the mass of the sea ice places a load on the sea surface, which if too large, can cause the surface cells of the model to become too thin so that the model eventually stops with an error message. Usually this problem occurs because of dynamical ice growth (i.e., convergence and ridging of ice) or simply too much net precipitation with insufficient summer surface melting. If the problem is dynamical in nature (e.g., caused by ridging in a deep inlet), the first step to try is to turn off the replacement pressure method (SEAICEpressReplFac = 0; in Section 8.6.2.5.1); turning this off provides resistance against additional growth due to further ridging, because the ice pressure \(P\) is no longer reduced as \(\Delta\rightarrow 0\) in nearly motionless thick ice (8.6). If this does not solve the problem, a somewhat more radical yet effective approach is simply to cap the sea ice load on the free surface by defining the CPP option SEAICE_CAP_ICELOAD. This option effectively limits the sea ice load (variable sIceLoad) to a mass of 1/5 of the the top grid cell depth. If desired, this limit can be changed in routine seaice_growth.F where variable heffTooHeavy is assigned.
8.6.2.8. Key subroutines¶
Toplevel routine: pkg/seaice/seaice_model.F
C !CALLING SEQUENCE:
c ...
c seaice_model (TOP LEVEL ROUTINE)
c 
c  #ifdef SEAICE_CGRID
c  SEAICE_DYNSOLVER
c  
c   < compute proxy for geostrophic velocity >
c  
c   < set up mass per unit area and Coriolis terms >
c  
c   < dynamic masking of areas with no ice >
c  
c  
c  #ELSE
c  DYNSOLVER
c  #ENDIF
c 
c  if ( useOBCS )
c  OBCS_APPLY_UVICE
c 
c  if ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow .OR. SEAICEadvSalt )
c  SEAICE_ADVDIFF
c 
c  SEAICE_REG_RIDGE
c 
c  if ( usePW79thermodynamics )
c  SEAICE_GROWTH
c 
c  if ( useOBCS )
c  if ( SEAICEadvHeff ) OBCS_APPLY_HEFF
c  if ( SEAICEadvArea ) OBCS_APPLY_AREA
c  if ( SEAICEadvSALT ) OBCS_APPLY_HSALT
c  if ( SEAICEadvSNOW ) OBCS_APPLY_HSNOW
c 
c  < do various exchanges >
c 
c  < do additional diagnostics >
c 
c o
8.6.2.9. SEAICE diagnostics¶
Diagnostics output is available via the diagnostics package (see Section 9.1). Available output fields are summarized in the following table:
+++
<Name>< grid >< Units >< Tile (max=80c)
+++
sIceLoadSM U1kg/m^2 seaice loading (in Mass of ice+snow / area unit)

SEA ICE STATE:

SIarea SM M1m^2/m^2 SEAICE fractional icecovered area [0 to 1]
SIheff SM M1m SEAICE effective ice thickness
SIhsnow SM M1m SEAICE effective snow thickness
SIhsalt SM M1g/m^2 SEAICE effective salinity
SIuice UU M1m/s SEAICE zonal ice velocity, >0 from West to East
SIvice VV M1m/s SEAICE merid. ice velocity, >0 from South to North

ATMOSPHERIC STATE AS SEEN BY SEA ICE:

SItices SM C M1K Surface Temperature over SeaIce (area weighted)
SIuwind UM U1m/s SEAICE zonal 10m wind speed, >0 increases uVel
SIvwind VM U1m/s SEAICE meridional 10m wind speed, >0 increases uVel
SIsnPrcpSM U1kg/m^2/s Snow precip. (+=dw) over SeaIce (area weighted)

FLUXES ACROSS ICEOCEAN INTERFACE (ATMOS to OCEAN FOR ICEFREE REGIONS):

SIfu UU U1N/m^2 SEAICE zonal surface wind stress, >0 increases uVel
SIfv VV U1N/m^2 SEAICE merid. surface wind stress, >0 increases vVel
SIqnet SM U1W/m^2 Ocean surface heatflux, turb+rad, >0 decreases theta
SIqsw SM U1W/m^2 Ocean surface shortwave radiat., >0 decreases theta
SIempmr SM U1kg/m^2/s Ocean surface freshwater flux, > 0 increases salt
SIqneto SM U1W/m^2 Open Ocean Part of SIqnet, turb+rad, >0 decr theta
SIqneti SM U1W/m^2 Ice Covered Part of SIqnet, turb+rad, >0 decr theta

FLUXES ACROSS ATMOSPHEREICE INTERFACE (ATMOS to OCEAN FOR ICEFREE REGIONS):

SIatmQntSM U1W/m^2 Net atmospheric heat flux, >0 decreases theta
SIatmFW SM U1kg/m^2/s Net freshwater flux from atmosphere & land (+=down)
SIfwSublSM U1kg/m^2/s Freshwater flux of sublimated ice, >0 decreases ice

THERMODYNAMIC DIAGNOSTICS:

SIareaPRSM M1m^2/m^2 SIarea preceeding ridging process
SIareaPTSM M1m^2/m^2 SIarea preceeding thermodynamic growth/melt
SIheffPTSM M1m SIheff preceeeding thermodynamic growth/melt
SIhsnoPTSM M1m SIhsnow preceeeding thermodynamic growth/melt
SIaQbOCNSM M1m/s Potential HEFF rate of change by ocean ice flux
SIaQbATCSM M1m/s Potential HEFF rate of change by atm flux over ice
SIaQbATOSM M1m/s Potential HEFF rate of change by open ocn atm flux
SIdHbOCNSM M1m/s HEFF rate of change by ocean ice flux
SIdSbATCSM M1m/s HSNOW rate of change by atm flux over sea ice
SIdSbOCNSM M1m/s HSNOW rate of change by ocean ice flux
SIdHbATCSM M1m/s HEFF rate of change by atm flux over sea ice
SIdHbATOSM M1m/s HEFF rate of change by open ocn atm flux
SIdHbFLOSM M1m/s HEFF rate of change by flooding snow
SIdAbATOSM M1m^2/m^2/s Potential AREA rate of change by open ocn atm flux
SIdAbATCSM M1m^2/m^2/s Potential AREA rate of change by atm flux over ice
SIdAbOCNSM M1m^2/m^2/s Potential AREA rate of change by ocean ice flux
SIdA SM M1m^2/m^2/s AREA rate of change (net)

DYNAMIC/RHEOLOGY DIAGNOSTICS:

SIpress SM M1N/m SEAICE strength (with upper and lower limit)
SIzeta SM M1kg/s SEAICE nonlinear bulk viscosity
SIeta SM M1kg/s SEAICE nonlinear shear viscosity
SIsig1 SM M1no units SEAICE normalized principle stress, component one
SIsig2 SM M1no units SEAICE normalized principle stress, component two
SIshear SM M11/s SEAICE shear deformation rate
SIdelta SM M11/s SEAICE Delta deformation rate
SItensilSM M1N/m SEAICE maximal tensile strength

ADVECTIVE/DIFFUSIVE FLUXES OF SEA ICE variables:

ADVxHEFFUU M1m.m^2/s Zonal Advective Flux of eff ice thickn
ADVyHEFFVV M1m.m^2/s Meridional Advective Flux of eff ice thickn
SIuheff UU M1m^2/s Zonal Transport of eff ice thickn (centered)
SIvheff VV M1m^2/s Meridional Transport of eff ice thickn (centered)
DFxEHEFFUU M1m^2/s Zonal Diffusive Flux of eff ice thickn
DFyEHEFFVV M1m^2/s Meridional Diffusive Flux of eff ice thickn
ADVxAREAUU M1m^2/m^2.m^2/s Zonal Advective Flux of fract area
ADVyAREAVV M1m^2/m^2.m^2/s Meridional Advective Flux of fract area
DFxEAREAUU M1m^2/m^2.m^2/s Zonal Diffusive Flux of fract area
DFyEAREAVV M1m^2/m^2.m^2/s Meridional Diffusive Flux of fract area
ADVxSNOWUU M1m.m^2/s Zonal Advective Flux of eff snow thickn
ADVySNOWVV M1m.m^2/s Meridional Advective Flux of eff snow thickn
DFxESNOWUU M1m.m^2/s Zonal Diffusive Flux of eff snow thickn
DFyESNOWVV M1m.m^2/s Meridional Diffusive Flux of eff snow thickn
ADVxSSLTUU M1(g/kg).m^2/s Zonal Advective Flux of seaice salinity
ADVySSLTVV M1(g/kg).m^2/s Meridional Advective Flux of seaice salinity
DFxESSLTUU M1(g/kg).m^2/s Zonal Diffusive Flux of seaice salinity
DFyESSLTVV M1(g/kg).m^2/s Meridional Diffusive Flux of seaice salinity
8.6.2.10. Experiments and tutorials that use seaice¶
verification/lab_sea: Labrador Sea experiment
verification/seaice_obcs, based on lab_sea
verification/offline_exf_seaice, idealized topography in a zonally reentrant channel, tests solvers and rheologies
verification/seaice_itd, based on offline_exf_seaice, tests ice thickness distribution
verification/global_ocean.cs32x15, global cubedsphereexperiment with combinations of pkg/seaice and pkg/thsice
verification/1D_ocean_ice_column, just thermodynamics