8.6.2. SEAICE Package

Authors: Martin Losch, Dimitris Menemenlis, An Nguyen, Jean-Michel Campin, Patrick Heimbach, Chris Hill and Jinlun Zhang

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). Run-time options, flags, filenames and field-related dates/times are set in data.seaice (Section 8.6.2.3). A description of key subroutines is given in Section 8.6.2.5. Available diagnostics output is listed in Section 8.6.2.6.

8.6.2.2. SEAICE configuration and compiling

8.6.2.2.1. Compile-time 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 adding seaice to it
  • or using genmake2 adding -enable=seaice or -disable=seaice switches
  • required 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.16 summarizes the most important ones. For more options see SEAICE_OPTIONS.h.

Table 8.16 Some of the most relevant CPP preporocessor flags in the seaice package.
CPP option Default Description”
SEAICE_DEBUG #undef enhance STDOUT for debugging
SEAICE_ALLOW_DYNAMICS #define sea ice dynamics code
SEAICE_CGRID #define LSR solver on C-grid (rather than original B-grid)
SEAICE_ALLOW_EVP #define enable use of EVP rheology solver
SEAICE_ALLOW_JFNK #define enable use of JFNK rheology solver
SEAICE_ALLOW_KRYLOV #define enable use of Krylov rheology solver
SEAICE_LSR_ZEBRA #undef use a coloring method for LSR solver
SEAICE_EXTERNAL_FLUXES #define use pkg/exf-computed fluxes as starting point
SEAICE_ZETA_SMOOTHREG #define use differentiable regularization for viscosities
SEAICE_DELTA_SMOOTHREG #undef use differentiable regularization for \(1/\Delta\)
SEAICE_ALLOW_BOTTOMDRAG #undef enable grounding parameterization for improved fastice in shallow seas
SEAICE_ITD #undef run with dynamical sea Ice Thickness Distribution (ITD)
SEAICE_VARIABLE_SALINITY #undef enable sea ice with variable salinity
ALLOW_SITRACER #undef enable sea ice tracer package
SEAICE_BICE_STRESS #undef B-grid only for backward compatiblity: turn on ice-stress on ocean
EXPLICIT_SSH_SLOPE #undef B-grid only for backward compatiblity: use ETAN for tilt computations rather than geostrophic velocities

8.6.2.3. Run-time parameters

Run-time parameters (see Table 8.17) 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.17 lists most run-time parameters.

Table 8.17 Run-time parameters and default values
Name Default value Description
SEAICEwriteState FALSE write sea ice state to file
SEAICEuseDYNAMICS TRUE use dynamics
SEAICEuseJFNK FALSE use the JFNK-solver
SEAICEuseTEM FALSE use truncated ellipse method
SEAICEuseStrImpCpl FALSE use strength implicit coupling in LSR/JFNK
SEAICEuseMetricTerms TRUE use metric terms in dynamics
SEAICEuseEVPpickup TRUE use EVP pickups
SEAICEuseFluxForm TRUE use flux form for 2nd central difference advection scheme
SEAICErestoreUnderIce FALSE enable restoring to climatology under ice
SEAICEupdateOceanStress TRUE update ocean surface stress accounting for sea ice cover
SEAICEscaleSurfStress TRUE scale atmosphere and ocean-surface stress on ice by concentration (AREA)
SEAICEaddSnowMass TRUE in computing seaiceMass, add snow contribution
useHB87stressCoupling FALSE turn on ice-ocean stress coupling following
usePW79thermodynamics TRUE flag to turn off zero-layer-thermodynamics for testing
SEAICEadvHeff TRUE flag to turn off advection of scalar variable HEFF
SEAICEadvArea TRUE flag to turn off advection of scalar variable AREA
SEAICEadvSnow TRUE flag to turn off advection of scalar variable HSNOW
SEAICEadvSalt TRUE flag to turn off advection of scalar variable HSALT
SEAICEadvScheme 77 set advection scheme for seaice scalar state variables
SEAICEuseFlooding TRUE use flood-freeze algorithm
SEAICE_no_slip FALSE use no-slip boundary conditions instead of free-slip
SEAICE_deltaTtherm dTtracerLev (1) time step for seaice thermodynamics (s)
SEAICE_deltaTdyn dTtracerLev (1) time step for seaice dynamics (s)
SEAICE_deltaTevp 0.0 EVP sub-cycling time step (s); values \(>\) 0 turn on EVP
SEAICEuseEVPstar FALSE use modified EVP* instead of EVP, following [lemieux12]
SEAICEuseEVPrev FALSE “revisited form” variation on EVP*, following [bouillon13]
SEAICEnEVPstarSteps unset number of modified EVP* iterations
SEAICE_evpAlpha unset EVP* parameter (non-dim.), to replace 2*SEAICE_evpTauRelax/SEAICE_deltaTevp
SEAICE_evpBeta unset EVP* parameter (non-dim.), to replace SEAICE_deltaTdyn/SEAICE_deltaTevp
SEAICEaEVPcoeff unset largest stabilized frequency for adaptive EVP (non-dim.)
SEAICEaEVPcStar 4.0 aEVP multiple of stability factor (non-dim.), see [kimmritz16] \(\alpha * \beta = c^\ast * \gamma\)
SEAICEaEVPalphaMin 5.0 aEVP lower limit of alpha and beta (non-dim.), see [kimmritz16]
SEAICE_elasticParm 0.33333333 EVP parameter \(E_0\) (non-dim.), sets relaxation timescale SEAICE_evpTauRelaxtau = SEAICE_elasticParm * SEAICE_deltaTdyn
SEAICE_evpTauRelax dTtracerLev (1) * SEAICE_elasticParm relaxation time scale \(T\) for EVP waves (s)
SEAICE_OLx OLx - 2 overlap for LSR-solver or preconditioner, \(x\)-dimension
SEAICE_OLy OLy - 2 overlap for LSR-solver or preconditioner, \(y\)-dimension
SEAICEnonLinIterMax 2/10 maximum number of non-linear (outer loop) iterations (LSR/JFNK)
SEAICElinearIterMax 1500/10 maximum number of linear iterations (LSR/JFNK)
SEAICE_JFNK_lsIter (off) start line search after “lsIter” Newton iterations
SEAICEnonLinTol 1.0E-05 non-linear tolerance parameter for JFNK solver
JFNKgamma_lin_min 0.10 minimum tolerance parameter for linear JFNK solver
JFNKgamma_lin_max 0.99 maximum tolerance parameter for linear JFNK solver
JFNKres_tFac unset tolerance parameter for FGMRES residual
SEAICE_JFNKepsilon 1.0E-06 step size for the FD-gradient in s/r seaice_jacvec
SEAICE_dumpFreq dumpFreq dump frequency (s)
SEAICE_dump_mdsio TRUE write snapshot using /pkg/mdsio
SEAICE_dump_mnc FALSE write snapshot using /pkg/mnc
SEAICE_initialHEFF 0.0 initial sea ice thickness averaged over grid cell (m)
SEAICE_drag 1.0E-03 air-ice drag coefficient (non-dim.)
OCEAN_drag 1.0E-03 air-ocean drag coefficient (non-dim.)
SEAICE_waterDrag 5.5E-03 water-ice drag coefficient (non-dim.)
SEAICE_dryIceAlb 0.75 winter sea ice albedo
SEAICE_wetIceAlb 0.66 summer sea ice albedo
SEAICE_drySnowAlb 0.84 dry snow albedo
SEAICE_wetSnowAlb 0.70 wet snow albedo
SEAICE_waterAlbedo 0.10 water albedo (not used if #define SEAICE_EXTERNAL_FLUXES)
SEAICE_strength 2.75E+04 sea ice strength constant \(P^{\ast}\) (N/m2)
SEAICE_cStar 20.0 sea ice strength constant \(C^{\ast}\) (non-dim.)
SEAICE_rhoAir 1.3 (or pkg/exf value) density of air (kg/m3)
SEAICE_cpAir 1004.0 (or pkg/exf value) specific heat of air (J/kg/K)
SEAICE_lhEvap 2.5E+06 (or pkg/exf value) latent heat of evaporation (J/kg)
SEAICE_lhFusion 3.34E+05 (or pkg/exf value) latent heat of fusion (J/kg)
SEAICE_dalton 1.75E-03 ice-ocean transfer coefficient for latent and sensible heat (non-dim.)
useMaykutSatVapPoly FALSE use Maykut polynomial to compute saturation vapor pressure
SEAICE_iceConduct 2.16560E+00 sea ice conductivity (W m-1 K-1)
SEAICE_snowConduct 3.10000E-01 snow conductivity (W m-1 K-1)
SEAICE_emissivity 0.970018 (or pkg/exf value) longwave ocean surface emissivity (non-dim.)
SEAICE_snowThick 0.15 cutoff snow thickness to use snow albedo (m)
SEAICE_shortwave 0.30 ice penetration shortwave radiation factor (non-dim.)
SEAICE_saltFrac 0.0 salinity newly formed ice (as fraction of ocean surface salinity)
SEAICE_frazilFrac 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)
SEAICEstressFactor 1.0 scaling factor for ice area in computing total ocean stress (non-dim.)
HeffFile unset filename for initial sea ice eff. thickness field HEFF (m)
AreaFile unset filename for initial fraction sea ice cover AREA (non-dim.)
HsnowFile unset filename for initial eff. snow thickness field HSNOW (m)
HsaltFile unset filename for initial eff. sea ice salinity field HSALT (g/m2)
LSR_ERROR 1.0E-04 sets accuracy of LSR solver
DIFF1 0.0 parameter used in advect.F
HO 0.5 lead closing parameter \(h_0\) (m); demarcation thickness between thick and thin ice which determines partition between vertical and lateral ice growth
MIN_ATEMP -50.0 minimum air temperature (oC)
MIN_LWDOWN 60.0 minimum downward longwave (W/m2)
MIN_TICE -50.0 minimum ice temperature (oC)
IMAX_TICE 10 number of iterations for ice surface temperature solution
SEAICE_EPS 1.0E-10 a “small number” used in various routines
SEAICE_area_reg 1.0E-5 minimum concentration to regularize ice thickness
SEAICE_hice_reg 0.05 minimum ice thickness (m) for regularization
SEAICE_multDim 1 number of ice categories for thermodynamics
SEAICE_useMultDimSnow TRUE use same fixed pdf for snow as for multi-thickness-category ice

The following dynamical ice thickness distribution and ridging parameters in Table 8.18 are only active with #define SEAICE_ITD. All parameters are non-dimensional unless indicated.

Table 8.18 Thickness distribution and ridging parameters
Name Default value Description
useHibler79IceStrength TRUE use [hib79] ice strength; do not use [rothrock:75] with #define SEAICE_ITD
SEAICEsimpleRidging TRUE use simple ridging a la [hib79]
SEAICE_cf 17.0 scaling parameter of [rothrock:75] ice strength parameterization
SEAICEpartFunc 0 use partition function of [thorndike:75]
SEAICEredistFunc 0 use redistribution function of [hib80]
SEAICEridgingIterMax 10 maximum number of ridging sweeps
SEAICEshearParm 0.5 fraction of shear to be used for ridging
SEAICEgStar 0.15 max. ice conc. that participates in ridging [thorndike:75]
SEAICEhStar 25.0 ridging parameter for [thorndike:75], [lipscomb:07]
SEAICEaStar 0.05 similar to SEAICEgStar for [lipscomb:07] participation function
SEAICEmuRidging 3.0 similar to SEAICEhStar for [lipscomb:07] ridging function
SEAICEmaxRaft 1.0 regularization parameter for rafting
SEAICEsnowFracRidge 0.5 fraction of snow that remains on ridged ice
SEAICEuseLinRemapITD TRUE use linear remapping scheme of [lipscomb:01]
Hlimit unset nITD+1-array of ice thickness category limits (m)
Hlimit_c1, Hlimit_c2, Hlimit_c3 3.0, 15.0, 3.0 when Hlimit is not set, then these parameters determine Hlimit from a simple function following [lipscomb:01]

8.6.2.4. Description

The MITgcm sea ice model is based on a variant of the viscous-plastic (VP) dynamic-thermodynamic sea ice model (Zhang and Hibler 1997 [zhang97]) first introduced in Hibler (1979) and Hibler (1980) [hib79][hib80]. In order to adapt this model to the requirements of coupled ice-ocean state estimation, many important aspects of the original code have been modified and improved, see Losch et al. (2010) [losch:10]:

  • the code has been rewritten for an Arakawa C-grid, both B- and C-grid variants are available; the C-grid code allows for no-slip and free-slip lateral boundary conditions;
  • three different solution methods for solving the nonlinear momentum equations have been adopted: LSOR (Zhang and Hibler 1997 [zhang97]), EVP (Hunke and Dukowicz 1997 [hun97]), JFNK (Lemieux et al. 2010 [lemieux10], Losch et al. 2014 [losch:14]);
  • ice-ocean stress can be formulated as in Hibler and Bryan (1987) [hibler87] or as in Campin et al. (2008) [cam:08];
  • 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 large-scale climate studies are the viscous-plastic (VP) model (Hilber 1979 [hib79]), the cavitating fluid (CF) model (Flato and Hibler 1992 [fla92]), and the elastic-viscous-plastic (EVP) model (Hunke and Dukowicz 1997 [hun97]). 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) [zhang97] for use in a parallel configuration. An EVP model and a free-drift implementation can be selected with run-time flags.

8.6.2.4.1. Compatibility with ice-thermodynamics package pkg/thsice

By default pkg/seaice includes the original so-called zero-layer thermodynamics with a snow cover as in the appendix of Semtner (1976) [sem76]. The zero-layer 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]) three-layer thermodynamic model that permits heat storage in ice.

The Winton (2000) sea-ice 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 zero-layer 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 sea-ice model is global: in ice-free regions bulk formulae (by default computed in package exf) are used to estimate oceanic forcing from the atmospheric fields.

8.6.2.4.3. Dynamics

The momentum equation of the sea-ice model is

(8.2)\[m \frac{D\mathbf{u}}{Dt} = -mf\mathbf{k}\times\mathbf{u} + \mathbf{\tau}_\mathrm{air} + \mathbf{\tau}_\mathrm{ocean} - m \nabla{\phi(0)} + \mathbf{F}\]

where \(m=m_{i}+m_{s}\) is the ice and snow mass per unit area; \(\mathbf{u}=u\mathbf{i}+v\mathbf{j}\) is the ice velocity vector; \(\mathbf{i}\), \(\mathbf{j}\), and \(\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 wind-ice and ocean-ice 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 sea-ice momentum is neglected. The wind and ice-ocean stress terms are given by

\[\begin{split}\begin{aligned} \mathbf{\tau}_\mathrm{air} = & \rho_\mathrm{air} C_\mathrm{air} |\mathbf{U}_\mathrm{air} -\mathbf{u}| R_\mathrm{air} (\mathbf{U}_\mathrm{air} -\mathbf{u}) \\ \mathbf{\tau}_\mathrm{ocean} = & \rho_\mathrm{ocean}C_\mathrm{ocean} |\mathbf{U}_\mathrm{ocean}-\mathbf{u}| R_\mathrm{ocean}(\mathbf{U}_\mathrm{ocean}-\mathbf{u}) \end{aligned}\end{split}\]

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.4.4. Viscous-Plastic (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 viscous-plastic (VP) constitutive law:

(8.3)\[\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} + \left[\zeta(\dot{\epsilon}_{ij},P) - \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij} - \frac{P}{2}\delta_{ij}\]

The ice strain rate is given by

\[\dot{\epsilon}_{ij} = \frac{1}{2}\left( \frac{\partial{u_{i}}}{\partial{x_{j}}} + \frac{\partial{u_{j}}}{\partial{x_{i}}}\right)\]

The maximum ice pressure \(P_{\max}\), a measure of ice strength, depends on both thickness \(h\) and compactness (concentration) \(c\):

(8.4)\[P_{\max} = P^{\ast}c\,h\,\exp\{-C^{\ast}\cdot(1-c)\},\]

with the constants \(P^{\ast}\) (run-time parameter SEAICE_strength) and \(C^{\ast}\) (run-time parameter SEAICE_cStar). The nonlinear bulk and shear viscosities \(\eta\) and \(\zeta\) 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\) equal to \(2\); they are given by:

\[\begin{split}\begin{aligned} \zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})}, \zeta_{\max}\right) \\ \eta =& \frac{\zeta}{e^2} \end{aligned}\end{split}\]

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, run-time 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 run-time 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 CPP-flag SEAICE_ZETA_SMOOTHREG in SEAICE_OPTIONS.h before compiling replaces the method for bounding \(\zeta\) by a smooth (differentiable) expression:

(8.5)\[\begin{split}\begin{split} \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min}) \,\zeta_{\max}}\right)\\ &= \frac{P}{2\Delta^\ast} \tanh\left(\frac{\Delta^\ast}{\min(\Delta,\Delta_{\min})}\right) \end{split}\end{split}\]

where \(\Delta_{\min}=10^{-20}\,\text{s}^{-1}\) should be chosen to avoid divisions by zero.

8.6.2.4.5. LSR and JFNK solver

In matrix notation, the discretized momentum equations can be written as

(8.6)\[ \mathbf{A}(\mathbf{x})\,\mathbf{x} = \mathbf{b}(\mathbf{x}).\]

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.6) in the sea ice component of MITgcm is an iterative Picard solver: in the \(k\)-th iteration a linearized form \(\mathbf{A}(\mathbf{x}^{k-1})\,\mathbf{x}^{k} = \mathbf{b}(\mathbf{x}^{k-1})\) 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 pseudo-time steps can be controlled by the run-time parameter SEAICEnonLinIterMax (default is 2).

In order to overcome the poor convergence of the Picard-solver, Lemieux et al. (2010) [lemieux10] introduced a Jacobian-free Newton-Krylov solver for the sea ice momentum equations. This solver is also implemented in MITgcm (see Losch et al. 2014 [losch:14]). 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 (\(k-1\)) estimate \(\mathbf{x}^{k-1}\):

(8.7)\[\mathbf{F}(\mathbf{x}^{k-1}+\delta\mathbf{x}^{k}) = \mathbf{F}(\mathbf{x}^{k-1}) + \mathbf{F}'(\mathbf{x}^{k-1}) \,\delta\mathbf{x}^{k}\]

with the Jacobian \(\mathbf{J}\equiv\mathbf{F}'\). The root \(\mathbf{F}(\mathbf{x}^{k-1}+\delta\mathbf{x}^{k})=0\) is found by solving

(8.8)\[\mathbf{J}(\mathbf{x}^{k-1})\,\delta\mathbf{x}^{k} = -\mathbf{F}(\mathbf{x}^{k-1})\]

for \(\delta\mathbf{x}^{k}\). The next (\(k\)-th) estimate is given by \(\mathbf{x}^{k}=\mathbf{x}^{k-1}+a\,\delta\mathbf{x}^{k}\). In order to avoid overshoots the factor \(a\) is iteratively reduced in a line search (\(a=1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \ldots\)) until \(\|\mathbf{F}(\mathbf{x}^k)\| < \|\mathbf{F}(\mathbf{x}^{k-1})\|\), where \(\|\cdot\|=\int\cdot\,dx^2\) is the \(L_2\)-norm. In practice, the line search is stopped at \(a=\frac{1}{8}\). The line search starts after SEAICE_JFNK_lsIter nonlinear Newton iterations (off by default).

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.8). The action of \(\mathbf{J}\) can be approximated by a first-order Taylor series expansion:

(8.9)\[\mathbf{J}(\mathbf{x}^{k-1})\,\mathbf{w} \approx \frac{\mathbf{F}(\mathbf{x}^{k-1}+\epsilon\mathbf{w}) - \mathbf{F}(\mathbf{x}^{k-1})} \epsilon\]

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 right-hand side preconditioning to solve (8.8) 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}^{k-1})\) where \(\mathbf{x}^{k-1}\) is the estimate of the previous Newton step \(k-1\). The transformed equation (8.8) becomes

(8.10)\[\mathbf{J}(\mathbf{x}^{k-1})\,\mathbf{P}^{-1}\delta\mathbf{z} = -\mathbf{F}(\mathbf{x}^{k-1}), \quad\text{with} \quad \delta{\mathbf{z}} = \mathbf{P}\delta\mathbf{x}^{k}\]

The Krylov method iteratively improves the approximate solution to (8.10) 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}^{k-1}) -\mathbf{J}(\mathbf{x}^{k-1})\,\delta\mathbf{x}^{k}_{0}\) is the initial residual of (8.8); \(\mathbf{r}_0=-\mathbf{F}(\mathbf{x}^{k-1})\) 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}^{k-1})\), 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) [losch:14]).

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}}\) (run-time parameter SEAICEnonLinTol = 1.E-4, 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}^{k-1})\|\) where \(\gamma_k\) is determined by

(8.11)\[\begin{split}\gamma_k = \begin{cases} \gamma_0 &\text{for $\|\mathbf{F}(\mathbf{x}^{k-1})\| \geq r$}, \\ \max\left(\gamma_{\min}, \frac{\|\mathbf{F}(\mathbf{x}^{k-1})\|} {\|\mathbf{F}(\mathbf{x}^{k-2})\|}\right) &\text{for $\|\mathbf{F}(\mathbf{x}^{k-1})\| < r$,} \end{cases}\end{split}\]

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 = .TRUE., turns on “strength implicit coupling” (see Hutchings et al. 2004 [hutchings04]) 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.4.6. Elastic-Viscous-Plastic (EVP) Dynamics

Hunke and Dukowicz (1997) [hun97] introduced an elastic contribution to the strain rate in order to regularize (8.3) in such a way that the resulting elastic-viscous-plastic (EVP) and VP models are identical at steady state,

(8.12)\[\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + \frac{1}{2\eta}\sigma_{ij} + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij} + \frac{P}{4\zeta}\delta_{ij} = \dot{\epsilon}_{ij}.\]

The EVP-model uses an explicit time stepping scheme with a short timestep. According to the recommendation in Hunke and Dukowicz (1997) [hun97], the EVP-model 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.12) can be written as:

(8.13)\[\frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} + \frac{P}{2T} = \frac{P}{2T\Delta} D_D\]
(8.14)\[\frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T} = \frac{P}{2T\Delta} D_T\]
(8.15)\[\frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T} = \frac{P}{4T\Delta} D_S\]

Here, the elastic parameter \(E\) is redefined in terms of a damping timescale \(T\) for elastic waves

\[E=\frac{\zeta}{T}\]

\(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 and recommend.

To use the 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). The solver is turned on by setting the sub-cycling time step SEAICE_deltaTevp to a value larger than zero. The choice of this time step is under debate. Hunke and Dukowicz (1997) [hun97] 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.4.7. More stable variants of Elastic-Viscous-Plastic 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 [hun01][lemieux12][bouillon13]). This has led to a modified EVP or EVP* (Lemieux et al. 2012, Bouillon et a1. 2013, Kimmritz et al. 2015 [lemieux12][bouillon13][kimmritz15]); here, we refer to these variants by modified EVP (mEVP) and adaptive EVP (aEVP). The main idea is to modify the “natural” time-discretization of the momentum equations:

(8.16)\[m\frac{D\mathbf{u}}{Dt} \approx m\frac{\mathbf{u}^{p+1}-\mathbf{u}^{n}}{\Delta{t}} + \beta^{\ast}\frac{\mathbf{u}^{p+1}-\mathbf{u}^{p}}{\Delta{t}_{\mathrm{EVP}}}\]

where \(n\) is the previous time step index, and \(p\) is the previous sub-cycling 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 re-interpreted as a pure iterative solver where the sub-cycling has no association with time-relation (through \(\Delta{t}_{\mathrm{EVP}}\)). Using the terminology of Kimmritz et al. 2015 [kimmritz15], the evolution equations of stress \(\sigma_{ij}\) and momentum \(\mathbf{u}\) can be written as:

(8.17)\[\sigma_{ij}^{p+1}=\sigma_{ij}^p+\frac{1}{\alpha} \Big(\sigma_{ij}(\mathbf{u}^p)-\sigma_{ij}^p\Big), \phantom{\int}\]
(8.18)\[\mathbf{u}^{p+1}=\mathbf{u}^p+\frac{1}{\beta} \Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+ \frac{\Delta t}{m}\mathbf{R}^{p}+\mathbf{u}_n -\mathbf{u}^p\Big)\]

\(\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 [kimmritz15]). Currently, there is no termination criterion and the number of mEVP iterations is fixed to SEAICEnEVPstarSteps.

In order to use mEVP in MITgcm, set SEAICEuseEVPstar = .TRUE., in data.seaice. If SEAICEuseEVPrev =.TRUE., the actual form of equations (8.17) and (8.18) is used with fewer implicit terms and the factor of \(e^{2}\) dropped in the stress equations (8.14) and (8.15). Although this modifies the original EVP-equations, it turns out to improve convergence (Bouillon et al. 2013 [bouillon13]).

Another variant is the aEVP scheme (Kimmritz et al. 2016 [kimmritz16]), where the value of \(\alpha\) is set dynamically based on the stability criterion

(8.19)\[\alpha = \beta = \max\left( \tilde{c} \pi\sqrt{c \frac{\zeta}{A_{c}} \frac{\Delta{t}}{\max(m,10^{-4}\,\text{kg})}},\alpha_{\min} \right)\]

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 [kimmritz16]). To use aEVP in MITgcm set SEAICEaEVPcoeff \(= \tilde{c}\); 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 [kimmritz16]): SEAICEaEVPcoeff \(= 0.5\), SEAICEnEVPstarSteps \(= 500\), SEAICEuseEVPstar = .TRUE., SEAICEuseEVPrev = .TRUE..

Note, that probably because of the C-grid staggering of velocities and stresses, mEVP may not converge as successfully as in Kimmritz et al. (2015) [kimmritz15], see also Kimmritz et al. (2016) [kimmritz16], and that convergence at very high resolution (order 5 km) has not been studied yet.

8.6.2.4.8. Truncated ellipse method (TEM) for yield curve

In the so-called truncated ellipse method the shear viscosity \(\eta\) is capped to suppress any tensile stress:

(8.20)\[\eta = \min\left(\frac{\zeta}{e^2}, \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} {\sqrt{\max(\Delta_{\min}^{2},(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2 +4\dot{\epsilon}_{12}^2})}\right).\]

To enable this method, set #define SEAICE_ALLOW_TEM in SEAICE_OPTIONS.h and turn it on with SEAICEuseTEM in data.seaice.

8.6.2.4.9. Ice-Ocean stress

Moving sea ice exerts a stress on the ocean which is the opposite of the stress \(\mathbf{\tau}_\mathrm{ocean}\) in (8.2). 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) [hibler87]. 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 ocean-ice 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) [hibler87]. 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) [hibler87], set useHB87StressCoupling =.TRUE., in data.seaice.

8.6.2.4.10. Finite-volume 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 C-points in the center of the grid cell. Discretization requires only averaging of \(\zeta\) and \(\eta\) to vorticity or Z-points (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 C-points (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 Z-points), the strain rates are discretized as:

\[\begin{split}\begin{aligned} \dot{\epsilon}_{11} &= \partial_{1}{u}_{1} + k_{2}u_{2} \\ \notag => (\epsilon_{11})_{i,j}^C &= \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} + k_{2,i,j}^{C}\frac{v_{i,j+1}+v_{i,j}}{2} \\ \dot{\epsilon}_{22} &= \partial_{2}{u}_{2} + k_{1}u_{1} \\\notag => (\epsilon_{22})_{i,j}^C &= \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} + k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \dot{\epsilon}_{12} = \dot{\epsilon}_{21} &= \frac{1}{2}\biggl( \partial_{1}{u}_{2} + \partial_{2}{u}_{1} - k_{1}u_{2} - k_{2}u_{1} \biggr) \\ \notag => (\epsilon_{12})_{i,j}^Z &= \frac{1}{2} \biggl( \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^V} + \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^U} \\\notag &\phantom{=\frac{1}{2}\biggl(} - k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} - k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \biggr),\end{aligned}\end{split}\]

so that the diagonal terms of the strain rate tensor are naturally defined at C-points and the symmetric off-diagonal term at Z-points. No-slip boundary conditions (\(u_{i,j-1}+u_{i,j}=0\) and \(v_{i-1,j}+v_{i,j}=0\) across boundaries) are implemented via “ghost-points”; 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:

\[\begin{split}\begin{aligned} k_{1,i,j}^{C} &= \frac{1}{\Delta{y}_{i,j}^{F}} \frac{\Delta{y}_{i+1,j}^{G}-\Delta{y}_{i,j}^{G}}{\Delta{x}_{i,j}^{F}} \\ k_{2,i,j}^{C} &= \frac{1}{\Delta{x}_{i,j}^{F}} \frac{\Delta{x}_{i,j+1}^{G}-\Delta{x}_{i,j}^{G}}{\Delta{y}_{i,j}^{F}} \\ k_{1,i,j}^{Z} &= \frac{1}{\Delta{y}_{i,j}^{U}} \frac{\Delta{y}_{i,j}^{C}-\Delta{y}_{i-1,j}^{C}}{\Delta{x}_{i,j}^{V}} \\ k_{2,i,j}^{Z} &= \frac{1}{\Delta{x}_{i,j}^{V}} \frac{\Delta{x}_{i,j}^{C}-\Delta{x}_{i,j-1}^{C}}{\Delta{y}_{i,j}^{U}}\end{aligned}\end{split}\]

The stress tensor is given by the constitutive viscous-plastic 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:

\[\begin{split}\begin{aligned} (\nabla\sigma)_{1}: \phantom{=}& \frac{1}{A_{i,j}^w} \int_{\mathrm{cell}}(\partial_1\sigma_{11}+\partial_2\sigma_{21})\,dx_1\,dx_2 \\\notag =& \frac{1}{A_{i,j}^w} \biggl\{ \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{11}dx_2\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}} + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{21}dx_1\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}} \biggr\} \\ \notag \approx& \frac{1}{A_{i,j}^w} \biggl\{ \Delta{x}_2\sigma_{11}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}} + \Delta{x}_1\sigma_{21}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}} \biggr\} \\ \notag =& \frac{1}{A_{i,j}^w} \biggl\{ (\Delta{x}_2\sigma_{11})_{i,j}^C - (\Delta{x}_2\sigma_{11})_{i-1,j}^C \\\notag \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{} + (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z \biggr\}\end{aligned}\end{split}\]

with

\[\begin{split}\begin{aligned} (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+} \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j} \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag &+ \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j} k_{2,i,j}^C \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j} \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j} k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag \phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\ (\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+} \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j} \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag & + \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j} \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j} k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j} k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}\end{aligned}\end{split}\]

Similarly, we have for the \(v\)-equation (\(\alpha=2\)):

\[\begin{split}\begin{aligned} (\nabla\sigma)_{2}: \phantom{=}& \frac{1}{A_{i,j}^s} \int_{\mathrm{cell}}(\partial_1\sigma_{12}+\partial_2\sigma_{22})\,dx_1\,dx_2 \\\notag =& \frac{1}{A_{i,j}^s} \biggl\{ \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{12}dx_2\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}} + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{22}dx_1\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}} \biggr\} \\ \notag \approx& \frac{1}{A_{i,j}^s} \biggl\{ \Delta{x}_2\sigma_{12}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}} + \Delta{x}_1\sigma_{22}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}} \biggr\} \\ \notag =& \frac{1}{A_{i,j}^s} \biggl\{ (\Delta{x}_2\sigma_{12})_{i+1,j}^Z - (\Delta{x}_2\sigma_{12})_{i,j}^Z \\ \notag \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{} + (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C \biggr\} \end{aligned}\end{split}\]

with

\[\begin{split}\begin{aligned} (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+} \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j} \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\\notag & + \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j} \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\\notag &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j} k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\\notag & - \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j} k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag (\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+} \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j} \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag &+ \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j} k_{2,i,j}^{C} \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j} \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j} k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag & -\Delta{x}_{i,j}^{F} \frac{P}{2}\end{aligned}\end{split}\]

Again, no-slip boundary conditions are realized via ghost points and \(u_{i,j-1}+u_{i,j}=0\) and \(v_{i-1,j}+v_{i,j}=0\) across boundaries. For free-slip 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.4.11. Thermodynamics

NOTE: THIS SECTION IS STILL NOT COMPLETE

In its original formulation the sea ice model uses simple thermodynamics following the appendix of Semtner (1976) [sem76]. This formulation does not allow storage of heat, that is, the heat capacity of ice is zero. Upward conductive heat flux is parameterized assuming a linear temperature profile and together with a constant ice conductivity. It is expressed as \((K/h)(T_{w}-T_{0})\), where \(K\) is the ice conductivity, \(h\) the ice thickness, and \(T_{w}-T_{0}\) the difference between water and ice surface temperatures. This type of model is often refered to as a “zero-layer” model. The surface heat flux is computed in a similar way to that of Parkinson and Washington (1979) [parkinson:79] and Manabe et al. (1979) [manabe:79]. All fluxes are assumed to instanteously reach a stationary balance (\(\partial T_{0}/\partial t = 0\)):

(8.21)\[\begin{split}\begin{aligned} \rho\,c_p \frac{\partial T_{0}}{\partial t} = 0 =& \frac{K}{h}(T_{w}-T_{0}) + Q_{SW\downarrow}(1-\mathrm{albedo}) + \epsilon Q_{LW\downarrow} - Q_{LW\uparrow}(T_{0}) \\ &+ Q_{LH}(T_{0}) + Q_{SH}(T_{0}) \end{aligned}\end{split}\]

where \(\epsilon\) is the emissivity of the surface (snow or ice), \(Q_{S/LW\downarrow}\) the downwelling shortwave and longwave radiation to be prescribed, and \(Q_{LW\uparrow}=\epsilon\sigma_B T_{0}^4\) the emitted long wave radiation with the Stefan-Boltzmann constant \(\sigma_B\). With explicit expressions in \(T_0\) for the turbulent fluxes of latent heat:

\[Q_{LH} = \rho_\mathrm{air} C_E (\Lambda_v + \Lambda_f) |\mathbf{U}_\mathrm{air}|\left[ q_\mathrm{air} - q_\mathrm{sat}(T_0)\right]\]

where \(\rho_\mathrm{air}\) is the air density (parameter SEAICE_rhoAir), \(C_E\) is the ice-ocean 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). Similarly, for sensible heat:

\[Q_{SH}= \rho_\mathrm{air} c_p C_E |\mathbf{U}_\mathrm{air}|\left[ T_\mathrm{10m} - T_{0} \right]\]

where \(c_p\) is the specific heat of air (parameter SEAICE_cpAir). (8.21) can be solved for \(T_0\) with an iterative Ralphson-Newton method which usually converges very quickly, in less that 10 iterations. For the latent heat \(Q_{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 zero-layer 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 sub-grid 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{2n-1}{7}\,h\) for \(n\in[1,N]\). The heat fluxes computed for each thickness category are area-averaged to give the total heat flux (see Hibler 1984 [hibler84]). 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 run-time parameter SEAICE_PDF (1D-array of lenght nITD) can be used to prescribe an arbitrary distribution of ice thicknesses, for example derived from observed distributions (Castro-Morales et al. 2014 [castro-morales14]). In order to include the ice thickness distribution also for snow, set SEAICE_useMultDimSnow = .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 Castro-Morales et al. 2014 [castro-morales14]).

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_{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 (run-time 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 so-called lead closing parameter \(h_{0}\) (run-time 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 [zha98a]). Snow modifies the effective conductivity according to

\[\frac{K}{h} \rightarrow \frac{1}{\frac{h_{s}}{K_{s}}+\frac{h}{K}},\]

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 flood-freeze algorithm following Archimedes’ principle) turns snow into ice until the ice surface is back at \(z=0\) (see Leppäranta 1983 [leppaeranta83]). The flood-freeze algorithm is turned on with run-time parameter SEAICEuseFlooding=.TRUE..

8.6.2.4.12. Advection of thermodynamic variables

Effective ice thickness (ice volume per unit area, \(c\cdot{h}\)), concentration \(c\) and effective snow thickness (\(c\cdot{h}_{s}\)) are advected by ice velocities:

(8.22)\[\frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\mathbf{u}\,X\right) + \Gamma_{X} + D_{X}\]

where \(\Gamma_X\) are the thermodynamic source terms and \(D_{X}\) the diffusive terms for quantities \(X=(c\cdot{h}), c, (c\cdot{h}_{s})\). From the various advection schemes that are available in MITgcm, we recommend flux-limited 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\). Run-time flags: SEAICEadvScheme (default=77, is a 2nd-order 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 3-layer 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 2nd-order advection scheme with a SuperBee flux limiter (Roe 1985 [roe:85]) should be used to advect all sea-ice-related quantities of the Winton (2000) [win00] thermodynamic model (run-time 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 ice-volume and enthalpy) but then false enthalpy extrema can occur, which then leads to unrealistic ice temperature. In the currently implemented solution, the sea-ice mass flux is used to advect the enthalpy in order to ensure conservation of enthalpy and to prevent false enthalpy extrema.

8.6.2.4.13. 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/CICE-Consortium/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) [thorndike:75], and the ice strength parameterization by Rothrock (1975) [rothrock:75] which uses this information. The following description is slightly modified from Ungermann et al. (2017) [ungermann:17]. 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 subgrid-scale (see Thorndike et al. 1975 [thorndike:75]), a probability density function for thickness \(h\) following the evolution equation

(8.23)\[\frac{\partial g}{\partial t} = - \nabla \cdot (\mathbf{u} g) - \frac{\partial}{\partial h}(fg) + \Psi.\]

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 (run-time 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 run-time parameter SEAICEsimpleRidging =.FALSE. in data.seaice). For the ITD model, the ridging mode in convergence

\[\omega_r(h)= \frac{-a(h)+n(h)}{N}\]

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) [thorndike:75] (run-time parameter SEAICEpartFunc =0) or Lipscomb et al. (2007) [lipscomb:07] (SEAICEpartFunc =1), and similarly the ridging function \(n(h)\) can be computed following Hilber (1980) [hib80] (run-time parameter SEAICEredistFunc =0) or Lipscomb et al. (2007) [lipscomb:07] (SEAICEredistFunc =1). As an example, we show here the functions that Lipscomb et al. (2007) [lipscomb:07] suggested to avoid noise in the solutions. These functions are smooth and avoid non-differentiable discontinuities, but so far we did not find any noise issues as in Lipscomb et al. (2007) [lipscomb:07].

With SEAICEpartFunc =1 in data.seaice, the participation function with the relative amount of ice of thickness \(h\) weighted by an exponential function

\[b(h) = b_0 \exp [ -G(h)/a^*]\]

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

\[n(h) = \int_0^\infty a(h_1)\gamma(h_1,h) \textrm{d} h_1\]

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

\[\gamma(h_1, h) = \frac{1}{k \lambda} \exp\left[{\frac{-(h-h_{\min})}{\lambda}}\right]\]

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 e-folding 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) [bitz:01]. For each thickness category in an ITD configuration, the volume conservation equation (8.22) is evaluated using the heat flux with the category-specific 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\) (run-time 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.23) that is discretized in thickness space with \(n+1\) category limits given by run-time parameter Hlimit. If Hlimit is not set in data.seaice, a simple recursive formula following Lipscomb (2001) [lipscomb:01] is used to compute Hlimit:

\[H_\mathrm{limit}(k) = H_\mathrm{limit}(k-1) + \frac{c_1}{n} + \frac{c_1 c_2}{n} [ 1 + \tanh c_3 (\frac{k-1}{n} - 1) ]\]

with \(H_\mathrm{limit}(0)=0\,\text{m}\) and \(H_\mathrm{limit}(n)=999.9\,\text{m}\). The three constants are the run-time 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.4), 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) [rothrock:75]

(8.24)\[P = C_f C_p \int_0^\infty h^2 \omega_r(h) \textrm{d}h\]

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.4) can be used for both ITD and single thickness category models. In contrast to (8.4), which is based on the plausible assumption that thick and compact ice is stronger than thin and loose drifting ice, this parameterization (8.24) 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.4), but in practice, the success is not so clear (Ungermann et al. 2007 [ungermann:17]). Ergo, the default is to use (8.4); set useHibler79IceStrength =.FALSE. in data.seaice to change this behavior.

8.6.2.5. Key subroutines

Top-level 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.6. 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)
---------+----------+----------------+-----------------
 sIceLoad|SM      U1|kg/m^2          |sea-ice loading (in Mass of ice+snow / area unit)
---
SEA ICE STATE:
---
 SIarea  |SM      M1|m^2/m^2         |SEAICE fractional ice-covered area [0 to 1]
 SIheff  |SM      M1|m               |SEAICE effective ice thickness
 SIhsnow |SM      M1|m               |SEAICE effective snow thickness
 SIhsalt |SM      M1|g/m^2           |SEAICE effective salinity
 SIuice  |UU      M1|m/s             |SEAICE zonal ice velocity, >0 from West to East
 SIvice  |VV      M1|m/s             |SEAICE merid. ice velocity, >0 from South to North
---
ATMOSPHERIC STATE AS SEEN BY SEA ICE:
---
 SItices |SM  C   M1|K               |Surface Temperature over Sea-Ice (area weighted)
 SIuwind |UM      U1|m/s             |SEAICE zonal 10-m wind speed, >0 increases uVel
 SIvwind |VM      U1|m/s             |SEAICE meridional 10-m wind speed, >0 increases uVel
 SIsnPrcp|SM      U1|kg/m^2/s        |Snow precip. (+=dw) over Sea-Ice (area weighted)
---
FLUXES ACROSS ICE-OCEAN INTERFACE (ATMOS to OCEAN FOR ICE-FREE REGIONS):
---
 SIfu    |UU      U1|N/m^2           |SEAICE zonal surface wind stress, >0 increases uVel
 SIfv    |VV      U1|N/m^2           |SEAICE merid. surface wind stress, >0 increases vVel
 SIqnet  |SM      U1|W/m^2           |Ocean surface heatflux, turb+rad, >0 decreases theta
 SIqsw   |SM      U1|W/m^2           |Ocean surface shortwave radiat., >0 decreases theta
 SIempmr |SM      U1|kg/m^2/s        |Ocean surface freshwater flux, > 0 increases salt
 SIqneto |SM      U1|W/m^2           |Open Ocean Part of SIqnet, turb+rad, >0 decr theta
 SIqneti |SM      U1|W/m^2           |Ice Covered Part of SIqnet, turb+rad, >0 decr theta
---
FLUXES ACROSS ATMOSPHERE-ICE INTERFACE (ATMOS to OCEAN FOR ICE-FREE REGIONS):
---
 SIatmQnt|SM      U1|W/m^2           |Net atmospheric heat flux, >0 decreases theta
 SIatmFW |SM      U1|kg/m^2/s        |Net freshwater flux from atmosphere & land (+=down)
 SIfwSubl|SM      U1|kg/m^2/s        |Freshwater flux of sublimated ice, >0 decreases ice
---
THERMODYNAMIC DIAGNOSTICS:
---
 SIareaPR|SM      M1|m^2/m^2         |SIarea preceeding ridging process
 SIareaPT|SM      M1|m^2/m^2         |SIarea preceeding thermodynamic growth/melt
 SIheffPT|SM      M1|m               |SIheff preceeeding thermodynamic growth/melt
 SIhsnoPT|SM      M1|m               |SIhsnow preceeeding thermodynamic growth/melt
 SIaQbOCN|SM      M1|m/s             |Potential HEFF rate of change by ocean ice flux
 SIaQbATC|SM      M1|m/s             |Potential HEFF rate of change by atm flux over ice
 SIaQbATO|SM      M1|m/s             |Potential HEFF rate of change by open ocn atm flux
 SIdHbOCN|SM      M1|m/s             |HEFF rate of change by ocean ice flux
 SIdSbATC|SM      M1|m/s             |HSNOW rate of change by atm flux over sea ice
 SIdSbOCN|SM      M1|m/s             |HSNOW rate of change by ocean ice flux
 SIdHbATC|SM      M1|m/s             |HEFF rate of change by atm flux over sea ice
 SIdHbATO|SM      M1|m/s             |HEFF rate of change by open ocn atm flux
 SIdHbFLO|SM      M1|m/s             |HEFF rate of change by flooding snow
 SIdAbATO|SM      M1|m^2/m^2/s       |Potential AREA rate of change by open ocn atm flux
 SIdAbATC|SM      M1|m^2/m^2/s       |Potential AREA rate of change by atm flux over ice
 SIdAbOCN|SM      M1|m^2/m^2/s       |Potential AREA rate of change by ocean ice flux
 SIdA    |SM      M1|m^2/m^2/s       |AREA rate of change (net)
---
DYNAMIC/RHEOLOGY DIAGNOSTICS:
---
 SIpress |SM      M1|N/m             |SEAICE strength (with upper and lower limit)
 SIzeta  |SM      M1|kg/s            |SEAICE nonlinear bulk viscosity
 SIeta   |SM      M1|kg/s            |SEAICE nonlinear shear viscosity
 SIsig1  |SM      M1|no units        |SEAICE normalized principle stress, component one
 SIsig2  |SM      M1|no units        |SEAICE normalized principle stress, component two
 SIshear |SM      M1|1/s             |SEAICE shear deformation rate
 SIdelta |SM      M1|1/s             |SEAICE Delta deformation rate
 SItensil|SM      M1|N/m             |SEAICE maximal tensile strength
---
ADVECTIVE/DIFFUSIVE FLUXES OF SEA ICE variables:
---
 ADVxHEFF|UU      M1|m.m^2/s         |Zonal      Advective Flux of eff ice thickn
 ADVyHEFF|VV      M1|m.m^2/s         |Meridional Advective Flux of eff ice thickn
 SIuheff |UU      M1|m^2/s           |Zonal      Transport of eff ice thickn (centered)
 SIvheff |VV      M1|m^2/s           |Meridional Transport of eff ice thickn (centered)
 DFxEHEFF|UU      M1|m^2/s           |Zonal      Diffusive Flux of eff ice thickn
 DFyEHEFF|VV      M1|m^2/s           |Meridional Diffusive Flux of eff ice thickn
 ADVxAREA|UU      M1|m^2/m^2.m^2/s   |Zonal      Advective Flux of fract area
 ADVyAREA|VV      M1|m^2/m^2.m^2/s   |Meridional Advective Flux of fract area
 DFxEAREA|UU      M1|m^2/m^2.m^2/s   |Zonal      Diffusive Flux of fract area
 DFyEAREA|VV      M1|m^2/m^2.m^2/s   |Meridional Diffusive Flux of fract area
 ADVxSNOW|UU      M1|m.m^2/s         |Zonal      Advective Flux of eff snow thickn
 ADVySNOW|VV      M1|m.m^2/s         |Meridional Advective Flux of eff snow thickn
 DFxESNOW|UU      M1|m.m^2/s         |Zonal      Diffusive Flux of eff snow thickn
 DFyESNOW|VV      M1|m.m^2/s         |Meridional Diffusive Flux of eff snow thickn
 ADVxSSLT|UU      M1|psu.m^2/s       |Zonal      Advective Flux of seaice salinity
 ADVySSLT|VV      M1|psu.m^2/s       |Meridional Advective Flux of seaice salinity
 DFxESSLT|UU      M1|psu.m^2/s       |Zonal      Diffusive Flux of seaice salinity
 DFyESSLT|VV      M1|psu.m^2/s       |Meridional Diffusive Flux of seaice salinity

8.6.2.7. Experiments and tutorials that use seaice