Warning, /doc/phys_pkgs/seaice.rst is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 204068be on 2025-11-16 04:31:18 UTC
adc83e5d7b Mart*0001 .. _sub_phys_pkg_seaice:
0002
0003 SEAICE Package
258fe29c91 Jeff*0004 **************
adc83e5d7b Mart*0005
0006 Authors: Martin Losch, Dimitris Menemenlis, An Nguyen, Jean-Michel
c512e371cc drin*0007 Campin, Patrick Heimbach, Chris Hill, Jinlun Zhang, and Damien Ringeisen
adc83e5d7b Mart*0008
0009 .. _ssub_phys_pkg_seaice_intro:
0010
0011 Introduction
258fe29c91 Jeff*0012 ============
adc83e5d7b Mart*0013
c512e371cc drin*0014 Package :filelink:`seaice <pkg/seaice>` provides a dynamic and thermodynamic
0015 interactive sea ice model.
adc83e5d7b Mart*0016
0017 CPP options enable or disable different aspects of the package
9986b4a53e Jeff*0018 (:numref:`ssub_phys_pkg_seaice_config`). Run-time options, flags, filenames and
c512e371cc drin*0019 field-related dates/times are set in ``data.seaice``
0020 (:numref:`ssub_phys_pkg_seaice_runtime`). A description of key subroutines is
0021 given in :numref:`ssub_phys_pkg_seaice_subroutines`. Available diagnostics
0022 output is listed in :numref:`ssub_phys_pkg_seaice_diagnostics`.
adc83e5d7b Mart*0023
61f2157921 Oliv*0024 .. _ssub_phys_pkg_seaice_config:
adc83e5d7b Mart*0025
61f2157921 Oliv*0026 SEAICE configuration and compiling
258fe29c91 Jeff*0027 ==================================
adc83e5d7b Mart*0028
61f2157921 Oliv*0029 Compile-time options
258fe29c91 Jeff*0030 --------------------
adc83e5d7b Mart*0031
c512e371cc drin*0032 As with all MITgcm packages, SEAICE can be turned on or off at compile time
0033 (see :numref:`building_code`)
adc83e5d7b Mart*0034
c512e371cc drin*0035 - using the ``packages.conf`` file by adding ``seaice`` to it
adc83e5d7b Mart*0036
c512e371cc drin*0037 - or using :filelink:`genmake2 <tools/genmake2>` adding ``-enable=seaice`` or
0038 ``-disable=seaice`` switches
adc83e5d7b Mart*0039
c512e371cc drin*0040 - **required packages and CPP options**:
0041 :filelink:`seaice <pkg/seaice>` requires the external forcing package
0042 :filelink:`pkg/exf` to be enabled; no additional CPP options are required.
adc83e5d7b Mart*0043
0044
c512e371cc drin*0045 Parts of the :filelink:`seaice <pkg/seaice>` code can be enabled or disabled at
0046 compile time via CPP preprocessor flags. These options are set in
0047 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>`.
258fe29c91 Jeff*0048 :numref:`tab_phys_pkg_seaice_cpp` summarizes the most important ones. For more
382462ccb5 Mart*0049 options see :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>`. Note
0050 that defining :varlink:`SEAICE_BGRID_DYNAMICS` turns on legacy code and thus
0051 automatically undefines more recent features, see :filelink:`SEAICE_OPTIONS.h
0052 <pkg/seaice/SEAICE_OPTIONS.h>` for details.
258fe29c91 Jeff*0053
0054 .. tabularcolumns:: |\Y{.375}|\Y{.1}|\Y{.55}|
adc83e5d7b Mart*0055
c512e371cc drin*0056 .. csv-table:: Some of the most relevant CPP preprocessor flags in the :filelink:`seaice <pkg/seaice>` package.
258fe29c91 Jeff*0057 :header: "CPP option", "Default", Description"
0058 :widths: 30, 10, 60
bc3b9fecef Mart*0059 :name: tab_phys_pkg_seaice_cpp
0060
258fe29c91 Jeff*0061 :varlink:`SEAICE_DEBUG`, #undef, enhance STDOUT for debugging
382462ccb5 Mart*0062 :varlink:`SEAICE_CGRID`, #define, use sea ice dynamics on C-grid
258fe29c91 Jeff*0063 :varlink:`SEAICE_ALLOW_EVP`, #define, enable use of EVP rheology solver
0064 :varlink:`SEAICE_ALLOW_JFNK`, #define, enable use of JFNK rheology solver
0065 :varlink:`SEAICE_ALLOW_KRYLOV`, #define, enable use of Krylov rheology solver
c512e371cc drin*0066 :varlink:`SEAICE_ALLOW_TEM`, #undef, enable use of the truncated ellipse method (TEM) and coulombic yield curve
0067 :varlink:`SEAICE_ALLOW_MCS`, #undef, enable use of Mohr-Coulomb yield curve with shear flow rule
0068 :varlink:`SEAICE_ALLOW_MCE`, #undef, enable use of Mohr-Coulomb yield curve with elliptical plastic potential
0069 :varlink:`SEAICE_ALLOW_TD`, #undef, enable use of teardrop and parabolic Lens yield curves with normal flow rules
258fe29c91 Jeff*0070 :varlink:`SEAICE_LSR_ZEBRA`, #undef, use a coloring method for LSR solver
a4e168e012 antn*0071 :varlink:`SEAICE_ALLOW_FREEDRIFT`, #undef, enable solve approximate sea ice momentum equation and bypass solving for sea ice internal stress
258fe29c91 Jeff*0072 :varlink:`SEAICE_EXTERNAL_FLUXES`, #define, use :filelink:`pkg/exf`-computed fluxes as starting point
0073 :varlink:`SEAICE_ZETA_SMOOTHREG`, #define, use differentiable regularization for viscosities
14673ec2d0 Mart*0074 :varlink:`SEAICE_DELTA_SMOOTHREG`, #undef, use differentiable regularization :math:`\Delta_{\mathrm{reg}}=\sqrt{\Delta^2+\Delta_{\min}}` instead of :math:`\max`-function for :math:`1/\Delta_{\mathrm{reg}}`
258fe29c91 Jeff*0075 :varlink:`SEAICE_ALLOW_BOTTOMDRAG`, #undef, enable grounding parameterization for improved fastice in shallow seas
5bb179ddc2 Mart*0076 :varlink:`SEAICE_ALLOW_SIDEDRAG`, #undef, enable lateral drag parameterization for improved fastice along coastlines and islands
382462ccb5 Mart*0077 :varlink:`SEAICE_BGRID_DYNAMICS`, #undef, use sea ice dynamics code on legacy B-grid; most of the previous flags are not available with B-grid
0078 :varlink:`SEAICE_BICE_STRESS`, #undef, B-grid only for backward compatiblity: turn on ice-stress on ocean; defined by default if :varlink:`SEAICE_BGRID_DYNAMICS` is defined
0079 :varlink:`EXPLICIT_SSH_SLOPE`, #undef, B-grid only for backward compatiblity: use ETAN for tilt computations rather than geostrophic velocities; defined by default if :varlink:`SEAICE_BGRID_DYNAMICS` is defined
0080 :varlink:`SEAICE_LSRBNEW`, #undef, FV discretization for B-grid
258fe29c91 Jeff*0081 :varlink:`SEAICE_ITD`, #undef, run with dynamical sea Ice Thickness Distribution (ITD)
0082 :varlink:`SEAICE_VARIABLE_SALINITY`, #undef, enable sea ice with variable salinity
e2fbc60f23 Jeff*0083 :varlink:`SEAICE_CAP_ICELOAD`, #undef, enable to limit seaice load (:varlink:`siceLoad`) on the sea surface
258fe29c91 Jeff*0084 :varlink:`ALLOW_SITRACER`, #undef, enable sea ice tracer package
a4e168e012 antn*0085 :varlink:`SEAICE_USE_GROWTH_ADX`, #undef, use of adjointable but more simplified sea ice thermodynamics model in :filelink:`seaice_growth_adx.F <pkg/seaice/seaice_growth_adx.F>` instead of :filelink:`seaice_growth.F <pkg/seaice/seaice_growth.F>`
adc83e5d7b Mart*0086
61f2157921 Oliv*0087 .. _ssub_phys_pkg_seaice_runtime:
adc83e5d7b Mart*0088
2c231b0ebd Mart*0089 Run-time parameters
9986b4a53e Jeff*0090 ===================
adc83e5d7b Mart*0091
9986b4a53e Jeff*0092 Run-time parameters (see :numref:`tab_phys_pkg_seaice_runtimeparms`) are set in
0bad585a21 Navi*0093 ``data.seaice`` (read in :filelink:`pkg/seaice/seaice_readparms.F`).
adc83e5d7b Mart*0094
0095 Enabling the package
258fe29c91 Jeff*0096 --------------------
adc83e5d7b Mart*0097
c512e371cc drin*0098 :filelink:`seaice <pkg/seaice>` package is switched on/off at run-time by
dc26f158aa Mart*0099 setting :varlink:`useSEAICE` ``= .TRUE.,`` in ``data.pkg``.
adc83e5d7b Mart*0100
0101 General flags and parameters
258fe29c91 Jeff*0102 ----------------------------
adc83e5d7b Mart*0103
9986b4a53e Jeff*0104 :numref:`tab_phys_pkg_seaice_runtimeparms` lists most run-time parameters.
adc83e5d7b Mart*0105
258fe29c91 Jeff*0106 .. tabularcolumns:: |\Y{.275}|\Y{.20}|\Y{.525}|
adc83e5d7b Mart*0107
9986b4a53e Jeff*0108 .. table:: Run-time parameters and default values
258fe29c91 Jeff*0109 :class: longtable
adc83e5d7b Mart*0110 :name: tab_phys_pkg_seaice_runtimeparms
0111
258fe29c91 Jeff*0112 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0113 | Name | Default value | Description |
0114 +====================================+==============================+=========================================================================+
c61841e2fd Jeff*0115 | :varlink:`SEAICEwriteState` | FALSE | write sea ice state to file |
258fe29c91 Jeff*0116 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0117 | :varlink:`SEAICEuseDYNAMICS` | TRUE | use dynamics |
0118 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0119 | :varlink:`SEAICEuseJFNK` | FALSE | use the JFNK-solver |
0120 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c512e371cc drin*0121 | :varlink:`SEAICEuseTEM` | FALSE | use truncated ellipse method or coulombic yield curve |
0122 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0123 | :varlink:`SEAICEuseMCS` | FALSE | use the Mohr-Coulomb yield curve with shear flow rule |
0124 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0125 | :varlink:`SEAICEuseMCE` | FALSE | use the Mohr-Coulomb yield curve with elliptical plastic potential |
0126 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0127 | :varlink:`SEAICEuseTD` | FALSE | use the teardrop yield curve with normal flow rule |
0128 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0129 | :varlink:`SEAICEusePL` | FALSE | use the parabolic Lens yield curve with normal flow rule |
258fe29c91 Jeff*0130 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0131 | :varlink:`SEAICEuseStrImpCpl` | FALSE | use strength implicit coupling in LSR/JFNK |
0132 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0133 | :varlink:`SEAICEuseMetricTerms` | TRUE | use metric terms in dynamics |
0134 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0135 | :varlink:`SEAICEuseEVPpickup` | TRUE | use EVP pickups |
0136 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
a4e168e012 antn*0137 | :varlink:`SEAICEuseFREEDRIFT` | FALSE | solve approximate momentum equation, bypassing rheology |
0138 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c61841e2fd Jeff*0139 | :varlink:`SEAICEuseFluxForm` | TRUE | use flux form for 2nd central difference advection scheme |
258fe29c91 Jeff*0140 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0141 | :varlink:`SEAICErestoreUnderIce` | FALSE | enable restoring to climatology under ice |
0142 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0143 | :varlink:`SEAICEupdateOceanStress` | TRUE | update ocean surface stress accounting for sea ice cover |
0144 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
dc26f158aa Mart*0145 | :varlink:`SEAICEscaleSurfStress` | TRUE | scale atmosphere and ocean-surface stress on ice by concentration (AREA)|
258fe29c91 Jeff*0146 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0147 | :varlink:`SEAICEaddSnowMass` | TRUE | in computing seaiceMass, add snow contribution |
0148 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0149 | :varlink:`useHB87stressCoupling` | FALSE | turn on ice-ocean stress coupling following |
0150 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0151 | :varlink:`usePW79thermodynamics` | TRUE | flag to turn off zero-layer-thermodynamics for testing |
0152 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0153 | :varlink:`SEAICEadvHeff` | TRUE | flag to turn off advection of scalar variable :varlink:`HEFF` |
0154 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0155 | :varlink:`SEAICEadvArea` | TRUE | flag to turn off advection of scalar variable :varlink:`AREA` |
0156 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0157 | :varlink:`SEAICEadvSnow` | TRUE | flag to turn off advection of scalar variable :varlink:`HSNOW` |
0158 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0159 | :varlink:`SEAICEadvSalt` | TRUE | flag to turn off advection of scalar variable :varlink:`HSALT` |
0160 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0161 | :varlink:`SEAICEadvScheme` | 77 | set advection scheme for seaice scalar state variables |
0162 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0163 | :varlink:`SEAICEuseFlooding` | TRUE | use flood-freeze algorithm |
0164 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
a4e168e012 antn*0165 | :varlink:`SINegFac` | 1.0 | over/undershoot factor for seaice advective term in forward/adjoint |
0166 | | | (SEAICE_USE_GROWTH_ADX only) |
0167 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0168 | :varlink:`SEAICE_no_slip` | FALSE | use no-slip boundary conditions instead of free-slip |
0169 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0170 | :varlink:`SEAICE_deltaTtherm` | :varlink:`dTtracerLev` (1) | time step for seaice thermodynamics (s) |
0171 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0172 | :varlink:`SEAICE_deltaTdyn` | :varlink:`dTtracerLev` (1) | time step for seaice dynamics (s) |
0173 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0174 | :varlink:`SEAICE_deltaTevp` | 0.0 | EVP sub-cycling time step (s); values :math:`>` 0 turn on EVP |
0175 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
dc26f158aa Mart*0176 | :varlink:`SEAICEuseEVPstar` | TRUE | use modified EVP\* instead of EVP, following :cite:`lemieux:12` |
258fe29c91 Jeff*0177 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
dc26f158aa Mart*0178 | :varlink:`SEAICEuseEVPrev` | TRUE | "revisited form" variation on EVP\*, following :cite:`bouillon:13` |
258fe29c91 Jeff*0179 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0180 | :varlink:`SEAICEnEVPstarSteps` | unset | number of modified EVP\* iterations |
0181 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0182 | :varlink:`SEAICE_evpAlpha` | unset | EVP\* parameter (non-dim.), to replace |
0183 | | | 2*\ :varlink:`SEAICE_evpTauRelax`\ /\ :varlink:`SEAICE_deltaTevp` |
0184 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0185 | :varlink:`SEAICE_evpBeta` | unset | EVP\* parameter (non-dim.), to replace |
0186 | | | :varlink:`SEAICE_deltaTdyn`\ /\ :varlink:`SEAICE_deltaTevp` |
0187 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0188 | :varlink:`SEAICEaEVPcoeff` | unset | largest stabilized frequency for adaptive EVP (non-dim.) |
0189 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c512e371cc drin*0190 | :varlink:`SEAICEaEVPcStar` | 4.0 | aEVP multiple of stability factor (non-dim.), see :cite:`kimmritz:16` |
258fe29c91 Jeff*0191 | | | :math:`\alpha * \beta = c^\ast * \gamma` |
0192 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c512e371cc drin*0193 | :varlink:`SEAICEaEVPalphaMin` | 5.0 | aEVP lower limit of alpha and beta (non-dim.), see :cite:`kimmritz:16` |
258fe29c91 Jeff*0194 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
204068be6a R. S*0195 | :varlink:`SEAICE_evpAreaReg` | -1.0 | a minimun ice fraction for regularizations of denomU/V in EVP; |
0196 | | | off by default (-1), 1.E-5 is a useful value for high-res simulation |
0197 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0198 | :varlink:`SEAICE_elasticParm` | 0.33333333 | EVP parameter :math:`E_0` (non-dim.), sets relaxation timescale |
c512e371cc drin*0199 | | | :varlink:`SEAICE_evpTauRelax` = |
258fe29c91 Jeff*0200 | | | :varlink:`SEAICE_elasticParm` * :varlink:`SEAICE_deltaTdyn` |
0201 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0202 | :varlink:`SEAICE_evpTauRelax` | :varlink:`dTtracerLev` (1) * | relaxation time scale :math:`T` for EVP waves (s) |
0203 | | :varlink:`SEAICE_elasticParm`| |
0204 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0205 | :varlink:`SEAICE_OLx` | :varlink:`OLx` - 2 | overlap for LSR-solver or preconditioner, :math:`x`-dimension |
0206 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0207 | :varlink:`SEAICE_OLy` | :varlink:`OLy` - 2 | overlap for LSR-solver or preconditioner, :math:`y`-dimension |
0208 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c61841e2fd Jeff*0209 | :varlink:`SEAICEnonLinIterMax` | 2/10 | maximum number of non-linear (outer loop) iterations (LSR/JFNK) |
258fe29c91 Jeff*0210 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c61841e2fd Jeff*0211 | :varlink:`SEAICElinearIterMax` | 1500/10 | maximum number of linear iterations (LSR/JFNK) |
258fe29c91 Jeff*0212 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0213 | :varlink:`SEAICE_JFNK_lsIter` | (off) | start line search after “lsIter†Newton iterations |
0214 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c704c5a1ef Mart*0215 | :varlink:`SEAICE_JFNK_lsLmax` | 4 | maximum number of line search steps |
0216 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0217 | :varlink:`SEAICE_JFNK_lsGamma` | 0.5 | line search step size parameter |
0218 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0219 | :varlink:`SEAICEnonLinTol` | 1.0E-05 | non-linear tolerance parameter for JFNK solver |
0220 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0221 | :varlink:`JFNKgamma_lin_min` | 0.10 | minimum tolerance parameter for linear JFNK solver |
0222 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0223 | :varlink:`JFNKgamma_lin_max` | 0.99 | maximum tolerance parameter for linear JFNK solver |
0224 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0225 | :varlink:`JFNKres_tFac` | unset | tolerance parameter for FGMRES residual |
0226 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0227 | :varlink:`SEAICE_JFNKepsilon` | 1.0E-06 | step size for the FD-gradient in s/r seaice_jacvec |
0228 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0229 | :varlink:`SEAICE_dumpFreq` | dumpFreq | dump frequency (s) |
0230 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0231 | :varlink:`SEAICE_dump_mdsio` | TRUE | write snapshot using :filelink:`/pkg/mdsio` |
0232 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0233 | :varlink:`SEAICE_dump_mnc` | FALSE | write snapshot using :filelink:`/pkg/mnc` |
0234 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0235 | :varlink:`SEAICE_initialHEFF` | 0.0 | initial sea ice thickness averaged over grid cell (m) |
0236 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0237 | :varlink:`SEAICE_drag` | 1.0E-03 | air-ice drag coefficient (non-dim.) |
0238 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0239 | :varlink:`OCEAN_drag` | 1.0E-03 | air-ocean drag coefficient (non-dim.) |
0240 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0241 | :varlink:`SEAICE_waterDrag` | 5.5E-03 | water-ice drag coefficient (non-dim.) |
0242 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0243 | :varlink:`SEAICE_dryIceAlb` | 0.75 | winter sea ice albedo |
0244 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0245 | :varlink:`SEAICE_wetIceAlb` | 0.66 | summer sea ice albedo |
0246 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0247 | :varlink:`SEAICE_drySnowAlb` | 0.84 | dry snow albedo |
0248 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0249 | :varlink:`SEAICE_wetSnowAlb` | 0.70 | wet snow albedo |
0250 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0251 | :varlink:`SEAICE_waterAlbedo` | 0.10 | water albedo (not used if #define :varlink:`SEAICE_EXTERNAL_FLUXES`) |
0252 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0253 | :varlink:`SEAICE_strength` | 2.75E+04 | sea ice strength constant :math:`P^{\ast}` (N/m\ :sup:`2`) |
0254 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0255 | :varlink:`SEAICE_cStar` | 20.0 | sea ice strength constant :math:`C^{\ast}` (non-dim.) |
0256 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
aa30c76f3a Dami*0257 | :varlink:`SEAICE_eccen` | 2.0 | VP rheology ellipse aspect ratio :math:`e` |
0258 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c512e371cc drin*0259 | :varlink:`SEAICE_eccfr` | = :varlink:`SEAICE_eccen` | sea ice plastic potential ellipse aspect ratio :math:`e_G` |
0260 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0261 | :varlink:`SEAICEmcMU` | 1.0 | slope of the Mohr-Coulomb yield curve |
0262 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0263 | :varlink:`SEAICEpressReplFac` | 1.0 | use replacement pressure (0.0-1.0) |
0264 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0265 | :varlink:`SEAICE_tensilFac` | 0.0 | tensile factor for the yield curve |
0266 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0267 | :varlink:`SEAICE_rhoAir` | 1.3 (or | density of air (kg/m\ :sup:`3`) |
0268 | | :filelink:`pkg/exf` value) | |
0269 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0270 | :varlink:`SEAICE_cpAir` | 1004.0 (or | specific heat of air (J/kg/K) |
0271 | | :filelink:`pkg/exf` value) | |
0272 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0273 | :varlink:`SEAICE_lhEvap` | 2.5E+06 (or | latent heat of evaporation (J/kg) |
0274 | | :filelink:`pkg/exf` value) | |
0275 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0276 | :varlink:`SEAICE_lhFusion` | 3.34E+05 (or | latent heat of fusion (J/kg) |
0277 | | :filelink:`pkg/exf` value) | |
0278 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c61841e2fd Jeff*0279 | :varlink:`SEAICE_dalton` | 1.75E-03 | ice-ocean transfer coefficient for latent and sensible heat (non-dim.) |
258fe29c91 Jeff*0280 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
de1b16b92a Jeff*0281 | :varlink:`useMaykutSatVapPoly` | FALSE | use Maykut polynomial to compute saturation vapor pressure |
0282 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0283 | :varlink:`SEAICE_iceConduct` | 2.16560E+00 | sea ice conductivity (W m\ :sup:`-1` K\ :sup:`-1`) |
0284 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0285 | :varlink:`SEAICE_snowConduct` | 3.10000E-01 | snow conductivity (W m\ :sup:`-1` K\ :sup:`-1`) |
0286 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0287 | :varlink:`SEAICE_emissivity` | 0.970018 (or | longwave ocean surface emissivity (non-dim.) |
0288 | | :filelink:`pkg/exf` value) | |
0289 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0290 | :varlink:`SEAICE_snowThick` | 0.15 | cutoff snow thickness to use snow albedo (m) |
0291 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0292 | :varlink:`SEAICE_shortwave` | 0.30 | ice penetration shortwave radiation factor (non-dim.) |
0293 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0294 | :varlink:`SEAICE_saltFrac` | 0.0 | salinity newly formed ice (as fraction of ocean surface salinity) |
0295 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0296 | :varlink:`SEAICE_frazilFrac` | 1.0 (or | frazil to sea ice conversion rate, as fraction |
0297 | | computed from other parms) | (relative to the local freezing point of sea ice water) |
0298 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0299 | :varlink:`SEAICEstressFactor` | 1.0 | scaling factor for ice area in computing total ocean stress (non-dim.) |
0300 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0301 | :varlink:`HeffFile` | unset | filename for initial sea ice eff. thickness field :varlink:`HEFF` (m) |
0302 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0303 | :varlink:`AreaFile` | unset | filename for initial fraction sea ice cover :varlink:`AREA` (non-dim.) |
0304 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0305 | :varlink:`HsnowFile` | unset | filename for initial eff. snow thickness field :varlink:`HSNOW` (m) |
0306 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0307 | :varlink:`HsaltFile` | unset | filename for initial eff. sea ice salinity field :varlink:`HSALT` |
0308 | | | (g/m\ :sup:`2`) |
0309 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
dc26f158aa Mart*0310 | :varlink:`LSR_ERROR` | 1.0E-05 | sets accuracy of LSR solver |
258fe29c91 Jeff*0311 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0312 | :varlink:`DIFF1` | 0.0 | parameter used in advect.F |
0313 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0314 | :varlink:`HO` | 0.5 | lead closing parameter :math:`h_0` (m); demarcation thickness between |
0315 | | | thick and thin ice which determines partition between vertical and |
0316 | | | lateral ice growth |
0317 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0318 | :varlink:`MIN_ATEMP` | -50.0 | minimum air temperature (:sup:`o`\ C) |
0319 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0320 | :varlink:`MIN_LWDOWN` | 60.0 | minimum downward longwave (W/m\ :sup:`2`) |
0321 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0322 | :varlink:`MIN_TICE` | -50.0 | minimum ice temperature (:sup:`o`\ C) |
0323 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0324 | :varlink:`IMAX_TICE` | 10 | number of iterations for ice surface temperature solution |
0325 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0326 | :varlink:`SEAICE_EPS` | 1.0E-10 | a "small number" used in various routines |
2c231b0ebd Mart*0327 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
14673ec2d0 Mart*0328 | :varlink:`SEAICE_deltaMin` | :varlink:`SEAICE_EPS` | minimum to regularize :math:`\Delta` |
0329 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0330 | :varlink:`SEAICE_area_reg` | 1.0E-5 | minimum concentration to regularize ice thickness |
0331 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0332 | :varlink:`SEAICE_hice_reg` | 0.05 | minimum ice thickness (m) for regularization |
0333 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0334 | :varlink:`SEAICE_multDim` | 1 | number of ice categories for thermodynamics |
0335 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0336 | :varlink:`SEAICE_useMultDimSnow` | TRUE | use same fixed pdf for snow as for multi-thickness-category ice |
0337 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
5bb179ddc2 Mart*0338 | :varlink:`SEAICEbasalDragK1` | 8.0 | basal drag parameter K\ :sub:`1` :cite:`lemieux:15` |
0339 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0340 | :varlink:`SEAICEbasalDragK2` | 0.0 | basal drag parameter K\ :sub:`2` |
0341 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0342 | :varlink:`SEAICE_cBasalStar` | :varlink:`SEAICE_cStar` value| basal drag parameter (no units) |
0343 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0344 | :varlink:`SEAICEbasalDragU0` | 5.E-5 | basal drag parameter (m/s) |
0345 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0346 | :varlink:`SEAICESideDrag` | 0.0 | lateral drag coefficient :cite:`liu:22` |
0347 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0348 | :varlink:`uCoastLineFile` | unset | filename for coastline length for u-equation |
0349 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0350 | :varlink:`vCoastLineFile` | unset | filename for coastline length for v-equation |
0351 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
258fe29c91 Jeff*0352
0353
c512e371cc drin*0354 The following dynamical ice thickness distribution and ridging parameters in
0355 :numref:`tab_phys_pkg_seaice_ridging` are only active with #define
0356 :varlink:`SEAICE_ITD`. All parameters are non-dimensional unless indicated.
258fe29c91 Jeff*0357
0358 .. tabularcolumns:: |\Y{.275}|\Y{.20}|\Y{.525}|
0359
0360 .. table:: Thickness distribution and ridging parameters
0361 :name: tab_phys_pkg_seaice_ridging
0362
0363
0364 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0365 | Name | Default value | Description |
0366 +====================================+==============================+=========================================================================+
c512e371cc drin*0367 | :varlink:`useHibler79IceStrength` | TRUE | use :cite:`hibler:79` ice strength; do not use :cite:`rothrock:75` |
258fe29c91 Jeff*0368 | | | with #define :varlink:`SEAICE_ITD` |
0369 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c512e371cc drin*0370 | :varlink:`SEAICEsimpleRidging` | TRUE | use simple ridging a la :cite:`hibler:79` |
258fe29c91 Jeff*0371 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0372 | :varlink:`SEAICE_cf` | 17.0 | scaling parameter of :cite:`rothrock:75` ice strength parameterization |
0373 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0374 | :varlink:`SEAICEpartFunc` | 0 | use partition function of :cite:`thorndike:75` |
0375 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c512e371cc drin*0376 | :varlink:`SEAICEredistFunc` | 0 | use redistribution function of :cite:`hibler:80` |
258fe29c91 Jeff*0377 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0378 | :varlink:`SEAICEridgingIterMax` | 10 | maximum number of ridging sweeps |
0379 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0380 | :varlink:`SEAICEshearParm` | 0.5 | fraction of shear to be used for ridging |
0381 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0382 | :varlink:`SEAICEgStar` | 0.15 | max. ice conc. that participates in ridging :cite:`thorndike:75` |
0383 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0384 | :varlink:`SEAICEhStar` | 25.0 | ridging parameter for :cite:`thorndike:75`, :cite:`lipscomb:07` |
0385 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0386 | :varlink:`SEAICEaStar` | 0.05 | similar to :varlink:`SEAICEgStar` for |
0387 | | | :cite:`lipscomb:07` participation function |
0388 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0389 | :varlink:`SEAICEmuRidging` | 3.0 | similar to :varlink:`SEAICEhStar` for |
0390 | | | :cite:`lipscomb:07` ridging function |
0391 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
2c231b0ebd Mart*0392 | :varlink:`SEAICEmaxRaft` | 1.0 | regularization parameter for rafting |
258fe29c91 Jeff*0393 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0394 | :varlink:`SEAICEsnowFracRidge` | 0.5 | fraction of snow that remains on ridged ice |
0395 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0396 | :varlink:`SEAICEuseLinRemapITD` | TRUE | use linear remapping scheme of :cite:`lipscomb:01` |
0397 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
c61841e2fd Jeff*0398 | :varlink:`Hlimit` | unset | nITD+1-array of ice thickness category limits (m) |
258fe29c91 Jeff*0399 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
0400 | :varlink:`Hlimit_c1`, | 3.0, | when :varlink:`Hlimit` is not set, then these parameters |
0401 | :varlink:`Hlimit_c2`, | 15.0, | determine :varlink:`Hlimit` from a simple function |
0402 | :varlink:`Hlimit_c3` | 3.0 | following :cite:`lipscomb:01` |
0403 +------------------------------------+------------------------------+-------------------------------------------------------------------------+
2c231b0ebd Mart*0404
adc83e5d7b Mart*0405
0406 .. _ssub_phys_pkg_seaice_descr:
0407
0408 Description
258fe29c91 Jeff*0409 ===========
adc83e5d7b Mart*0410
c512e371cc drin*0411 The MITgcm sea ice model is based on a variant of the viscous-plastic (VP)
0412 dynamic-thermodynamic sea ice model (Zhang and Hibler 1997 :cite:`zhang:97`)
0413 first introduced in Hibler (1979) and Hibler (1980)
0414 :cite:`hibler:79,hibler:80`. In order to adapt this model to the requirements
0415 of coupled ice-ocean state estimation, many important aspects of the original
0416 code have been modified and improved, see Losch et al. (2010) :cite:`losch:10`:
adc83e5d7b Mart*0417
0418 - the code has been rewritten for an Arakawa C-grid, both B- and C-grid
c512e371cc drin*0419 variants are available; the C-grid code allows for no-slip and free-slip
0420 lateral boundary conditions;
adc83e5d7b Mart*0421
0422 - three different solution methods for solving the nonlinear momentum
c512e371cc drin*0423 equations have been adopted: LSOR (Zhang and Hibler 1997 :cite:`zhang:97`),
0424 EVP (Hunke and Dukowicz 1997 :cite:`hunke:97`),
0425 JFNK (Lemieux et al. 2010 :cite:`lemieux:10`, Losch et al. 2014
0426 :cite:`losch:14`);
adc83e5d7b Mart*0427
258fe29c91 Jeff*0428 - ice-ocean stress can be formulated as in Hibler and Bryan (1987)
c512e371cc drin*0429 :cite:`hibler:87` or as in Campin et al. (2008) :cite:`campin:08`;
adc83e5d7b Mart*0430
0431 - ice variables are advected by sophisticated, conservative advection
0432 schemes with flux limiting;
0433
0434 - growth and melt parameterizations have been refined and extended in
0435 order to allow for more stable automatic differentiation of the code.
0436
0437 The sea ice model is tightly coupled to the ocean compontent of the
c512e371cc drin*0438 MITgcm. Heat, fresh water fluxes and surface stresses are computed from the
0439 atmospheric state and, by default, modified by the ice model at every time
0440 step.
0441
0442 The ice dynamics models that are most widely used for large-scale climate
0443 studies are the viscous-plastic (VP) model (Hilber 1979 :cite:`hibler:79`), the
0444 cavitating fluid (CF) model (Flato and Hibler 1992 :cite:`flato:92`), and the
0445 elastic-viscous-plastic (EVP) model (Hunke and Dukowicz 1997 :cite:`hunke:97`).
0446 Compared to the VP model, the CF model does not allow ice shear in calculating
0447 ice motion, stress, and deformation. EVP models approximate VP by adding an
0448 elastic term to the equations for easier adaptation to parallel
0449 computers. Because of its higher accuracy in plastic solution and relatively
0450 simpler formulation, compared to the EVP model, we decided to use the VP model
0451 as the default dynamic component of our ice model. To do this we extended the
0452 line successive over relaxation (LSOR) method of Zhang and Hibler (1997)
0453 :cite:`zhang:97` for use in a parallel configuration. An EVP model and a
9986b4a53e Jeff*0454 free-drift implementation can be selected with run-time flags.
adc83e5d7b Mart*0455
a4e168e012 antn*0456 :filelink:`pkg/seaice` includes the original so-called zero-layer
0457 thermodynamics with a snow cover as in the appendix of Semtner (1976)
0458 :cite:`semtner:76`. Two versions of this zero-layer thermodynamic code exist,
0459 with a more developed version :filelink:`seaice_growth.F
0460 <pkg/seaice/seaice_growth.F>` and a simplified version
0461 :filelink:`seaice_growth_adx.F <pkg/seaice/seaice_growth_adx.F>` based on
0462 Fenty (2013) :cite:`fenty:13` that excludes physics such as ITD, treatment for
0463 sublimation, and frazil ice but provides a stable sea ice adjointable with
0464 physical sensitivity. When the seaice_growth_adx code is enabled (by defining
0465 :varlink:`SEAICE_USE_GROWTH_ADX` in :filelink:`SEAICE_OPTIONS.h
0466 <pkg/seaice/SEAICE_OPTIONS.h>`), the regularization parameter
0467 :varlink:`SINegFac` is set to zero in adjoint mode to disable the potential
0468 propagation of unphysical terms associated with sea ice dynamics.
0469
adc83e5d7b Mart*0470
0471 .. _para_phys_pkg_seaice_thsice:
0472
258fe29c91 Jeff*0473 Compatibility with ice-thermodynamics package :filelink:`pkg/thsice`
0474 --------------------------------------------------------------------
adc83e5d7b Mart*0475
a4e168e012 antn*0476 The zero-layer thermodynamic model assumes that ice does
c512e371cc drin*0477 not store heat and, therefore, tends to exaggerate the seasonal variability in
0478 ice thickness. This exaggeration can be significantly reduced by using Winton's
0479 (Winton 2000 :cite:`winton:00`) three-layer thermodynamic model that permits
0480 heat storage in ice.
0481
0482 The Winton (2000) sea-ice thermodynamics have been ported to MITgcm; they
0483 currently reside under :filelink:`pkg/thsice`, described in
0484 :numref:`sub_phys_pkg_thsice`. It is fully compatible with the packages
0485 :filelink:`seaice <pkg/seaice>` and :filelink:`exf <pkg/exf>`. When turned on
0486 together with :filelink:`seaice <pkg/seaice>`, the zero-layer thermodynamics
0487 are replaced by the Winton thermodynamics. In order to use package
0488 :filelink:`seaice <pkg/seaice>` with the thermodynamics of
0489 :filelink:`pkg/thsice`, compile both packages and turn both package on in
0490 ``data.pkg``; see an example in
0491 :filelink:`verification/global_ocean.cs32x15/input.icedyn`. Note, that once
258fe29c91 Jeff*0492 :filelink:`thsice <pkg/thsice>` is turned on, the variables and diagnostics
0493 associated to the default thermodynamics are meaningless, and the diagnostics
0494 of :filelink:`thsice <pkg/thsice>` must be used instead.
adc83e5d7b Mart*0495
0496 .. _para_phys_pkg_seaice_surfaceforcing:
0497
0498 Surface forcing
258fe29c91 Jeff*0499 ---------------
adc83e5d7b Mart*0500
c512e371cc drin*0501 The sea ice model requires the following input fields: 10 m winds, 2 m air
0502 temperature and specific humidity, downward longwave and shortwave radiations,
0503 precipitation, evaporation, and river and glacier runoff. The sea ice model
0504 also requires surface temperature from the ocean model and the top level
0505 horizontal velocity. Output fields are surface wind stress, evaporation minus
0506 precipitation minus runoff, net surface heat flux, and net shortwave flux. The
0507 sea-ice model is global: in ice-free regions bulk formulae (by default computed
0508 in package :filelink:`exf <pkg/exf>`) are used to estimate oceanic forcing from
0509 the atmospheric fields.
adc83e5d7b Mart*0510
a4e168e012 antn*0511 .. _ssub_phys_pkg_seaice_dynamics:
adc83e5d7b Mart*0512
0513 Dynamics
a4e168e012 antn*0514 ========
adc83e5d7b Mart*0515
0516 The momentum equation of the sea-ice model is
0517
0518 .. math::
0bad585a21 Navi*0519 m \frac{D\mathbf{u}}{Dt} = -mf\hat{\mathbf{k}}\times\mathbf{u} +
9c29098ece Jeff*0520 \mathbf{\tau}_\mathrm{air} + \mathbf{\tau}_\mathrm{ocean}
258fe29c91 Jeff*0521 - m \nabla{\phi(0)} + \mathbf{F}
adc83e5d7b Mart*0522 :label: eq_momseaice
0523
14673ec2d0 Mart*0524 where :math:`m=m_{i}+m_{s}` is the ice and snow mass per unit area. The ice
0525 mass per grid cell is :math:`m_i=\rho_{\mathrm{ice}} h\,c` with the mean ice
0526 density :math:`\rho_{\mathrm{ice}}` and the mean thickness :math:`h\,c = `
0527 volume per grid cell area that is the product of the actual thickness :math:`h`
0528 of the ice covered part of the cell and the fractional ice cover :math:`c =
0529 [0,1]`, sloppily also called ice concentration. A similar relationship defines
0530 the snow mass per grid cell :math:`m_s`.
0531 :math:`\mathbf{u}=u\hat{\mathbf{i}}+v\hat{\mathbf{j}}` is the ice velocity
0532 vector; :math:`\hat{\mathbf{i}}`, :math:`\hat{\mathbf{j}}`, and
0533 :math:`\hat{\mathbf{k}}` are unit vectors in the :math:`x`, :math:`y`, and
0534 :math:`z` directions, respectively; :math:`f` is the Coriolis parameter;
0535 :math:`\mathbf{\tau}_\mathrm{air}` and :math:`\mathbf{\tau}_\mathrm{ocean}` are
0536 the wind-ice and ocean-ice stresses, respectively; :math:`g` is the gravity
0537 accelation; :math:`\nabla\phi(0)` is the gradient (or tilt) of the sea surface
0538 height; :math:`\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}` is the sea
0539 surface height potential in response to ocean dynamics (:math:`g\eta`),
0540 atmospheric pressure loading (:math:`p_{a}/\rho_{0}`, where :math:`\rho_{0}` is
0541 a reference density), and a term due to snow and ice loading; and
0542 :math:`\mathbf{F}= \nabla \cdot\sigma` is the divergence of the internal ice
0543 stress tensor :math:`\sigma_{ij}`. Advection of sea-ice momentum is
0544 neglected. The wind and ice-ocean stress terms are given by
adc83e5d7b Mart*0545
0546 .. math::
0547 \begin{aligned}
9c29098ece Jeff*0548 \mathbf{\tau}_\mathrm{air} = & \rho_\mathrm{air} C_\mathrm{air}
c512e371cc drin*0549 |\mathbf{U}_\mathrm{air} -\mathbf{u}| R_\mathrm{air}
0550 (\mathbf{U}_\mathrm{air} - \mathbf{u}) \\
9c29098ece Jeff*0551 \mathbf{\tau}_\mathrm{ocean} = & \rho_\mathrm{ocean}C_\mathrm{ocean}
0552 |\mathbf{U}_\mathrm{ocean}-\mathbf{u}|
c512e371cc drin*0553 R_\mathrm{ocean}(\mathbf{U}_\mathrm{ocean} - \mathbf{u})
adc83e5d7b Mart*0554 \end{aligned}
0555
c512e371cc drin*0556 where :math:`\mathbf{U}_\mathrm{air/ocean}` are the surface winds of the
0557 atmosphere and surface currents of the ocean, respectively;
0558 :math:`C_\mathrm{air/ocean}` are air and ocean drag coefficients;
9c29098ece Jeff*0559 :math:`\rho_\mathrm{air/ocean}` are reference densities; and
0560 :math:`R_\mathrm{air/ocean}` are rotation matrices that act on the wind/current
adc83e5d7b Mart*0561 vectors.
0562
0563 .. _para_phys_pkg_seaice_VPrheology:
0564
0565 Viscous-Plastic (VP) Rheology
258fe29c91 Jeff*0566 -----------------------------
adc83e5d7b Mart*0567
c512e371cc drin*0568 For an isotropic system the stress tensor :math:`\sigma_{ij}` (:math:`i,j=1,2`)
0569 can be related to the ice strain rate and strength by a nonlinear
0570 viscous-plastic (VP) constitutive law:
adc83e5d7b Mart*0571
0572 .. math::
2c231b0ebd Mart*0573 \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
258fe29c91 Jeff*0574 + \left[\zeta(\dot{\epsilon}_{ij},P) -
2c231b0ebd Mart*0575 \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}
258fe29c91 Jeff*0576 - \frac{P}{2}\delta_{ij}
0577 :label: eq_vpequation
adc83e5d7b Mart*0578
0579 The ice strain rate is given by
0580
0581 .. math::
2c231b0ebd Mart*0582 \dot{\epsilon}_{ij} = \frac{1}{2}\left(
adc83e5d7b Mart*0583 \frac{\partial{u_{i}}}{\partial{x_{j}}} +
258fe29c91 Jeff*0584 \frac{\partial{u_{j}}}{\partial{x_{i}}}\right)
adc83e5d7b Mart*0585
c512e371cc drin*0586 The maximum ice pressure :math:`P_{\max}` (variable :varlink:`PRESS0` in the
0587 code), a measure of ice strength, depends on both thickness :math:`h` and
0588 compactness (concentration) :math:`c`:
adc83e5d7b Mart*0589
0590 .. math::
0452697f42 Oliv*0591 :label: eq_icestrength
adc83e5d7b Mart*0592
3b6b5ca15d Mart*0593 P_{\max} = P^{\ast}c\,h\,\exp\{-C^{\ast}\cdot(1-c)\},
adc83e5d7b Mart*0594
2c231b0ebd Mart*0595 with the constants :math:`P^{\ast}` (run-time parameter
9986b4a53e Jeff*0596 :varlink:`SEAICE_strength`) and :math:`C^{\ast}` (run-time parameter
14673ec2d0 Mart*0597 :varlink:`SEAICE_cStar`). Note that Hibler (1979) :cite:`hibler:79` defines
0598 :math:`h` as the "mean thickness" or an "equivalent ice thickness" for mass,
0599 which is :math:`c\,h` with our definitions. By default, :math:`P` (variable
0600 :varlink:`PRESS` in the code) is the replacement pressure
c512e371cc drin*0601
0602 .. math::
0603 :label: eq_pressrepl
0604
0605 P = (1-k_t)\,P_{\max} \left( (1 - f_{r})
0bad585a21 Navi*0606 + f_{r} \frac{\Delta}{\Delta_{\rm reg}} \right)
c512e371cc drin*0607
14673ec2d0 Mart*0608 where :math:`f_{r}` is a run-time parameter :varlink:`SEAICEpressReplFac`
0bad585a21 Navi*0609 (default = 1.0), and :math:`\Delta_{\rm reg}` is a regularized form of
c512e371cc drin*0610 :math:`\Delta = \left[ \left(\dot{\epsilon}_{11}+\dot{\epsilon}_{22}\right)^2 +
0611 e^{-2}\left( \left(\dot{\epsilon}_{11}-\dot{\epsilon}_{22} \right)^2 +
14673ec2d0 Mart*0612 4\,\dot{\epsilon}_{12}^2 \right) \right]^{\frac{1}{2}}`. By default
0613 :math:`\Delta_{\mathrm{reg}}=\max(\Delta,\Delta_{\min})`. If CPP-flag
0614 :varlink:`SEAICE_DELTA_SMOOTHREG` is defined,
0615 :math:`\Delta_{\mathrm{reg}}=\sqrt{\Delta^2+\Delta^2_{\min}}`. Run-time
0616 parameter :varlink:`SEAICE_deltaMin` :math:`= \Delta_{\min} = 10^{-10}` by
0617 default.
c512e371cc drin*0618
0619 The tensile strength factor :math:`k_t` (run-time parameter
0620 :varlink:`SEAICE_tensilFac`) determines the ice tensile strength :math:`T =
0621 k_t\cdot P_{\max}`, as defined by König Beatty and Holland (2010)
0622 :cite:`konig:10`. :varlink:`SEAICE_tensilFac` is zero by default.
0623
0624 Different VP rheologies can be used to model sea ice dynamics. The different
0625 rheologies are characterized by different definitions of the bulk and shear
0626 viscosities :math:`\zeta` and :math:`\eta` in :eq:`eq_vpequation`. The
0627 following :numref:`tab_phys_pkg_seaice_rheologies` is a summary of the
0628 available choices with recommended (sensible) parameter values. All the
0629 rheologies presented here depend on the ice strength :math:`P`
0630 :eq:`eq_pressrepl`.
0631
0632 .. tabularcolumns:: |\Y{.275}|\Y{.450}|\Y{.275}|
0633
0634 .. table:: Overview over availabe sea ice viscous-plastic rheologies
0635 :class: longtable
0636 :name: tab_phys_pkg_seaice_rheologies
0637
0638 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0639 | Name | CPP flags | Run-time flags (recommended value) |
0640 +=======================================+=======================================+====================================================+
0641 | :ref:`rheologies_ellnfr` | None (default) | - :varlink:`SEAICE_eccen` (= 2.0) |
0642 | | | - :varlink:`SEAICE_tensilFac` (= 0.0) |
0643 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0644 | :ref:`rheologies_ellnnfr` | None | - :varlink:`SEAICE_eccen` (= 2.0) |
0645 | | | - :varlink:`SEAICE_eccfr` (< 2.0) |
0646 | | | - :varlink:`SEAICE_tensilFac` (= 0.0) |
0647 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0648 | :ref:`rheologies_TEM` | :varlink:`SEAICE_ALLOW_TEM` | - :varlink:`SEAICEuseTEM` (=.TRUE.) |
0649 | | | - :varlink:`SEAICE_eccen` (= 1.4) |
0650 | | | - :varlink:`SEAICE_eccfr` (< 1.4) |
0651 | | | - :varlink:`SEAICE_tensilFac` (= 0.05) |
0652 | | | - :varlink:`SEAICEmcMU` (= 0.6 to 0.8) |
0653 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0654 | :ref:`rheologies_MCE` | :varlink:`SEAICE_ALLOW_MCE` | - :varlink:`SEAICEuseMCE` (=.TRUE.) |
0655 | | | - :varlink:`SEAICE_eccen` (= 1.4) |
0656 | | | - :varlink:`SEAICE_eccfr` (< 1.4) |
0657 | | | - :varlink:`SEAICE_tensilFac` (= 0.05) |
0658 | | | - :varlink:`SEAICEmcMU` (= 0.6 to 0.8) |
0659 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0660 | :ref:`rheologies_MCS` | :varlink:`SEAICE_ALLOW_MCS` | - :varlink:`SEAICEuseMCS` (=.TRUE.) |
0661 | | | - :varlink:`SEAICE_tensilFac` (= 0.05) |
0662 | | | - :varlink:`SEAICEmcMU` (= 0.6 to 0.8) |
0663 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0664 | :ref:`rheologies_TD` | :varlink:`SEAICE_ALLOW_TD` | - :varlink:`SEAICEuseTD` (=.TRUE.) |
0665 | | | - :varlink:`SEAICE_tensilFac` (= 0.025) |
0666 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0667 | :ref:`rheologies_PL` | :varlink:`SEAICE_ALLOW_TD` | - :varlink:`SEAICEusePL` (=.TRUE.) |
0668 | | | - :varlink:`SEAICE_tensilFac` (= 0.025) |
0669 +---------------------------------------+---------------------------------------+----------------------------------------------------+
0670
0671
0672 **Note:** With the exception of the default rheology and the TEM (with
0673 :varlink:`SEAICEmcMU` : :math:`\mu=1.0`), these rheologies are not implemented
0674 in EVP (:numref:`para_phys_pkg_seaice_EVPdynamics`).
0675
0676 .. _rheologies_ellnfr:
0677
0678 Elliptical yield curve with normal flow rule
0679 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0680
0681 The default rheology in the sea ice module of the MITgcm implements the widely
0682 used elliptical yield curve with a normal flow rule :cite:`hibler:79`. For
0683 this yield curve, the nonlinear bulk and shear viscosities :math:`\zeta` and
0684 :math:`\eta` are functions of ice strain rate invariants and ice strength such
0685 that the principal components of the stress lie on an elliptical yield curve
0686 with the ratio of major to minor axis :math:`e = 2.0` (run-time parameter
0687 :varlink:`SEAICE_eccen`); they are given by:
adc83e5d7b Mart*0688
0689 .. math::
0690 \begin{aligned}
14673ec2d0 Mart*0691 \zeta =& \min\left(\frac{(1+k_t)P_{\max}}{2\Delta_\mathrm{reg}},
adc83e5d7b Mart*0692 \zeta_{\max}\right) \\
c512e371cc drin*0693 \eta =& \frac{\zeta}{e^2}
0694 \end{aligned}
0695 :label: eq_zetareg
258fe29c91 Jeff*0696
0697
0698 with the abbreviation
0699
0700 .. math::
0701 \Delta = \left[
c61841e2fd Jeff*0702 \left(\dot{\epsilon}_{11}+\dot{\epsilon}_{22}\right)^2
0703 + e^{-2}\left( \left(\dot{\epsilon}_{11}-\dot{\epsilon}_{22} \right)^2
14673ec2d0 Mart*0704 + 4\,\dot{\epsilon}_{12}^2 \right)
c61841e2fd Jeff*0705 \right]^{\frac{1}{2}}
adc83e5d7b Mart*0706
0707 The bulk viscosities are bounded above by imposing both a minimum
14673ec2d0 Mart*0708 :math:`\Delta_{\min}` and replacing :math:`\Delta` by the regularized version
0709 :math:`\Delta_\mathrm{reg}` (for historical reasons, run-time parameter
c512e371cc drin*0710 :varlink:`SEAICE_deltaMin` is set to a default value of
0711 :math:`10^{-10}\,\text{s}^{-1}`, the value of :varlink:`SEAICE_EPS`) and a
0712 maximum :math:`\zeta_{\max} = P_{\max}/(2\Delta^\ast)`, where
14673ec2d0 Mart*0713 :math:`\Delta^\ast=(2\times10^4/5\times10^{12})\,\text{s}^{-1} =
0714 2\times10^{-9}\,\text{s}^{-1}` (:varlink:`SEAICE_zetaMaxFac`
0715 :math:`=\frac{1}{2\Delta^\ast}`). Obviously, this corresponds to regularizing
c512e371cc drin*0716 :math:`\Delta` with the typical value of :varlink:`SEAICE_deltaMin` :math:`=
0717 2\times10^{-9}`. Clearly, some of this regularization is redundant. (There is
0718 also the option of bounding :math:`\zeta` from below by setting run-time
0719 parameter :varlink:`SEAICE_zetaMin` :math:`>0`, but this is generally not
0720 recommended). For stress tensor computation the replacement pressure :math:`P =
0721 2\,\Delta\zeta` is used so that the stress state always lies on the elliptic
0722 yield curve by definition.
adc83e5d7b Mart*0723
258fe29c91 Jeff*0724 Defining the CPP-flag :varlink:`SEAICE_ZETA_SMOOTHREG` in
c512e371cc drin*0725 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` before compiling
0726 replaces the method for bounding :math:`\zeta` by a smooth (differentiable)
0727 expression:
adc83e5d7b Mart*0728
0729 .. math::
258fe29c91 Jeff*0730 \begin{split}
c512e371cc drin*0731 \zeta &= \zeta_{\max}\tanh\left(\frac{(1+k_t)P_{\max}}{2\,
14673ec2d0 Mart*0732 \Delta_\mathrm{reg} \,\zeta_{\max}}\right)\\
c512e371cc drin*0733 &= \frac{(1+k_t)P_{\max}}{2\Delta^\ast}
14673ec2d0 Mart*0734 \tanh\left(\frac{\Delta^\ast}{\Delta_\mathrm{reg}}\right)
258fe29c91 Jeff*0735 \end{split}
0736 :label: eq_zetaregsmooth
adc83e5d7b Mart*0737
c61841e2fd Jeff*0738 where :math:`\Delta_{\min}=10^{-20}\,\text{s}^{-1}` should be chosen to avoid
adc83e5d7b Mart*0739 divisions by zero.
0740
c512e371cc drin*0741 In this default formulation the yield curve does not allow isotropic tensile
0742 stress, that is, sea ice can be "pulled apart" without any effort. Setting the
0743 parameter :math:`k_t` (:varlink:`SEAICE_tensilFac`) to a small value larger
0744 than zero, extends the yield curve into a region where the divergence of the
0745 stress :math:`\sigma_{11}+\sigma_{22} > 0` to allow some tensile stress.
0746
0747 Besides this commonly used default rheology, a number of a alternative
0748 rheologies are implemented. Some of these are experiemental and should be used
0749 with caution.
0750
0751 .. _rheologies_ellnnfr:
0752
0753 Elliptical yield curve with non-normal flow rule
0754 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0755
0756 Defining the run-time parameter :varlink:`SEAICE_eccfr` with a value different
0757 from :varlink:`SEAICE_eccen` allows one to use an elliptical yield curve with a
0758 non-normal flow rule as described in Ringeisen et al. (2020)
0759 :cite:`ringeisen:20`. In this case the viscosities are functions of
0760 :math:`e_F` (:varlink:`SEAICE_eccen`) and :math:`e_G`
0761 (:varlink:`SEAICE_eccfr`):
0762
0763 .. math::
0764 \begin{aligned}
0765 \zeta &= \frac{P_{\max}(1+k_t)}{2\Delta} \\
0766 \eta &= \frac{\zeta}{e_G^2} = \frac{P_{\max}(1+k_t)}{2e_G^2\Delta}
0767 \end{aligned}
0768
0769 with the abbreviation
0770
0771 .. math::
0772 \Delta = \sqrt{(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2
0773 +\frac{e_F^2}{e_G^4}((\dot{\epsilon}_{11}
14673ec2d0 Mart*0774 -\dot{\epsilon}_{22})^2+4\,\dot{\epsilon}_{12}^2)}.
c512e371cc drin*0775
0776 Note that if :math:`e_G=e_F=e`, these formulae reduce to the normal flow rule.
0777
0778 .. _rheologies_TEM:
0779
0780 Truncated ellipse method (TEM) for elliptical yield curve
0781 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0782
0783 In the so-called truncated ellipse method, the shear viscosity :math:`\eta` is
0784 capped to suppress any tensile stress:
0785
0786 .. math::
0787 \eta = \min\left(\frac{\zeta}{e^2},
0788 \frac{\frac{(1+k_t)\,P_{\max}}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
0789 {\sqrt{\max(\Delta_{\min}^{2},(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2
0790 +4\dot{\epsilon}_{12}^2})}\right).
0791 :label: eq_etatem
0792
0793 To enable this method, set ``#define`` :varlink:`SEAICE_ALLOW_TEM` in
0794 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` and turn it on with
dc26f158aa Mart*0795 :varlink:`SEAICEuseTEM` ``=.TRUE.,`` in ``data.seaice``. This parameter
c512e371cc drin*0796 combination implies the default of :varlink:`SEAICEmcMU` :math:`= 1.0`.
0797
0798 Instead of an ellipse that is truncated by constant slope coulombic limbs, this
0799 yield curve can also be seen as a Mohr-Coulomb yield curve with elliptical flow
0800 rule that is truncated for high :math:`P` by an ellipse. As a consequence, the
0801 Mohr-Coulomb slope :varlink:`SEAICEmcMU` can be set in ``data.seaice`` to
0802 values :math:`\ne 1.0`. This defines a coulombic yield curve similar to the
0803 ones shown in Hibler and Schulson (2000) :cite:`hibler:00` and Ringeisen et
0804 al. (2019) :cite:`ringeisen:19`.
0805
0806 For this rheology, it is recommended to use a non-zero tensile strength, so set
0807 :varlink:`SEAICE_tensilFac` :math:`=k_{t}>0` in ``data.seaice``, e.g., :math:`=
0808 0.05` or 5%.
0809
0810 .. _rheologies_MCE:
0811
0812 Mohr-Coulomb yield curve with elliptical plastic potential
0813 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0814
0815 To use a Mohr-Coulomb rheology, set ``#define`` :varlink:`SEAICE_ALLOW_MCE` in
0816 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` and
dc26f158aa Mart*0817 :varlink:`SEAICEuseMCE` ``= .TRUE.,`` in ``data.seaice``. This Mohr-Coulomb
c512e371cc drin*0818 yield curve uses an elliptical plastic potential to define the flow rule. The
0819 slope of the Mohr-Coulomb yield curve is defined by :varlink:`SEAICEmcMU` in
0820 ``data.seaice``, and the plastic potential ellipse aspect ratio is set by
0821 :varlink:`SEAICE_eccfr` in ``data.seaice``. For details of this rheology, see
0822 https://doi.org/10.26092/elib/380, Chapter 2.
0823
0824 For this rheology, it is recommended to use a non-zero tensile strength, so set
0825 :varlink:`SEAICE_tensilFac` :math:`>0` in ``data.seaice``, e.g., :math:`= 0.05`
0826 or 5%.
0827
0828 .. _rheologies_MCS:
0829
0830 Mohr-Coulomb yield curve with shear flow rule
0831 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0832
0833 To use the specifc Mohr-Coulomb rheology as defined first by Ip et al. (1991)
0834 :cite:`ip:91`, set ``#define`` :varlink:`SEAICE_ALLOW_MCS` in
0835 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` and
dc26f158aa Mart*0836 :varlink:`SEAICEuseMCS` ``= .TRUE.,`` in ``data.seaice``. The slope of the
c512e371cc drin*0837 Mohr-Coulomb yield curve is defined by :varlink:`SEAICEmcMU` in
0838 ``data.seaice``. For details of this rheology, including the tensile strength,
0839 see https://doi.org/10.26092/elib/380, Chapter 2.
0840
0841 For this rheology, it is recommended to use a non-zero tensile strength, so set
0842 :varlink:`SEAICE_tensilFac` :math:`>0` in ``data.seaice``, e.g., :math:`= 0.05`
0843 or 5%.
0844
0845 **WARNING: This rheology is known to be unstable. Use with caution!**
0846
0847 .. _rheologies_TD:
0848
0849 Teardrop yield curve with normal flow rule
0850 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0851
0852 The teardrop rheology was first described in Zhang and Rothrock (2005)
0853 :cite:`zha:05`. Here we implement a slightly modified version (See
0854 https://doi.org/10.26092/elib/380, Chapter 2).
0855
0856 To use this rheology, set ``#define`` :varlink:`SEAICE_ALLOW_TEARDROP` in
0857 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` and
dc26f158aa Mart*0858 :varlink:`SEAICEuseTD` ``= .TRUE.,`` in ``data.seaice``. The size of the yield
c512e371cc drin*0859 curve can be modified by changing the tensile strength, using
0860 :varlink:`SEAICE_tensFac` in ``data.seaice``.
0861
0862 For this rheology, it is recommended to use a non-zero tensile strength, so set
0863 :varlink:`SEAICE_tensilFac` :math:`>0` in ``data.seaice``, e.g., :math:`=
0864 0.025` or 2.5%.
0865
0866 .. _rheologies_PL:
0867
0868 Parabolic lens yield curve with normal flow rule
0869 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0870
0871 The parabolic lens rheology was first described in Zhang and Rothrock (2005)
0872 :cite:`zha:05`. Here we implement a slightly modified version (See
0873 https://doi.org/10.26092/elib/380, Chapter 2).
0874
0875 To use this rheology, set ``#define`` :varlink:`SEAICE_ALLOW_TEARDROP` in
0876 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` and
dc26f158aa Mart*0877 :varlink:`SEAICEusePL` ``= .TRUE.,`` in ``data.seaice``. The size of the yield
c512e371cc drin*0878 curve can be modified by changing the tensile strength, using
0879 :varlink:`SEAICE_tensFac` in ``data.seaice``.
0880
0881 For this rheology, it is recommended to use a non-zero tensile strength, so set
0882 :varlink:`SEAICE_tensilFac` :math:`>0` in ``data.seaice``, e.g., :math:`=
0883 0.025` or 2.5%.
0884
adc83e5d7b Mart*0885 .. _para_phys_pkg_seaice_LSRJFNK:
0886
0887 LSR and JFNK solver
258fe29c91 Jeff*0888 -------------------
adc83e5d7b Mart*0889
c512e371cc drin*0890 In matrix notation, the discretized momentum equations can be written as
adc83e5d7b Mart*0891
0892 .. math::
0893 :label: eq_matrixmom
2c231b0ebd Mart*0894
c61841e2fd Jeff*0895 \mathbf{A}(\mathbf{x})\,\mathbf{x} = \mathbf{b}(\mathbf{x}).
adc83e5d7b Mart*0896
c512e371cc drin*0897 The solution vector :math:`\mathbf{x}` consists of the two velocity components
0898 :math:`u` and :math:`v` that contain the velocity variables at all grid points
0899 and at one time level. The standard (and default) method for solving
0900 Eq. :eq:`eq_matrixmom` in the sea ice component of MITgcm is an iterative
0901 Picard solver: in the :math:`k`-th iteration a linearized form
adc83e5d7b Mart*0902 :math:`\mathbf{A}(\mathbf{x}^{k-1})\,\mathbf{x}^{k} =
c512e371cc drin*0903 \mathbf{b}(\mathbf{x}^{k-1})` is solved (in the case of MITgcm it is a Line
0904 Successive (over) Relaxation (LSR) algorithm). Picard solvers converge slowly,
0905 but in practice the iteration is generally terminated after only a few
0906 nonlinear steps and the calculation continues with the next time level. This
0907 method is the default method in MITgcm. The number of nonlinear iteration steps
0908 or pseudo-time steps can be controlled by the run-time parameter
dc26f158aa Mart*0909 :varlink:`SEAICEnonLinIterMax`. This parameter's default is 2, but using a
0910 number of at least 10 is recommended for better solutions that are converged at
0911 least in an energy norm sense (Zhang and Hibler 1997) :cite:`zhang:97`.
c512e371cc drin*0912
dc26f158aa Mart*0913 In order to overcome the poor convergence of the Picard solver, Lemieux et
c512e371cc drin*0914 al. (2010) :cite:`lemieux:10` introduced a Jacobian-free Newton-Krylov solver
0915 for the sea ice momentum equations. This solver is also implemented in MITgcm
0916 (see Losch et al. 2014 :cite:`losch:14`). The Newton method transforms
0917 minimizing the residual :math:`\mathbf{F}(\mathbf{x}) =
0918 \mathbf{A}(\mathbf{x})\,\mathbf{x} - \mathbf{b}(\mathbf{x})` to finding the
0919 roots of a multivariate Taylor expansion of the residual :math:`\mathbf{F}`
0920 around the previous (:math:`k-1`) estimate :math:`\mathbf{x}^{k-1}`:
adc83e5d7b Mart*0921
0922 .. math::
258fe29c91 Jeff*0923 \mathbf{F}(\mathbf{x}^{k-1}+\delta\mathbf{x}^{k}) =
0924 \mathbf{F}(\mathbf{x}^{k-1}) + \mathbf{F}'(\mathbf{x}^{k-1})
0925 \,\delta\mathbf{x}^{k}
adc83e5d7b Mart*0926 :label: eq_jfnktaylor
0927
c512e371cc drin*0928 with the Jacobian :math:`\mathbf{J}\equiv\mathbf{F}'`. The root
0929 :math:`\mathbf{F}(\mathbf{x}^{k-1}+\delta\mathbf{x}^{k})=0` is found by solving
adc83e5d7b Mart*0930
0931 .. math::
258fe29c91 Jeff*0932 \mathbf{J}(\mathbf{x}^{k-1})\,\delta\mathbf{x}^{k} =
0933 -\mathbf{F}(\mathbf{x}^{k-1})
adc83e5d7b Mart*0934 :label: eq_jfnklin
0935
c512e371cc drin*0936 for :math:`\delta\mathbf{x}^{k}`. The next (:math:`k`-th) estimate is given by
c704c5a1ef Mart*0937 :math:`\mathbf{x}^{k}=\mathbf{x}^{k-1}+(1-\gamma_{\mathrm{LS}})^{l}
0938 \,\delta\mathbf{x}^{k}`.
0939
0940 By default :math:`l=0`, but in order to avoid overshoots, the step size factor
0941 :math:`(1-\gamma_{\mathrm{LS}})^{l}` with :math:`\gamma_{\mathrm{LS}}<1` can be
0942 iteratively reduced in a line search with :math:`l=0,1,2,\ldots` until
c512e371cc drin*0943 :math:`\|\mathbf{F}(\mathbf{x}^k)\| < \|\mathbf{F}(\mathbf{x}^{k-1})\|`, where
c704c5a1ef Mart*0944 :math:`\|\cdot\|=\int\cdot\,dx^2` is the :math:`L_2`-norm. The line search
0945 starts after :varlink:`SEAICE_JFNK_lsIter` nonlinear Newton iterations (off by
0946 default) to allow for full Newton steps at the beginning of the iteration. If
0947 the line search is turned on by setting :varlink:`SEAICE_JFNK_lsIter` to a
0948 non-negative value in ``data.seaice``, by default, the line search with
0949 :math:`\gamma_\mathrm{LS}=\frac{1}{2}` (runtime parameter
0950 :varlink:`SEAICE_JFNK_lsGamma`) is stopped after :math:`L_{\max}=4` (runtime
0951 parameter :varlink:`SEAICE_JFNK_lsLmax`) steps.
c512e371cc drin*0952
0953 Forming the Jacobian :math:`\mathbf{J}` explicitly is often avoided as “too
0954 error prone and time consumingâ€. Instead, Krylov methods only require the
0955 action of :math:`\mathbf{J}` on an arbitrary vector :math:`\mathbf{w}` and
0956 hence allow a matrix free algorithm for solving :eq:`eq_jfnklin`. The action of
0957 :math:`\mathbf{J}` can be approximated by a first-order Taylor series
0958 expansion:
adc83e5d7b Mart*0959
0960 .. math::
258fe29c91 Jeff*0961 \mathbf{J}(\mathbf{x}^{k-1})\,\mathbf{w} \approx
0962 \frac{\mathbf{F}(\mathbf{x}^{k-1}+\epsilon\mathbf{w})
0963 - \mathbf{F}(\mathbf{x}^{k-1})} \epsilon
adc83e5d7b Mart*0964 :label: eq_jfnkjacvecfd
0965
0966 or computed exactly with the help of automatic differentiation (AD)
258fe29c91 Jeff*0967 tools. :varlink:`SEAICE_JFNKepsilon` sets the step size :math:`\epsilon`.
adc83e5d7b Mart*0968
258fe29c91 Jeff*0969 We use the Flexible Generalized Minimum RESidual (FMGRES) method with
c512e371cc drin*0970 right-hand side preconditioning to solve :eq:`eq_jfnklin` iteratively starting
0971 from a first guess of :math:`\delta\mathbf{x}^{k}_{0} = 0`. For the
0972 preconditioning matrix :math:`\mathbf{P}` we choose a simplified form of the
0973 system matrix :math:`\mathbf{A}(\mathbf{x}^{k-1})` where
0974 :math:`\mathbf{x}^{k-1}` is the estimate of the previous Newton step
0975 :math:`k-1`. The transformed equation :eq:`eq_jfnklin` becomes
adc83e5d7b Mart*0976
0977 .. math::
0978 \mathbf{J}(\mathbf{x}^{k-1})\,\mathbf{P}^{-1}\delta\mathbf{z} =
0979 -\mathbf{F}(\mathbf{x}^{k-1}), \quad\text{with} \quad
258fe29c91 Jeff*0980 \delta{\mathbf{z}} = \mathbf{P}\delta\mathbf{x}^{k}
0981 :label: eq_jfnklinpc
adc83e5d7b Mart*0982
c512e371cc drin*0983 The Krylov method iteratively improves the approximate solution to
0984 :eq:`eq_jfnklinpc` in subspace (:math:`\mathbf{r}_0`,
0985 :math:`\mathbf{J}\mathbf{P}^{-1}\mathbf{r}_0`,
2c231b0ebd Mart*0986 :math:`(\mathbf{J}\mathbf{P}^{-1})^2\mathbf{r}_0`, :math:`\dots`,
c512e371cc drin*0987 :math:`(\mathbf{J}\mathbf{P}^{-1})^m\mathbf{r}_0`) with increasing :math:`m`;
0988 :math:`\mathbf{r}_0 = -\mathbf{F}(\mathbf{x}^{k-1})
0989 -\mathbf{J}(\mathbf{x}^{k-1})\,\delta\mathbf{x}^{k}_{0}` is the initial
0990 residual of :eq:`eq_jfnklin`;
0991 :math:`\mathbf{r}_0=-\mathbf{F}(\mathbf{x}^{k-1})` with the first guess
dc26f158aa Mart*0992 :math:`\delta\mathbf{x}^{k}_{0}=0`. We allow a Krylov subspace of dimension \
c512e371cc drin*0993 :math:`m=50` and we do allow restarts for more than 50 Krylov iterations. The
0994 preconditioning operation involves applying :math:`\mathbf{P}^{-1}` to the
0995 basis vectors :math:`\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2, \ldots,
0996 \mathbf{v}_m` of the Krylov subspace. This operation is approximated by solving
0997 the linear system :math:`\mathbf{P}\,\mathbf{w}=\mathbf{v}_i`. Because
0998 :math:`\mathbf{P} \approx \mathbf{A}(\mathbf{x}^{k-1})`, we can use the
dc26f158aa Mart*0999 LSR algorithm already implemented in the Picard solver. Each preconditioning
1000 operation uses a fixed number of 10 LSR iterations avoiding any termination
c512e371cc drin*1001 criterion. More details and results can be found in Losch et al. (2014)
1002 :cite:`losch:14`).
1003
dc26f158aa Mart*1004 To use the JFNK solver set :varlink:`SEAICEuseJFNK` ``= .TRUE.,`` in the
c512e371cc drin*1005 namelist file ``data.seaice``; ``#define`` :varlink:`SEAICE_ALLOW_JFNK` in
1006 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` and we recommend
1007 using a smooth regularization of :math:`\zeta` by ``#define``
1008 :varlink:`SEAICE_ZETA_SMOOTHREG` (see above) for better convergence. The
1009 nonlinear Newton iteration is terminated when the :math:`L_2`-norm of the
1010 residual is reduced by :math:`\gamma_{\mathrm{nl}}` (run-time parameter
1011 :varlink:`SEAICEnonLinTol` ``= 1.E-4,`` will already lead to expensive
1012 simulations) with respect to the initial norm:
1013 :math:`\|\mathbf{F}(\mathbf{x}^k)\| <
1014 \gamma_{\mathrm{nl}}\|\mathbf{F}(\mathbf{x}^0)\|`. Within a nonlinear
1015 iteration, the linear FGMRES solver is terminated when the residual is smaller
1016 than :math:`\gamma_k\|\mathbf{F}(\mathbf{x}^{k-1})\|` where :math:`\gamma_k` is
1017 determined by
adc83e5d7b Mart*1018
1019 .. math::
2c231b0ebd Mart*1020 \gamma_k =
1021 \begin{cases}
1022 \gamma_0 &\text{for $\|\mathbf{F}(\mathbf{x}^{k-1})\| \geq r$}, \\
258fe29c91 Jeff*1023 \max\left(\gamma_{\min},
1024 \frac{\|\mathbf{F}(\mathbf{x}^{k-1})\|}
2c231b0ebd Mart*1025 {\|\mathbf{F}(\mathbf{x}^{k-2})\|}\right)
258fe29c91 Jeff*1026 &\text{for $\|\mathbf{F}(\mathbf{x}^{k-1})\| < r$,}
1027 \end{cases}
1028 :label: eq_jfnkgammalin
adc83e5d7b Mart*1029
c512e371cc drin*1030 so that the linear tolerance parameter :math:`\gamma_k` decreases with the
1031 nonlinear Newton step as the nonlinear solution is approached. This inexact
1032 Newton method is generally more robust and computationally more efficient than
1033 exact methods. Typical parameter choices are :math:`\gamma_0 =`
1034 :varlink:`JFNKgamma_lin_max` :math:`= 0.99`, :math:`\gamma_{\min} =`
1035 :varlink:`JFNKgamma_lin_min` :math:`= 0.1`, and :math:`r =`
258fe29c91 Jeff*1036 :varlink:`JFNKres_tFac` :math:`\times\|\mathbf{F}(\mathbf{x}^{0})\|` with
1037 :varlink:`JFNKres_tFac` :math:`= 0.5`. We recommend a maximum number of
c512e371cc drin*1038 nonlinear iterations :varlink:`SEAICEnewtonIterMax` :math:`= 100` and a maximum
1039 number of Krylov iterations :varlink:`SEAICEkrylovIterMax` :math:`= 50`,
1040 because the Krylov subspace has a fixed dimension of 50 (but restarts are
1041 allowed for :varlink:`SEAICEkrylovIterMax` :math:`> 50`).
adc83e5d7b Mart*1042
dc26f158aa Mart*1043 Setting :varlink:`SEAICEuseStrImpCpl` to ``.TRUE.`` turns on “strength implicit
1044 coupling†(see Hutchings et al. 2004 :cite:`hutchings:04`) in the LSR solver
1045 and in the LSR preconditioner for the JFNK solver. In this mode, the different
c512e371cc drin*1046 contributions of the stress divergence terms are reordered so as to increase
1047 the diagonal dominance of the system matrix. Unfortunately, the convergence
dc26f158aa Mart*1048 rate of the LSR solver is increased only slightly, while the JFNK convergence
c512e371cc drin*1049 appears to be unaffected.
adc83e5d7b Mart*1050
1051 .. _para_phys_pkg_seaice_EVPdynamics:
1052
1053 Elastic-Viscous-Plastic (EVP) Dynamics
258fe29c91 Jeff*1054 --------------------------------------
adc83e5d7b Mart*1055
c512e371cc drin*1056 Hunke and Dukowicz (1997) :cite:`hunke:97` introduced an elastic contribution
1057 to the strain rate in order to regularize :eq:`eq_vpequation` in such a way
1058 that the resulting elastic-viscous-plastic (EVP) and VP models are identical at
1059 steady state,
adc83e5d7b Mart*1060
1061 .. math::
258fe29c91 Jeff*1062 \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
2c231b0ebd Mart*1063 \frac{1}{2\eta}\sigma_{ij}
1064 + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}
258fe29c91 Jeff*1065 + \frac{P}{4\zeta}\delta_{ij}
1066 = \dot{\epsilon}_{ij}.
adc83e5d7b Mart*1067 :label: eq_evpequation
1068
dc26f158aa Mart*1069 The EVP model uses an explicit time stepping scheme with a short timestep.
c512e371cc drin*1070 According to the recommendation in Hunke and Dukowicz (1997) :cite:`hunke:97`,
258fe29c91 Jeff*1071 the EVP-model should be stepped forward in time 120 times
1072 (:varlink:`SEAICE_deltaTevp` = :varlink:`SEAICE_deltaTdyn` /120) within the
c512e371cc drin*1073 physical ocean model time step (although this parameter is under debate), to
1074 allow for elastic waves to disappear. Because the scheme does not require a
258fe29c91 Jeff*1075 matrix inversion it is fast in spite of the small internal timestep and simple
1076 to implement on parallel computers. For completeness, we repeat the equations
c512e371cc drin*1077 for the components of the stress tensor :math:`\sigma_{1} =
1078 \sigma_{11}+\sigma_{22}`, :math:`\sigma_{2}= \sigma_{11}-\sigma_{22}`, and
1079 :math:`\sigma_{12}`. Introducing the divergence :math:`D_D =
1080 \dot{\epsilon}_{11}+\dot{\epsilon}_{22}`, and the horizontal tension and
1081 shearing strain rates, :math:`D_T = \dot{\epsilon}_{11}-\dot{\epsilon}_{22}`
1082 and :math:`D_S = 2\dot{\epsilon}_{12}`, respectively, and using the above
1083 abbreviations, the equations :eq:`eq_evpequation` can be written as:
adc83e5d7b Mart*1084
1085 .. math::
258fe29c91 Jeff*1086 \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
1087 \frac{P}{2T} = \frac{P}{2T\Delta} D_D
1088 :label: eq_evpstresstensor1
0452697f42 Oliv*1089
1090 .. math::
258fe29c91 Jeff*1091 \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
1092 = \frac{P}{2T\Delta} D_T
1093 :label: eq_evpstresstensor2
0452697f42 Oliv*1094
1095 .. math::
258fe29c91 Jeff*1096 \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
1097 = \frac{P}{4T\Delta} D_S
1098 :label: eq_evpstresstensor12
adc83e5d7b Mart*1099
1100 Here, the elastic parameter :math:`E` is redefined in terms of a damping
1101 timescale :math:`T` for elastic waves
1102
258fe29c91 Jeff*1103 .. math:: E=\frac{\zeta}{T}
adc83e5d7b Mart*1104
c512e371cc drin*1105 :math:`T=E_{0}\Delta{t}` with the tunable parameter :math:`E_0<1` and the
1106 external (long) timestep :math:`\Delta{t}`. :math:`E_{0} = \frac{1}{3}` is the
dc26f158aa Mart*1107 default value in the code and close to what Hunke and Dukowicz (1997)
1108 :cite:`hunke:97` recommend.
c512e371cc drin*1109
dc26f158aa Mart*1110 We do not recommend to use the EVP solver in its original form. Instead, use
1111 mEVP or aEVP instead (see :numref:`para_phys_pkg_seaice_EVPstar`). If you
1112 really need to use the original EVP solver, make sure that both ``#define``
1113 :varlink:`SEAICE_CGRID` and ``#define`` :varlink:`SEAICE_ALLOW_EVP` are set in
c512e371cc drin*1114 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` (both are defined by
dc26f158aa Mart*1115 default). By default, the runtime parameters :varlink:`SEAICEuseEVPstar` and
1116 :varlink:`SEAICEuseEVPrev` are set to ``.TRUE.``, which already improves the
1117 behavoir of EVP, but for the original EVP they should be set to ``.FALSE.``. The
1118 solver is turned on by setting the sub-cycling time step
c512e371cc drin*1119 :varlink:`SEAICE_deltaTevp` to a value larger than zero. The choice of this
dc26f158aa Mart*1120 time step is under debate. Hunke and Dukowicz (1997) :cite:`hunke:97` recommend
1121 order 120 time steps for the EVP solver within one model time step
258fe29c91 Jeff*1122 :math:`\Delta{t}` (:varlink:`deltaTmom`). One can also choose order 120 time
c512e371cc drin*1123 steps within the forcing time scale, but then we recommend adjusting the
1124 damping time scale :math:`T` accordingly, by setting either
1125 :varlink:`SEAICE_elasticParm` (:math:`E_{0}`), so that :math:`E_{0}\Delta{t}=`
1126 forcing time scale, or directly :varlink:`SEAICE_evpTauRelax` (:math:`T`) to
1127 the forcing time scale. (**NOTE**: with the improved EVP variants of the next
1128 section, the above recommendations are obsolete. Use mEVP or aEVP instead.)
adc83e5d7b Mart*1129
1130 .. _para_phys_pkg_seaice_EVPstar:
1131
1c8cebb321 Jeff*1132 More stable variants of Elastic-Viscous-Plastic Dynamics: EVP\*, mEVP, and aEVP
1133 -------------------------------------------------------------------------------
adc83e5d7b Mart*1134
c512e371cc drin*1135 The genuine EVP scheme appears to give noisy solutions (see Hunke 2001, Lemieux
1136 et al. 2012, Bouillon et a1. 2013
1137 :cite:`hunke:01,lemieux:12,bouillon:13`). This has led to a modified EVP or
1138 EVP\* (Lemieux et al. 2012, Bouillon et a1. 2013, Kimmritz et al. 2015
1139 :cite:`lemieux:12,bouillon:13,kimmritz:15`); here, we refer to these variants
1140 by modified EVP (mEVP) and adaptive EVP (aEVP). The main idea is to modify the
1141 “natural†time-discretization of the momentum equations:
adc83e5d7b Mart*1142
1143 .. math::
258fe29c91 Jeff*1144 m\frac{D\mathbf{u}}{Dt} \approx
1145 m\frac{\mathbf{u}^{p+1}-\mathbf{u}^{n}}{\Delta{t}} +
1146 \beta^{\ast}\frac{\mathbf{u}^{p+1}-\mathbf{u}^{p}}{\Delta{t}_{\mathrm{EVP}}}
adc83e5d7b Mart*1147 :label: eq_evpstar
1148
c512e371cc drin*1149 where :math:`n` is the previous time step index, and :math:`p` is the previous
1150 sub-cycling index. The extra “intertial†term
1151 :math:`m\,(\mathbf{u}^{p+1}-\mathbf{u}^{n})/\Delta{t})` allows the definition
1152 of a residual :math:`|\mathbf{u}^{p+1}-\mathbf{u}^{p}|` that, as
1153 :math:`\mathbf{u}^{p+1} \rightarrow \mathbf{u}^{n+1}`, converges to
1154 :math:`0`. In this way EVP can be re-interpreted as a pure iterative solver
1155 where the sub-cycling has no association with time-relation (through
dc26f158aa Mart*1156 :math:`\Delta{t}_{\mathrm{EVP}}`). With the setting of
1157 :varlink:`SEAICEuseEVPstar` to ``.TRUE.`` (default), this form of EVP is used.
1158 Using the terminology of Kimmritz et al. 2015 :cite:`kimmritz:15`, the evolution
1159 equations of stress :math:`\sigma_{ij}` and momentum :math:`\mathbf{u}` can be
1160 written as:
adc83e5d7b Mart*1161
1162 .. math::
258fe29c91 Jeff*1163 \sigma_{ij}^{p+1}=\sigma_{ij}^p+\frac{1}{\alpha}
1164 \Big(\sigma_{ij}(\mathbf{u}^p)-\sigma_{ij}^p\Big),
1165 \phantom{\int}
1166 :label: eq_evpstarsigma
0452697f42 Oliv*1167
1168 .. math::
258fe29c91 Jeff*1169 \mathbf{u}^{p+1}=\mathbf{u}^p+\frac{1}{\beta}
0bad585a21 Navi*1170 \Big(\frac{\Delta t}{m} \nabla \cdot\boldsymbol{\sigma}^{p+1}+
258fe29c91 Jeff*1171 \frac{\Delta t}{m}\mathbf{R}^{p}+\mathbf{u}_n
c61841e2fd Jeff*1172 -\mathbf{u}^p\Big)
258fe29c91 Jeff*1173 :label: eq_evpstarmom
adc83e5d7b Mart*1174
c512e371cc drin*1175 :math:`\mathbf{R}` contains all terms in the momentum equations except for the
1176 rheology terms and the time derivative; :math:`\alpha` and :math:`\beta` are
1177 free parameters (:varlink:`SEAICE_evpAlpha`, :varlink:`SEAICE_evpBeta`) that
1178 replace the time stepping parameters :varlink:`SEAICE_deltaTevp`
1179 (:math:`\Delta{t}_{\mathrm{EVP}}`), :varlink:`SEAICE_elasticParm`
1180 (:math:`E_{0}`), or :varlink:`SEAICE_evpTauRelax` (:math:`T`). :math:`\alpha`
1181 and :math:`\beta` determine the speed of convergence and the
1182 stability. Usually, it makes sense to use :math:`\alpha = \beta`, and
1183 :varlink:`SEAICEnEVPstarSteps` :math:`\gg (\alpha,\,\beta)` (Kimmritz et
1184 al. 2015 :cite:`kimmritz:15`). Currently, there is no termination criterion and
1185 the number of mEVP iterations is fixed to :varlink:`SEAICEnEVPstarSteps`.
adc83e5d7b Mart*1186
dc26f158aa Mart*1187 In order to use mEVP in MITgcm, compile with both ``#define``
1188 :varlink:`SEAICE_CGRID` and ``#define`` :varlink:`SEAICE_ALLOW_EVP` in
1189 :filelink:`SEAICE_OPTIONS.h <pkg/seaice/SEAICE_OPTIONS.h>` (default) and make
1190 sure that :varlink:`SEAICEuseEVPstar` ``= .TRUE.,`` (default) in ``data.seaice``.
1191 By default :varlink:`SEAICEuseEVPrev` is set to ``.TRUE.`` and the
1192 actual form of equations :eq:`eq_evpstarsigma` and :eq:`eq_evpstarmom` is used
1193 with fewer implicit terms and the factor of :math:`e^{2}` dropped in the stress
1194 equations :eq:`eq_evpstresstensor2` and :eq:`eq_evpstresstensor12`. Although
1195 this modifies the original EVP equations, it turns out to improve convergence
c512e371cc drin*1196 (Bouillon et al. 2013 :cite:`bouillon:13`).
adc83e5d7b Mart*1197
dc26f158aa Mart*1198 The aEVP scheme is an enhanced variant of mEVP (Kimmritz et al. 2016
1199 :cite:`kimmritz:16`), where the value of :math:`\alpha` is set dynamically based
1200 on the stability criterion
adc83e5d7b Mart*1201
1202 .. math::
c61841e2fd Jeff*1203 \alpha = \beta = \max\left( \tilde{c} \pi\sqrt{c \frac{\zeta}{A_{c}}
258fe29c91 Jeff*1204 \frac{\Delta{t}}{\max(m,10^{-4}\,\text{kg})}},\alpha_{\min} \right)
0452697f42 Oliv*1205 :label: eq_aevpalpha
adc83e5d7b Mart*1206
c512e371cc drin*1207 with the grid cell area :math:`A_c` and the ice and snow mass :math:`m`. This
1208 choice sacrifices speed of convergence for stability with the result that aEVP
1209 converges quickly to VP where :math:`\alpha` can be small and more slowly in
1210 areas where the equations are stiff. In practice, aEVP leads to an overall
dc26f158aa Mart*1211 better convergence than mEVP (Kimmritz et al. 2016 :cite:`kimmritz:16`). To use
1212 aEVP in MITgcm set :varlink:`SEAICEaEVPcoeff` :math:`= \tilde{c}`
1213 (see :eq:`eq_aevpalpha`; default is unset); this also
1214 sets the default values of :varlink:`SEAICEaEVPcStar` (:math:`c=4`) and
c512e371cc drin*1215 :varlink:`SEAICEaEVPalphaMin` (:math:`\alpha_{\min}=5`). Good convergence has
1216 been obtained with these values (Kimmritz et al. 2016 :cite:`kimmritz:16`):
adc83e5d7b Mart*1217
dc26f158aa Mart*1218 ::
1219
1220 SEAICEaEVPcoeff = 0.5,
1221 SEAICEnEVPstarSteps = 500,
1222 # The following two parameters are required by mEVP and aEVP,
1223 # but they are TRUE by default:
1224 SEAICEuseEVPstar = .TRUE.,
1225 SEAICEuseEVPrev = .TRUE.,
1226
1227 Because of the C-grid staggering of velocities and
c512e371cc drin*1228 stresses, mEVP may not converge as successfully as in Kimmritz et al. (2015)
dc26f158aa Mart*1229 :cite:`kimmritz:15`, see also Kimmritz et al. (2016) :cite:`kimmritz:16`.
1230 Convergence at very high resolution (order 5 km) has not yet been studied.
adc83e5d7b Mart*1231
1232 .. _para_phys_pkg_seaice_iceoceanstress:
1233
1234 Ice-Ocean stress
258fe29c91 Jeff*1235 ----------------
adc83e5d7b Mart*1236
c512e371cc drin*1237 Moving sea ice exerts a stress on the ocean which is the opposite of the stress
1238 :math:`\mathbf{\tau}_\mathrm{ocean}` in :eq:`eq_momseaice`. This stress is
1239 applied directly to the surface layer of the ocean model. An alternative ocean
1240 stress formulation is given by Hibler and Bryan (1987)
1241 :cite:`hibler:87`. Rather than applying :math:`\mathbf{\tau}_\mathrm{ocean}`
1242 directly, the stress is derived from integrating over the ice thickness to the
1243 bottom of the oceanic surface layer. In the resulting equation for the
1244 *combined* ocean-ice momentum, the interfacial stress cancels and the total
1245 stress appears as the sum of windstress and divergence of internal ice
1246 stresses: :math:`\delta(z) (\mathbf{\tau}_\mathrm{air} + \mathbf{F})/\rho_0`,
1247 see also Eq. (2) of Hibler and Bryan (1987) :cite:`hibler:87`. The disadvantage
1248 of this formulation is that now the velocity in the surface layer of the ocean
1249 that is used to advect tracers, is really an average over the ocean surface
adc83e5d7b Mart*1250 velocity and the ice velocity leading to an inconsistency as the ice
c512e371cc drin*1251 temperature and salinity are different from the oceanic variables. To turn on
1252 the stress formulation of Hibler and Bryan (1987) :cite:`hibler:87`, set
dc26f158aa Mart*1253 :varlink:`useHB87StressCoupling` ``=.TRUE.,``, in ``data.seaice``.
adc83e5d7b Mart*1254
1255 .. _para_phys_pkg_seaice_discretization:
1256
1257
1258 Finite-volume discretization of the stress tensor divergence
258fe29c91 Jeff*1259 ------------------------------------------------------------
adc83e5d7b Mart*1260
c512e371cc drin*1261 On an Arakawa C grid, ice thickness and concentration and thus ice strength
1262 :math:`P` and bulk and shear viscosities :math:`\zeta` and :math:`\eta` are
1263 naturally defined a C-points in the center of the grid cell. Discretization
1264 requires only averaging of :math:`\zeta` and :math:`\eta` to vorticity or
1265 Z-points (or :math:`\zeta`-points, but here we use Z in order avoid confusion
1266 with the bulk viscosity) at the bottom left corner of the cell to give
1267 :math:`\overline{\zeta}^{Z}` and :math:`\overline{\eta}^{Z}`. In the following,
1268 the superscripts indicate location at Z or C points, distance across the cell
1269 (F), along the cell edge (G), between :math:`u`-points (U), :math:`v`-points
1270 (V), and C-points (C). The control volumes of the :math:`u`- and
adc83e5d7b Mart*1271 :math:`v`-equations in the grid cell at indices :math:`(i,j)` are
1272 :math:`A_{i,j}^{w}` and :math:`A_{i,j}^{s}`, respectively. With these
1273 definitions (which follow the model code documentation except that
c512e371cc drin*1274 :math:`\zeta`-points have been renamed to Z-points), the strain rates are
1275 discretized as:
adc83e5d7b Mart*1276
1277 .. math::
1278 \begin{aligned}
1279 \dot{\epsilon}_{11} &= \partial_{1}{u}_{1} + k_{2}u_{2} \\ \notag
2c231b0ebd Mart*1280 => (\epsilon_{11})_{i,j}^C &= \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}}
1281 + k_{2,i,j}^{C}\frac{v_{i,j+1}+v_{i,j}}{2} \\
adc83e5d7b Mart*1282 \dot{\epsilon}_{22} &= \partial_{2}{u}_{2} + k_{1}u_{1} \\\notag
2c231b0ebd Mart*1283 => (\epsilon_{22})_{i,j}^C &= \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}}
1284 + k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\
adc83e5d7b Mart*1285 \dot{\epsilon}_{12} = \dot{\epsilon}_{21} &= \frac{1}{2}\biggl(
1286 \partial_{1}{u}_{2} + \partial_{2}{u}_{1} - k_{1}u_{2} - k_{2}u_{1}
1287 \biggr) \\ \notag
1288 => (\epsilon_{12})_{i,j}^Z &= \frac{1}{2}
2c231b0ebd Mart*1289 \biggl( \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^V}
adc83e5d7b Mart*1290 + \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^U} \\\notag
1291 &\phantom{=\frac{1}{2}\biggl(}
1292 - k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
1293 - k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
c512e371cc drin*1294 \biggr),
1295 \end{aligned}
adc83e5d7b Mart*1296
c512e371cc drin*1297 so that the diagonal terms of the strain rate tensor are naturally defined at
1298 C-points and the symmetric off-diagonal term at Z-points. No-slip boundary
1299 conditions (:math:`u_{i,j-1}+u_{i,j}=0` and :math:`v_{i-1,j}+v_{i,j}=0` across
1300 boundaries) are implemented via “ghost-pointsâ€; for free slip boundary
1301 conditions :math:`(\epsilon_{12})^Z=0` on boundaries.
adc83e5d7b Mart*1302
1303 For a spherical polar grid, the coefficients of the metric terms are
1304 :math:`k_{1}=0` and :math:`k_{2}=-\tan\phi/a`, with the spherical radius
c512e371cc drin*1305 :math:`a` and the latitude :math:`\phi`; :math:`\Delta{x}_1 = \Delta{x} =
1306 a\cos\phi \Delta\lambda`, and :math:`\Delta{x}_2 = \Delta{y}=a\Delta\phi`. For
1307 a general orthogonal curvilinear grid, :math:`k_{1}` and :math:`k_{2}` can be
1308 approximated by finite differences of the cell widths:
adc83e5d7b Mart*1309
1310 .. math::
1311 \begin{aligned}
1312 k_{1,i,j}^{C} &= \frac{1}{\Delta{y}_{i,j}^{F}}
1313 \frac{\Delta{y}_{i+1,j}^{G}-\Delta{y}_{i,j}^{G}}{\Delta{x}_{i,j}^{F}} \\
1314 k_{2,i,j}^{C} &= \frac{1}{\Delta{x}_{i,j}^{F}}
1315 \frac{\Delta{x}_{i,j+1}^{G}-\Delta{x}_{i,j}^{G}}{\Delta{y}_{i,j}^{F}} \\
1316 k_{1,i,j}^{Z} &= \frac{1}{\Delta{y}_{i,j}^{U}}
1317 \frac{\Delta{y}_{i,j}^{C}-\Delta{y}_{i-1,j}^{C}}{\Delta{x}_{i,j}^{V}} \\
1318 k_{2,i,j}^{Z} &= \frac{1}{\Delta{x}_{i,j}^{V}}
c512e371cc drin*1319 \frac{\Delta{x}_{i,j}^{C}-\Delta{x}_{i,j-1}^{C}}{\Delta{y}_{i,j}^{U}}
1320 \end{aligned}
adc83e5d7b Mart*1321
1322 The stress tensor is given by the constitutive viscous-plastic relation
1323 :math:`\sigma_{\alpha\beta} = 2\eta\dot{\epsilon}_{\alpha\beta} +
c512e371cc drin*1324 [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2 ]\delta_{\alpha\beta}` . The
1325 stress tensor divergence :math:`(\nabla\sigma)_{\alpha} =
1326 \partial_\beta\sigma_{\beta\alpha}`, is discretized in finite volumes . This
1327 conveniently avoids dealing with further metric terms, as these are “hidden†in
1328 the differential cell widths. For the :math:`u`-equation (:math:`\alpha=1`) we
1329 have:
adc83e5d7b Mart*1330
1331 .. math::
1332 \begin{aligned}
1333 (\nabla\sigma)_{1}: \phantom{=}&
1334 \frac{1}{A_{i,j}^w}
c512e371cc drin*1335 \int_{\mathrm{cell}}(\partial_1\sigma_{11}+\partial_2\sigma_{21})
1336 \,dx_1\,dx_2 \\\notag
adc83e5d7b Mart*1337 =& \frac{1}{A_{i,j}^w} \biggl\{
c512e371cc drin*1338 \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{11}dx_2\biggl|_{x_{1}}^{x_{1}
1339 +\Delta{x}_{1}}
1340 + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{21}dx_1\biggl|_{x_{2}}^{x_{2}
1341 +\Delta{x}_{2}}
adc83e5d7b Mart*1342 \biggr\} \\ \notag
1343 \approx& \frac{1}{A_{i,j}^w} \biggl\{
1344 \Delta{x}_2\sigma_{11}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
1345 + \Delta{x}_1\sigma_{21}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
1346 \biggr\} \\ \notag
1347 =& \frac{1}{A_{i,j}^w} \biggl\{
1348 (\Delta{x}_2\sigma_{11})_{i,j}^C -
2c231b0ebd Mart*1349 (\Delta{x}_2\sigma_{11})_{i-1,j}^C
adc83e5d7b Mart*1350 \\\notag
1351 \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{}
1352 + (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z
c512e371cc drin*1353 \biggr\}
1354 \end{aligned}
adc83e5d7b Mart*1355
1356 with
1357
1358 .. math::
1359 \begin{aligned}
1360 (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+}
1361 \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
1362 \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
1363 &+ \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
1364 k_{2,i,j}^C \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
1365 \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
1366 \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
1367 \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
1368 k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
1369 \phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\
1370 (\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+}
1371 \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
1372 \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag
1373 & + \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
1374 \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag
2c231b0ebd Mart*1375 & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
adc83e5d7b Mart*1376 k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag
2c231b0ebd Mart*1377 & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
c512e371cc drin*1378 k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
1379 \end{aligned}
adc83e5d7b Mart*1380
1381 Similarly, we have for the :math:`v`-equation (:math:`\alpha=2`):
1382
1383 .. math::
1384 \begin{aligned}
1385 (\nabla\sigma)_{2}: \phantom{=}&
1386 \frac{1}{A_{i,j}^s}
c512e371cc drin*1387 \int_{\mathrm{cell}}(\partial_1\sigma_{12}+\partial_2\sigma_{22})
1388 \,dx_1\,dx_2 \\\notag
adc83e5d7b Mart*1389 =& \frac{1}{A_{i,j}^s} \biggl\{
c512e371cc drin*1390 \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{12}dx_2\biggl|_{x_{1}}^{x_{1}
1391 +\Delta{x}_{1}}
1392 + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{22}dx_1\biggl|_{x_{2}}^{x_{2}
1393 +\Delta{x}_{2}}
adc83e5d7b Mart*1394 \biggr\} \\ \notag
1395 \approx& \frac{1}{A_{i,j}^s} \biggl\{
1396 \Delta{x}_2\sigma_{12}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
1397 + \Delta{x}_1\sigma_{22}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
1398 \biggr\} \\ \notag
1399 =& \frac{1}{A_{i,j}^s} \biggl\{
1400 (\Delta{x}_2\sigma_{12})_{i+1,j}^Z - (\Delta{x}_2\sigma_{12})_{i,j}^Z
1401 \\ \notag
1402 \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{}
1403 + (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C
1404 \biggr\} \end{aligned}
1405
1406 with
1407
1408 .. math::
1409 \begin{aligned}
1410 (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+}
1411 \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
2c231b0ebd Mart*1412 \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}}
adc83e5d7b Mart*1413 \\\notag &
1414 + \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
1415 \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\\notag
1416 &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
2c231b0ebd Mart*1417 k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
adc83e5d7b Mart*1418 \\\notag &
1419 - \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
1420 k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag
1421 (\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+}
1422 \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
1423 \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
1424 &+ \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
1425 k_{2,i,j}^{C} \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
1426 & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
1427 \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
1428 & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
1429 k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
1430 & -\Delta{x}_{i,j}^{F} \frac{P}{2}\end{aligned}
1431
258fe29c91 Jeff*1432 Again, no-slip boundary conditions are realized via ghost points and
adc83e5d7b Mart*1433 :math:`u_{i,j-1}+u_{i,j}=0` and :math:`v_{i-1,j}+v_{i,j}=0` across
c512e371cc drin*1434 boundaries. For free-slip boundary conditions the lateral stress is set to
1435 zeros. In analogy to :math:`(\epsilon_{12})^Z=0` on boundaries, we set
1436 :math:`\sigma_{21}^{Z}=0`, or equivalently :math:`\eta_{i,j}^{Z}=0`, on
1437 boundaries.
adc83e5d7b Mart*1438
a4e168e012 antn*1439 .. _ssub_phys_pkg_seaice_thermodynamics:
adc83e5d7b Mart*1440
1441 Thermodynamics
a4e168e012 antn*1442 ==============
adc83e5d7b Mart*1443
c61841e2fd Jeff*1444 **NOTE: THIS SECTION IS STILL NOT COMPLETE**
adc83e5d7b Mart*1445
c512e371cc drin*1446 In its original formulation the sea ice model uses simple 0-layer
1447 thermodynamics following the appendix of Semtner (1976)
1448 :cite:`semtner:76`. This formulation neglects storage of heat, that is, the
1449 heat capacity of ice is zero, and all internal heat sources so that the heat
1450 equation reduces to a constant conductive heat flux. This constant upward
1451 conductive heat flux together with a constant ice conductivity implies a linear
1452 temperature profile. The boundary conditions for the heat equations are: at the
0bad585a21 Navi*1453 bottom of the ice :math:`T|_{\rm bottom} = T_{\rm fr}` (freezing point temperature of
1454 sea water), and at the surface: :math:`Q_{\rm top} =
1455 \frac{\partial{T}}{\partial{z}} = (K/h)(T_{0}-T_{\rm fr})`, where :math:`K` is the
1456 ice conductivity, :math:`h` the ice thickness, and :math:`T_{0}-T_{\rm fr}` the
c512e371cc drin*1457 difference between the ice surface temperature and the water temperature at the
1458 bottom of the ice (at the freezing point). The surface heat flux
0bad585a21 Navi*1459 :math:`Q_{\rm top}` is computed in a similar way to that of Parkinson and
c512e371cc drin*1460 Washington (1979) :cite:`parkinson:79` and Manabe et al. (1979)
1461 :cite:`manabe:79`. The resulting equation for surface temperature is
adc83e5d7b Mart*1462
c61841e2fd Jeff*1463 .. math::
1464 \begin{aligned}
0bad585a21 Navi*1465 \frac{K}{h}(T_{0}-T_{\rm fr}) &= Q_{\rm SW\downarrow}(1-\mathrm{albedo}) \\
1466 & + \epsilon Q_{\rm LW\downarrow} - Q_{\rm LW\uparrow}(T_{0}) \\
1467 & + Q_{\rm LH}(T_{0}) + Q_{\rm SH}(T_{0}),
c61841e2fd Jeff*1468 \end{aligned}
1469 :label: eq_zerolayerheatbalance
2c231b0ebd Mart*1470
c61841e2fd Jeff*1471 where :math:`\epsilon` is the emissivity of the surface (snow or ice),
0bad585a21 Navi*1472 :math:`Q_{\rm S/LW\downarrow}` the downwelling shortwave and longwave radiation to
1473 be prescribed, and :math:`Q_{\rm LW\uparrow}=\epsilon\sigma_B T_{0}^4` the emitted
c512e371cc drin*1474 long wave radiation with the Stefan-Boltzmann constant :math:`\sigma_B`. With
1475 explicit expressions in :math:`T_0` for the turbulent fluxes of latent and
1476 sensible heat
c61841e2fd Jeff*1477
1478 .. math::
2c231b0ebd Mart*1479 \begin{aligned}
0bad585a21 Navi*1480 Q_{\rm LH} &= \rho_\mathrm{air} C_E (\Lambda_v + \Lambda_f)
2c231b0ebd Mart*1481 |\mathbf{U}_\mathrm{air}|
1482 \left[ q_\mathrm{air} - q_\mathrm{sat}(T_0)\right] \\
0bad585a21 Navi*1483 Q_{\rm SH} &= \rho_\mathrm{air} c_p C_E |\mathbf{U}_\mathrm{air}|
2c231b0ebd Mart*1484 \left[ T_\mathrm{10m} - T_{0} \right],
1485 \end{aligned}
c61841e2fd Jeff*1486
0bad585a21 Navi*1487 :eq:`eq_zerolayerheatbalance` can be solved for :math:`T_0` with an iterative
1488 Ralphson-Newton method, which usually converges very quickly in less that 10
1489 iterations. In these equations, :math:`\rho_\mathrm{air}` is the air density
1490 (parameter :varlink:`SEAICE_rhoAir`), :math:`C_E` is the ice-ocean transfer
1491 coefficient for sensible and latent heat (parameter :varlink:`SEAICE_dalton`),
c512e371cc drin*1492 :math:`\Lambda_v` and :math:`\Lambda_f` are the latent heat of vaporization and
0bad585a21 Navi*1493 fusion, respectively (parameters :varlink:`SEAICE_lhEvap` and
1494 :varlink:`SEAICE_lhFusion`), and :math:`c_p` is the specific heat of air
1495 (parameter :varlink:`SEAICE_cpAir`). For the latent heat :math:`Q_{\rm LH}` a
1496 choice can be made between the old polynomial expression for saturation
1497 humidity :math:`q_\mathrm{sat}(T_0)` (by setting
1498 :varlink:`useMaykutSatVapPoly` to ``.TRUE.``) and the default exponential
1499 relation approximation that is more accurate at low temperatures.
c512e371cc drin*1500
1501 In the zero-layer model of Semtner (1976) :cite:`semtner:76`, the conductive
1502 heat flux depends strongly on the ice thickness :math:`h`. However, the ice
1503 thickness in the model represents a mean over a potentially very heterogeneous
1504 thickness distribution. In order to parameterize a sub-grid scale distribution
14673ec2d0 Mart*1505 for heat flux computations, the ice thickness :math:`h` is split into
c512e371cc drin*1506 :math:`N` thickness categories :math:`H_{n}` that are equally distributed
1507 between :math:`2h` and a minimum imposed ice thickness of :math:`5\,\text{cm}`
1508 by :math:`H_n= \frac{2n-1}{7}\,h` for :math:`n\in[1,N]`. The heat fluxes
1509 computed for each thickness category are area-averaged to give the total heat
1510 flux (see Hibler 1984 :cite:`hibler:84`). To use this thickness category
1511 parameterization set :varlink:`SEAICE_multDim` to the number of desired
1512 categories in ``data.seaice`` (7 is a good guess, for anything larger than 7
1513 modify :filelink:`SEAICE_SIZE.h <pkg/seaice/SEAICE_SIZE.h>`). Note that this
1514 requires different restart files and switching this flag on in the middle of an
1515 integration is not advised. As an alternative to the flat distribution, the
1516 run-time parameter :varlink:`SEAICE_PDF` (1D-array of lenght :varlink:`nITD`)
1517 can be used to prescribe an arbitrary distribution of ice thicknesses, for
1518 example derived from observed distributions (Castro-Morales et al. 2014
1519 :cite:`castro-morales:14`). In order to include the ice thickness distribution
dc26f158aa Mart*1520 also for snow, set :varlink:`SEAICE_useMultDimSnow` to ``.TRUE.`` (this is the
c512e371cc drin*1521 default); only then, the parameterization of always having a fraction of thin
1522 ice is efficient and generally thicker ice is produced (see Castro-Morales et
1523 al. 2014 :cite:`castro-morales:14`).
1524
1525 The atmospheric heat flux is balanced by an oceanic heat flux from below. The
1526 oceanic flux is proportional to :math:`\rho\,c_{p}\left(T_{w}-T_{fr}\right)`
1527 where :math:`\rho` and :math:`c_{p}` are the density and heat capacity of sea
0bad585a21 Navi*1528 water and :math:`T_{\rm fr}` is the local freezing point temperature that is a
c512e371cc drin*1529 function of salinity. This flux is not assumed to instantaneously melt or
1530 create ice, but a time scale of three days (run-time parameter
1531 :varlink:`SEAICE_gamma_t`) is used to relax :math:`T_{w}` to the freezing
1532 point. The parameterization of lateral and vertical growth of sea ice follows
1533 that of Hibler (1979) and Hibler (1980) :cite:`hibler:79,hibler:80`; the
1534 so-called lead closing parameter :math:`h_{0}` (run-time parameter
1535 :varlink:`HO`) has a default value of 0.5 meters.
1536
1537 On top of the ice there is a layer of snow that modifies the heat flux and the
1538 albedo (Zhang et al. 1998 :cite:`zha:98`). Snow modifies the effective
1539 conductivity according to
adc83e5d7b Mart*1540
1541 .. math:: \frac{K}{h} \rightarrow \frac{1}{\frac{h_{s}}{K_{s}}+\frac{h}{K}},
1542
1543 where :math:`K_s` is the conductivity of snow and :math:`h_s` the snow
c512e371cc drin*1544 thickness. If enough snow accumulates so that its weight submerges the ice and
1545 the snow is flooded, a simple mass conserving parameterization of snowice
1546 formation (a flood-freeze algorithm following Archimedes’ principle) turns snow
1547 into ice until the ice surface is back at :math:`z=0` (see Leppäranta 1983
1548 :cite:`leppaeranta:83`). The flood-freeze algorithm is turned on with run-time
dc26f158aa Mart*1549 parameter :varlink:`SEAICEuseFlooding` set to ``.TRUE.``.
adc83e5d7b Mart*1550
1551 .. _para_phys_pkg_seaice_advection:
1552
1553 Advection of thermodynamic variables
258fe29c91 Jeff*1554 ------------------------------------
adc83e5d7b Mart*1555
14673ec2d0 Mart*1556 Mean ice thickness (ice volume per unit area, :math:`c h`, model variable
1557 :varlink:`HEFF`, which implies the misleading name "effective thickness"),
1558 concentration :math:`c` (model variable :varlink:`AREA`) and mean snow
1559 thickness (:math:`c h_s`, model variable :varlink:`HSNOW`) are advected by ice
1560 velocities:
adc83e5d7b Mart*1561
1562 .. math::
258fe29c91 Jeff*1563 \frac{\partial{X}}{\partial{t}} =
0bad585a21 Navi*1564 - \nabla \cdot\left(\mathbf{u}\,X\right) + \Gamma_{X} + D_{X}
0452697f42 Oliv*1565 :label: eq_advection
adc83e5d7b Mart*1566
c512e371cc drin*1567 where :math:`\Gamma_X` are the thermodynamic source terms and :math:`D_{X}` the
14673ec2d0 Mart*1568 diffusive terms for quantities :math:`X= c h, c, c h_s` or any other tracer,
1569 such as sea ice salinity. From the various advection schemes that are available
1570 in MITgcm, we recommend flux-limited schemes (runtime flag
1571 :varlink:`SEAICEadvScheme`; default=77, a 2nd-order flux limited scheme) to
1572 preserve sharp gradients and edges that are typical of sea ice distributions
1573 and to rule out unphysical over- and undershoots (negative thickness or
1574 concentration). These schemes conserve volume and horizontal area and are
1575 unconditionally stable, so that we can set :math:`D_{X}=0` (runtime flag
1576 :varlink:`DIFF1` = :math:`D_{X}/\Delta{x}`; default=0).
adc83e5d7b Mart*1577
c512e371cc drin*1578 The MITgcm sea ice model provides the option to use the thermodynamics model of
1579 Winton (2000) :cite:`winton:00`, which in turn is based on the 3-layer model of
1580 Semtner (1976) :cite:`semtner:76` which treats brine content by means of
1581 enthalpy conservation; the corresponding package :filelink:`thsice
1582 <pkg/thsice>` is described in section :numref:`sub_phys_pkg_thsice`. This
1583 scheme requires additional state variables, namely the enthalpy of the two ice
1584 layers (instead of effective ice salinity), to be advected by ice
1585 velocities. The internal sea ice temperature is inferred from ice enthalpy. To
1586 avoid unphysical (negative) values for ice thickness and concentration, a
258fe29c91 Jeff*1587 positive 2nd-order advection scheme with a SuperBee flux limiter (Roe 1985
c512e371cc drin*1588 :cite:`roe:85`) should be used to advect all sea-ice-related quantities of the
1589 Winton (2000) :cite:`winton:00` thermodynamic model (run-time flag
1590 :varlink:`thSIceAdvScheme` :math:`= 77` and :varlink:`thSIce_diffK` :math:`=
1591 D_{X} = 0` in ``data.ice``, defaults are 0). Because of the nonlinearity of the
1592 advection scheme, care must be taken in advecting these quantities: when simply
1593 using ice velocity to advect enthalpy, the total energy (i.e., the volume
1594 integral of enthalpy) is not conserved. Alternatively, one can advect the
1595 energy content (i.e., product of ice-volume and enthalpy) but then false
1596 enthalpy extrema can occur, which then leads to unrealistic ice temperature. In
1597 the currently implemented solution, the sea-ice mass flux is used to advect the
1598 enthalpy in order to ensure conservation of enthalpy and to prevent false
1599 enthalpy extrema.
adc83e5d7b Mart*1600
dce651c1fb Mart*1601 .. _para_phys_pkg_seaice_itd:
1602
1603 Dynamical Ice Thickness Distribution (ITD)
258fe29c91 Jeff*1604 ------------------------------------------
dce651c1fb Mart*1605
258fe29c91 Jeff*1606 The ice thickness distribution model used by MITgcm follows the implementation
1607 in the Los Alamos sea ice model CICE (https://github.com/CICE-Consortium/CICE).
c512e371cc drin*1608 There are two parts to it that are closely connected: the participation and
1609 ridging functions that determine which thickness classes take part in ridging
1610 and which thickness classes receive ice during ridging based on Thorndike et
1611 al. (1975) :cite:`thorndike:75`, and the ice strength parameterization by
1612 Rothrock (1975) :cite:`rothrock:75` which uses this information. The following
1613 description is slightly modified from Ungermann et al. (2017)
1614 :cite:`ungermann:17`. Verification experiment :filelink:`seaice_itd
1615 <verification/seaice_itd>` uses the ITD model.
dce651c1fb Mart*1616
1617 Distribution, participation and redistribution functions in ridging
258fe29c91 Jeff*1618 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
dce651c1fb Mart*1619
c512e371cc drin*1620 When :varlink:`SEAICE_ITD` is defined in :filelink:`SEAICE_OPTIONS.h
1621 <pkg/seaice/SEAICE_OPTIONS.h>`, the ice thickness is described by the ice
1622 thickness distribution :math:`g(h,\mathbf{x},t)` for the subgrid-scale (see
1623 Thorndike et al. 1975 :cite:`thorndike:75`), a probability density function for
1624 thickness :math:`h` following the evolution equation
dce651c1fb Mart*1625
1626
1627 .. math::
0bad585a21 Navi*1628 \frac{\partial g}{\partial t} = - \nabla \cdot (\mathbf{u} g) - \frac{\partial}{\partial h}(fg) + \Psi.
dce651c1fb Mart*1629 :label: eq_itd
1630
1631
c512e371cc drin*1632 Here :math:`f=\frac{\mathrm{d} h}{\mathrm{d} t}` is the thermodynamic growth
1633 rate and :math:`\Psi` a function describing the mechanical redistribution of
1634 sea ice during ridging or lead opening.
dce651c1fb Mart*1635
c512e371cc drin*1636 The mechanical redistribution function :math:`\Psi` generates open water in
1637 divergent motion and creates ridged ice during convergent motion. The ridging
1638 process depends on total strain rate and on the ratio between shear (run-time
1639 parameter :varlink:`SEAICEshearParm`) and divergent strain. In the single
1640 category model, ridge formation is treated implicitly by limiting the ice
1641 concentration to a maximum of one (see Hibler 1979 :cite:`hibler:79`), so that
1642 further volume increase in convergent motion leads to thicker ice. (This is
1643 also the default for ITD models; to change from the default, set run-time
dc26f158aa Mart*1644 parameter :varlink:`SEAICEsimpleRidging` ``=.FALSE.,`` in ``data.seaice``). For
c512e371cc drin*1645 the ITD model, the ridging mode in convergence
dce651c1fb Mart*1646
1647 .. math::
1648 \omega_r(h)= \frac{-a(h)+n(h)}{N}
2c231b0ebd Mart*1649
c512e371cc drin*1650 gives the effective change for the ice volume with thickness between :math:`h`
1651 and :math:`h+\textrm{d} h` as the normalized difference between the ice
1652 :math:`n(h)` generated by ridging and the ice :math:`a(h)` participating in
1653 ridging.
1654
1655 The participation function :math:`a(h) = b(h)g(h)` can be computed either
1656 following Thorndike et al. (1975) :cite:`thorndike:75` (run-time parameter
1657 :varlink:`SEAICEpartFunc` =0) or Lipscomb et al. (2007) :cite:`lipscomb:07`
1658 (:varlink:`SEAICEpartFunc` =1), and similarly the ridging function :math:`n(h)`
1659 can be computed following Hilber (1980) :cite:`hibler:80` (run-time parameter
1660 :varlink:`SEAICEredistFunc` =0) or Lipscomb et al. (2007) :cite:`lipscomb:07`
1661 (:varlink:`SEAICEredistFunc` =1). As an example, we show here the functions
1662 that Lipscomb et al. (2007) :cite:`lipscomb:07` suggested to avoid noise in the
1663 solutions. These functions are smooth and avoid non-differentiable
1664 discontinuities, but so far we did not find any noise issues as in Lipscomb et
1665 al. (2007) :cite:`lipscomb:07`.
1666
dc26f158aa Mart*1667 With :varlink:`SEAICEpartFunc` ``= 1,`` in ``data.seaice``, the participation
c512e371cc drin*1668 function with the relative amount of ice of thickness :math:`h` weighted by an
1669 exponential function
dce651c1fb Mart*1670
1671 .. math::
1672 b(h) = b_0 \exp [ -G(h)/a^*]
2c231b0ebd Mart*1673
c512e371cc drin*1674 where :math:`G(h)=\int_0^h g(h) \textrm{d} h` is the cumulative thickness
1675 distribution function, :math:`b_0` a normalization factor, and :math:`a^*`
1676 (:varlink:`SEAICEaStar`) the exponential constant that determines which
1677 relative amount of thicker and thinner ice take part in ridging.
dce651c1fb Mart*1678
dc26f158aa Mart*1679 With :varlink:`SEAICEredistFunc` ``= 1,`` in ``data.seaice``, the ice generated by
c512e371cc drin*1680 ridging is calculated as
dce651c1fb Mart*1681
1682 .. math::
1683 n(h) = \int_0^\infty a(h_1)\gamma(h_1,h) \textrm{d} h_1
1684
c512e371cc drin*1685 where the density function :math:`\gamma(h_1,h)` of resulting thickness
1686 :math:`h` for ridged ice with an original thickness of :math:`h_1` is taken as
dce651c1fb Mart*1687
1688 .. math::
c512e371cc drin*1689 \gamma(h_1, h) = \frac{1}{k \lambda}
1690 \exp\left[{\frac{-(h-h_{\min})}{\lambda}}\right]
1691
1692 for :math:`h \geq h_{\min}`, with :math:`\gamma(h_1,h)=0` for :math:`h <
1693 h_{\min}`. In this parameterization, the normalization factor
1694 :math:`k=\frac{h_{\min} + \lambda}{h_1}`, the e-folding scale :math:`\lambda =
1695 \mu h_1^{1/2}` and the minimum ridge thickness :math:`h_{\min}=\min(2h_1,h_1 +
1696 h_{\textrm{raft}})` all depend on the original thickness :math:`h_1`. The
1697 maximal ice thickness allowed to raft :math:`h_{\textrm{raft}}` is constant
1698 (:varlink:`SEAICEmaxRaft`, default =1 m) and :math:`\mu`
1699 (:varlink:`SEAICEmuRidging`) is a tunable parameter.
258fe29c91 Jeff*1700
1701 In the numerical model these equations are discretized into a set of :math:`n`
c512e371cc drin*1702 (:varlink:`nITD` defined in :filelink:`SEAICE_SIZE.h
1703 <pkg/seaice/SEAICE_SIZE.h>`) thickness categories employing the delta function
1704 scheme of Bitz et al. (2001) :cite:`bitz:01`. For each thickness category in
1705 an ITD configuration, the volume conservation equation :eq:`eq_advection` is
1706 evaluated using the heat flux with the category-specific values for ice and
1707 snow thickness, so there are no conceptual differences in the thermodynamics
1708 between the single category and ITD configurations. The only difference is
1709 that only in the thinnest category the creation of new ice of thickness
1710 :math:`H_0` (run-time parameter :varlink:`HO`) is possible, all other
1711 categories are limited to basal growth. The conservation of ice area is
1712 replaced by the evolution equation of the ITD :eq:`eq_itd` that is discretized
1713 in thickness space with :math:`n+1` category limits given by run-time parameter
1714 :varlink:`Hlimit`. If :varlink:`Hlimit` is not set in ``data.seaice``, a
1715 simple recursive formula following Lipscomb (2001) :cite:`lipscomb:01` is used
1716 to compute :varlink:`Hlimit`:
258fe29c91 Jeff*1717
1718 .. math::
c512e371cc drin*1719 H_\mathrm{limit}(k) = H_\mathrm{limit}(k-1) + \frac{c_1}{n}
1720 + \frac{c_1 c_2}{n} [ 1 + \tanh c_3 (\frac{k-1}{n} - 1) ]
258fe29c91 Jeff*1721
0bad585a21 Navi*1722 with :math:`H_\mathrm{limit}(0)=0` m and
1723 :math:`H_\mathrm{limit}(n)=999.9` m. The three constants are the
c512e371cc drin*1724 run-time parameters :varlink:`Hlimit_c1`, :varlink:`Hlimit_c2`, and
1725 :varlink:`Hlimit_c3`. The total ice concentration and volume can then be
1726 calculated by summing up the values for each category.
dce651c1fb Mart*1727
1728 Ice strength parameterization
258fe29c91 Jeff*1729 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
dce651c1fb Mart*1730
c512e371cc drin*1731 In the default approach of equation :eq:`eq_icestrength`, the ice strength is
1732 parameterized following Hibler (1979) :cite:`hibler:79` and :math:`P` depends
1733 only on average ice concentration and thickness per grid cell and the constant
1734 ice strength parameters :math:`P^{\ast}` (:varlink:`SEAICE_strength`) and
1735 :math:`C^{\ast}` (:varlink:`SEAICE_cStar`). With an ice thickness
1736 distribution, it is possible to use a different parameterization following
1737 Rothrock (1975) :cite:`rothrock:75`
dce651c1fb Mart*1738
1739 .. math::
258fe29c91 Jeff*1740 P = C_f C_p \int_0^\infty h^2 \omega_r(h) \textrm{d}h
dce651c1fb Mart*1741 :label: eq_rothrock
1742
c512e371cc drin*1743 by considering the production of potential energy and the frictional energy
1744 loss in ridging. The physical constant :math:`C_p = \rho_i (\rho_w - \rho_i)
1745 \hat{g} / (2 \rho_w)` is a combination of the gravitational acceleration
1746 :math:`\hat{g}` and the densities :math:`\rho_i`, :math:`\rho_w` of ice and
1747 water, and :math:`C_f` (:varlink:`SEAICE_cf`) is a scaling factor relating the
1748 amount of work against gravity necessary for ridging to the amount of work
1749 against friction. To calculate the integral, this parameterization needs
1750 information about the ITD in each grid cell, while the default
1751 parameterization :eq:`eq_icestrength` can be used for both ITD and single
1752 thickness category models. In contrast to :eq:`eq_icestrength`, which is based
1753 on the plausible assumption that thick and compact ice is stronger than thin
1754 and loose drifting ice, this parameterization :eq:`eq_rothrock` clearly
1755 contains the more physical assumptions about energy conservation. For that
1756 reason alone this parameterization is often considered to be more physically
1757 realistic than :eq:`eq_icestrength`, but in practice, the success is not so
1758 clear (Ungermann et al. 2007 :cite:`ungermann:17`). Ergo, the default is to
dc26f158aa Mart*1759 use :eq:`eq_icestrength`; set :varlink:`useHibler79IceStrength` ``=.FALSE.,`` in
c512e371cc drin*1760 ``data.seaice`` to change this behavior.
dce651c1fb Mart*1761
5b18319545 Jeff*1762 Known issues and work-arounds
1763 =============================
1764
1765 - An often encountered problem in long simulations with sea ice models is
1766 (local) perpetually increasing sea ice (plus snow) height; this is
1767 problematic when using a non-linear free surface and
dc26f158aa Mart*1768 :varlink:`useRealFreshWaterFlux` set to ``.TRUE.``, because the mass of the sea ice
5b18319545 Jeff*1769 places a load on the sea surface, which if too large, can cause the surface
1770 cells of the model to become too thin so that the model eventually stops with
1771 an error message. Usually this problem occurs because of dynamical ice growth
1772 (i.e., convergence and ridging of ice) or simply too much net precipitation
1773 with insufficient summer surface melting. If the problem is dynamical in
1774 nature (e.g., caused by ridging in a deep inlet), the first step to try is to
1775 turn off the replacement pressure method (:varlink:`SEAICEpressReplFac` = 0;
1776 in :numref:`para_phys_pkg_seaice_VPrheology`); turning this off provides
1777 resistance against additional growth due to further ridging, because the ice
1778 pressure :math:`P` is no longer reduced as :math:`\Delta\rightarrow 0` in
1779 nearly motionless thick ice :eq:`eq_pressrepl`. If this does not solve the
1780 problem, a somewhat more radical yet effective approach is simply to cap the
1781 sea ice load on the free surface by defining the CPP option
1782 :varlink:`SEAICE_CAP_ICELOAD`. This option effectively limits the sea ice
1783 load (variable :varlink:`sIceLoad`) to a mass of 1/5 of the the top grid cell
1784 depth. If desired, this limit can be changed in routine
1785 :filelink:`seaice_growth.F <pkg/seaice/seaice_growth.F>` where variable
1786 :varlink:`heffTooHeavy` is assigned.
dc26f158aa Mart*1787
61f2157921 Oliv*1788 .. _ssub_phys_pkg_seaice_subroutines:
adc83e5d7b Mart*1789
1790 Key subroutines
258fe29c91 Jeff*1791 ===============
adc83e5d7b Mart*1792
258fe29c91 Jeff*1793 Top-level routine: :filelink:`pkg/seaice/seaice_model.F`
adc83e5d7b Mart*1794
1795 ::
1796
1797
1798 C !CALLING SEQUENCE:
1799 c ...
1800 c seaice_model (TOP LEVEL ROUTINE)
1801 c |
1802 c |-- #ifdef SEAICE_CGRID
1803 c | SEAICE_DYNSOLVER
1804 c | |
1805 c | |-- < compute proxy for geostrophic velocity >
1806 c | |
1807 c | |-- < set up mass per unit area and Coriolis terms >
1808 c | |
1809 c | |-- < dynamic masking of areas with no ice >
1810 c | |
1811 c | |
1812 c | #ELSE
1813 c | DYNSOLVER
1814 c | #ENDIF
1815 c |
2c231b0ebd Mart*1816 c |-- if ( useOBCS )
adc83e5d7b Mart*1817 c | OBCS_APPLY_UVICE
1818 c |
1819 c |-- if ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow .OR. SEAICEadvSalt )
1820 c | SEAICE_ADVDIFF
1821 c |
1822 c | SEAICE_REG_RIDGE
1823 c |
2c231b0ebd Mart*1824 c |-- if ( usePW79thermodynamics )
adc83e5d7b Mart*1825 c | SEAICE_GROWTH
1826 c |
2c231b0ebd Mart*1827 c |-- if ( useOBCS )
adc83e5d7b Mart*1828 c | if ( SEAICEadvHeff ) OBCS_APPLY_HEFF
1829 c | if ( SEAICEadvArea ) OBCS_APPLY_AREA
1830 c | if ( SEAICEadvSALT ) OBCS_APPLY_HSALT
1831 c | if ( SEAICEadvSNOW ) OBCS_APPLY_HSNOW
1832 c |
1833 c |-- < do various exchanges >
1834 c |
1835 c |-- < do additional diagnostics >
1836 c |
1837 c o
1838
61f2157921 Oliv*1839 .. _ssub_phys_pkg_seaice_diagnostics:
adc83e5d7b Mart*1840
1841 SEAICE diagnostics
258fe29c91 Jeff*1842 ==================
adc83e5d7b Mart*1843
c512e371cc drin*1844 Diagnostics output is available via the diagnostics package (see
1845 :numref:`sub_outp_pkg_diagnostics`). Available output fields are summarized in
1846 the following table:
d25560575e Oliv*1847
1848 .. code-block:: text
1849
1850 ---------+----------+----------------+-----------------
1851 <-Name->|<- grid ->|<-- Units -->|<- Tile (max=80c)
1852 ---------+----------+----------------+-----------------
1853 sIceLoad|SM U1|kg/m^2 |sea-ice loading (in Mass of ice+snow / area unit)
1854 ---
1855 SEA ICE STATE:
1856 ---
1857 SIarea |SM M1|m^2/m^2 |SEAICE fractional ice-covered area [0 to 1]
1858 SIheff |SM M1|m |SEAICE effective ice thickness
1859 SIhsnow |SM M1|m |SEAICE effective snow thickness
1860 SIhsalt |SM M1|g/m^2 |SEAICE effective salinity
1861 SIuice |UU M1|m/s |SEAICE zonal ice velocity, >0 from West to East
1862 SIvice |VV M1|m/s |SEAICE merid. ice velocity, >0 from South to North
1863 ---
1864 ATMOSPHERIC STATE AS SEEN BY SEA ICE:
1865 ---
1866 SItices |SM C M1|K |Surface Temperature over Sea-Ice (area weighted)
1867 SIuwind |UM U1|m/s |SEAICE zonal 10-m wind speed, >0 increases uVel
1868 SIvwind |VM U1|m/s |SEAICE meridional 10-m wind speed, >0 increases uVel
1869 SIsnPrcp|SM U1|kg/m^2/s |Snow precip. (+=dw) over Sea-Ice (area weighted)
1870 ---
1871 FLUXES ACROSS ICE-OCEAN INTERFACE (ATMOS to OCEAN FOR ICE-FREE REGIONS):
1872 ---
1873 SIfu |UU U1|N/m^2 |SEAICE zonal surface wind stress, >0 increases uVel
1874 SIfv |VV U1|N/m^2 |SEAICE merid. surface wind stress, >0 increases vVel
1875 SIqnet |SM U1|W/m^2 |Ocean surface heatflux, turb+rad, >0 decreases theta
1876 SIqsw |SM U1|W/m^2 |Ocean surface shortwave radiat., >0 decreases theta
1877 SIempmr |SM U1|kg/m^2/s |Ocean surface freshwater flux, > 0 increases salt
1878 SIqneto |SM U1|W/m^2 |Open Ocean Part of SIqnet, turb+rad, >0 decr theta
1879 SIqneti |SM U1|W/m^2 |Ice Covered Part of SIqnet, turb+rad, >0 decr theta
1880 ---
1881 FLUXES ACROSS ATMOSPHERE-ICE INTERFACE (ATMOS to OCEAN FOR ICE-FREE REGIONS):
1882 ---
1883 SIatmQnt|SM U1|W/m^2 |Net atmospheric heat flux, >0 decreases theta
1884 SIatmFW |SM U1|kg/m^2/s |Net freshwater flux from atmosphere & land (+=down)
1885 SIfwSubl|SM U1|kg/m^2/s |Freshwater flux of sublimated ice, >0 decreases ice
1886 ---
1887 THERMODYNAMIC DIAGNOSTICS:
1888 ---
1889 SIareaPR|SM M1|m^2/m^2 |SIarea preceeding ridging process
1890 SIareaPT|SM M1|m^2/m^2 |SIarea preceeding thermodynamic growth/melt
1891 SIheffPT|SM M1|m |SIheff preceeeding thermodynamic growth/melt
1892 SIhsnoPT|SM M1|m |SIhsnow preceeeding thermodynamic growth/melt
1893 SIaQbOCN|SM M1|m/s |Potential HEFF rate of change by ocean ice flux
1894 SIaQbATC|SM M1|m/s |Potential HEFF rate of change by atm flux over ice
1895 SIaQbATO|SM M1|m/s |Potential HEFF rate of change by open ocn atm flux
1896 SIdHbOCN|SM M1|m/s |HEFF rate of change by ocean ice flux
1897 SIdSbATC|SM M1|m/s |HSNOW rate of change by atm flux over sea ice
1898 SIdSbOCN|SM M1|m/s |HSNOW rate of change by ocean ice flux
1899 SIdHbATC|SM M1|m/s |HEFF rate of change by atm flux over sea ice
1900 SIdHbATO|SM M1|m/s |HEFF rate of change by open ocn atm flux
1901 SIdHbFLO|SM M1|m/s |HEFF rate of change by flooding snow
1902 SIdAbATO|SM M1|m^2/m^2/s |Potential AREA rate of change by open ocn atm flux
1903 SIdAbATC|SM M1|m^2/m^2/s |Potential AREA rate of change by atm flux over ice
1904 SIdAbOCN|SM M1|m^2/m^2/s |Potential AREA rate of change by ocean ice flux
1905 SIdA |SM M1|m^2/m^2/s |AREA rate of change (net)
1906 ---
1907 DYNAMIC/RHEOLOGY DIAGNOSTICS:
1908 ---
b8665dacca Mart*1909 SIpress |SM M1|N/m |SEAICE strength (with upper and lower limit)
1910 SIzeta |SM M1|kg/s |SEAICE nonlinear bulk viscosity
1911 SIeta |SM M1|kg/s |SEAICE nonlinear shear viscosity
1912 SIsig1 |SM M1|no units |SEAICE normalized principle stress, component one
1913 SIsig2 |SM M1|no units |SEAICE normalized principle stress, component two
1914 SIshear |SM M1|1/s |SEAICE shear deformation rate
1915 SIdelta |SM M1|1/s |SEAICE Delta deformation rate
1916 SItensil|SM M1|N/m |SEAICE maximal tensile strength
d25560575e Oliv*1917 ---
1918 ADVECTIVE/DIFFUSIVE FLUXES OF SEA ICE variables:
1919 ---
1920 ADVxHEFF|UU M1|m.m^2/s |Zonal Advective Flux of eff ice thickn
1921 ADVyHEFF|VV M1|m.m^2/s |Meridional Advective Flux of eff ice thickn
1922 SIuheff |UU M1|m^2/s |Zonal Transport of eff ice thickn (centered)
1923 SIvheff |VV M1|m^2/s |Meridional Transport of eff ice thickn (centered)
1924 DFxEHEFF|UU M1|m^2/s |Zonal Diffusive Flux of eff ice thickn
1925 DFyEHEFF|VV M1|m^2/s |Meridional Diffusive Flux of eff ice thickn
1926 ADVxAREA|UU M1|m^2/m^2.m^2/s |Zonal Advective Flux of fract area
1927 ADVyAREA|VV M1|m^2/m^2.m^2/s |Meridional Advective Flux of fract area
1928 DFxEAREA|UU M1|m^2/m^2.m^2/s |Zonal Diffusive Flux of fract area
1929 DFyEAREA|VV M1|m^2/m^2.m^2/s |Meridional Diffusive Flux of fract area
1930 ADVxSNOW|UU M1|m.m^2/s |Zonal Advective Flux of eff snow thickn
1931 ADVySNOW|VV M1|m.m^2/s |Meridional Advective Flux of eff snow thickn
1932 DFxESNOW|UU M1|m.m^2/s |Zonal Diffusive Flux of eff snow thickn
1933 DFyESNOW|VV M1|m.m^2/s |Meridional Diffusive Flux of eff snow thickn
ba0b047096 Mart*1934 ADVxSSLT|UU M1|(g/kg).m^2/s |Zonal Advective Flux of seaice salinity
1935 ADVySSLT|VV M1|(g/kg).m^2/s |Meridional Advective Flux of seaice salinity
1936 DFxESSLT|UU M1|(g/kg).m^2/s |Zonal Diffusive Flux of seaice salinity
1937 DFyESSLT|VV M1|(g/kg).m^2/s |Meridional Diffusive Flux of seaice salinity
d25560575e Oliv*1938
adc83e5d7b Mart*1939
1940 Experiments and tutorials that use seaice
258fe29c91 Jeff*1941 =========================================
1942
1943 - :filelink:`verification/lab_sea`: Labrador Sea experiment
1944 - :filelink:`verification/seaice_obcs`, based on :filelink:`lab_sea <verification/lab_sea>`
c512e371cc drin*1945 - :filelink:`verification/offline_exf_seaice`, idealized topography in a zonally re-entrant channel, tests solvers and rheologies
258fe29c91 Jeff*1946 - :filelink:`verification/seaice_itd`, based on :filelink:`offline_exf_seaice <verification/offline_exf_seaice>`, tests ice thickness distribution
1947 - :filelink:`verification/global_ocean.cs32x15`, global cubed-sphere-experiment with combinations of :filelink:`pkg/seaice` and :filelink:`pkg/thsice`
1948 - :filelink:`verification/1D_ocean_ice_column`, just thermodynamics
