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