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