Sea-Ice Model

Sea-Ice Model

The sea-ice component of ROMS is a combination of the elastic-viscous-plastic (EVP) rheology ( Hunke and Dukowicz, 1997, Hunke, 2001) and simple one-layer ice and snow thermodynamics with a molecular sublayer under the ice ( Mellor and Kantha, 1989). It is tightly coupled, having the same grid (Arakawa-C) and timestep as the ocean and sharing the same parallel coding structure for use with MPI or OpenMP.

Dynamics

The momentum equations describe the change in ice/snow velocity due to the combined effects of the Coriolis force, surface ocean tilt, air and water stress, and internal ice stress:

$${\begin{aligned}{\frac {d}{dt}}(Ah_{i}u)&=Ah_{i}fv-Ah_{i}g{\frac {\partial \zeta _{w}}{\partial x}}+{\frac {A}{\rho _{i}}}(\tau _{a}^{x}+\tau _{w}^{x})+{\frac {1}{\rho _{i}}}{\cal {F}}_{x}\\{\frac {d}{dt}}(Ah_{i}v)&=-Ah_{i}fu-Ah_{i}g{\frac {\partial \zeta _{w}}{\partial y}}+{\frac {A}{\rho _{i}}}(\tau _{a}^{y}+\tau _{w}^{y})+{\frac {1}{\rho _{i}}}{\cal {F}}_{y}\,.\end{aligned}}$$

In this model, we neglect the nonlinear advection terms as well as the curvilinear terms in the internal ice stress. Nonlinear formulae are used for both the ocean-ice and air-ice surface stress:

$${\begin{aligned}{\vec {\tau }}_{a}&=\rho _{a}C_{a}|{\vec {V}}_{10}|{\vec {V}}_{10}\\C_{a}&={1 \over 2}C_{d}\left[1-\cos(2\pi \min(h_{i}+.1,.5)\right]\\{\vec {\tau }}_{w}&=\rho _{w}C_{w}|{\vec {v}}_{w}-{\vec {v}}|({\vec {v}}_{w}-{\vec {v}})\,.\end{aligned}}$$

The force due to the internal ice stress is given by the divergence of the stress tensor $\sigma$. The rheology is given by the stress-strain relation of the medium. We would like to emulate the viscous-plastic rheology of Hibler (1979):

$${\begin{aligned}\sigma _{ij}&=2\eta {\dot {\epsilon }}_{ij}+(\zeta -\eta ){\dot {\epsilon }}_{kk}\delta _{ij}-{P \over 2}\delta _{ij}\\{\dot {\epsilon }}_{ij}&\equiv {1 \over 2}\left({\partial u_{i} \over \partial x_{j}}+{\partial u_{j} \over \partial x_{i}}\right)\end{aligned}}$$

where the nonlinear viscosities are given by

$${\begin{aligned}\zeta &={P \over 2\left[(\epsilon _{11}^{2}+\epsilon _{22}^{2})(1+1/e^{2})+4e^{-2}\epsilon _{12}^{2}+2\epsilon _{11}\epsilon _{22}(1-1/e^{2})\right]^{1/2}}\\\eta &={\zeta \over e^{2}}\,.\end{aligned}}$$

Hibler's ice strength is given by:

$$P=P^{*}Ah_{i}e^{-C(1-A)}$$

while Overland and Pease (1988): advocate a nonlinear strength:

$$P=P^{*}Ah_{i}^{2}e^{-C(1-A)}$$

in which $P^{*}$ now depends on the grid spacing. Both options are available in ROMS.

We would also like to have an explicit model that can be solved efficiently on parallel computers. The EVP rheology has a tunable coefficient $E$ (the Young's modulus) which can be chosen to make the elastic term small compared to the other terms. We rearrange the VP rheology:

$${1 \over 2\eta }\sigma _{ij}+{\eta -\zeta \over 4\eta \zeta }\sigma _{kk}\delta _{ij}+{P \over 4\zeta }\delta _{ij}={\dot {\epsilon }}_{ij}$$

then add the elastic term:

$${1 \over E}{\partial \sigma _{ij} \over \partial t}+{1 \over 2\eta }\sigma _{ij}+{\eta -\zeta \over 4\eta \zeta }\sigma _{kk}\delta _{ij}+{P \over 4\zeta }\delta _{ij}={\dot {\epsilon }}_{ij}$$

Much like the ocean model, the ice model has a split timestep. The internal ice stress term is updated on a shorter timestep so as to allow the elastic wave velocity to be resolved.

Once the new ice velocities are computed, the ice tracers can be advected using the MPDATA scheme ( Smolarkiewicz). The tracers in this case are the ice thickness, ice concentration, snow thickness, internal ice temperature, and surface melt ponds. The continuity equations describing the evolution of these parameters also include thermodynamic terms ($S_{h}$, $S_{s}$ and $S_{A}$), described below:

$${\begin{aligned}{\partial Ah_{i} \over \partial t}&=-{\partial (uAh_{i}) \over \partial x}-{\partial (vAh_{i}) \over \partial y}+S_{h}+{\cal {D}}_{h}\\{\partial Ah_{s} \over \partial t}&=-{\partial (uAh_{s}) \over \partial x}-{\partial (vAh_{s}) \over \partial y}+S_{s}+{\cal {D}}_{s}\\{\partial A \over \partial t}&=-{\partial (uA) \over \partial x}-{\partial (vA) \over \partial y}+S_{A}+{\cal {D}}_{A}\qquad \qquad 0\leq A\leq 1\,.\end{aligned}}$$

The first two equations represent the conservation of ice and snow. The third is discussed in some detail in Mellor and Kantha, 1989, but represents the advection of ice blocks in which no ridging occurs as long as there is any open water.

The ice model variables:

Name	Description
$A(x,y,t)$	ice concentration, $0.0\leq A\leq 1.0$
$C_{a}$	nonlinear air drag coefficient
$C_{d}=2.2\times 10^{-3}$	air drag coefficient
$C_{w}=10\times 10^{-3}$	water drag coefficient
${\cal {D}}_{h},{\cal {D}}_{s},{\cal {D}}_{A}$	diffusion terms
$\delta _{ij}$	Kronecker delta function
$E$	Young's modulus
$e=2$	eccentricity of the elliptical yield curve
$\epsilon _{ij}(x,y,t)$	strain rate tensor
$\eta (x,y,t)$	nonlinear shear viscosity
$f(x,y)$	Coriolis parameter
${\cal {F}}_{x},{\cal {F}}_{y}$	forces due to internal ice stress
$g$	gravity
$H$	Heaviside function
$h_{i}(x,y,t)$	ice thickness
$h_{o}$	ice cutoff thickness
$h_{s}(x,y,t)$	snow thickness on ice-covered fraction
$M$	ice mass, $\rho _{i}Ah_{i}$
$P$	ice strength
$P^{*},C$	ice strength parameters
$S_{h},S_{s},S_{A}$	thermodynamic terms
$\sigma _{ij}$	internal ice stress tensor
${\vec {\tau }}_{a}$	air stress
${\vec {\tau }}_{w}$	water stress
$u,v$	the ($x,y$) components of ice velocity ${\vec {v}}$
${\vec {V}}_{10},{\vec {v}}_{w}$	10 meter air and surface water velocities
$\zeta$	nonlinear bulk viscosity
$\zeta _{w}$	surface elevation of the underlying water

Note that Hibler's $h_{I}$ variable is equivalent to our $Ah_{i}$ combination - his $h_{I}$ is the average thickness over the whole gridbox while our $h_{i}$ is the average thickness over the ice-covered fraction of the gridbox.

Thermodynamics

The thermodynamics is based on calculating how much ice grows and melts on each of the surface, bottom, and sides of the ice floes, as well as frazil ice formation:

Once the ice tracers are advected, the ice concentration and thickness are timestepped according to the terms on the right:

$${\begin{aligned}{\frac {DAh_{i}}{Dt}}&={\frac {\rho _{o}}{\rho _{i}}}\left[A(W_{io}-W_{ai})+(1-A)W_{ao}+W_{fr}\right]\\{\frac {DA}{Dt}}&={\frac {\rho _{o}}{\rho _{i}h_{i}}}\left[\Phi (1-A)w_{ao}+(1-A)W_{fr}\right]\;\;\;\;\;0\leq A\leq 1\end{aligned}}$$

The term $Ah_{i}$ is the "effective thickness", a measure of the ice volume. Its evolution equation is simply quantifying the change in the amount of ice. The ice concentration equation is more interesting in that it provides the partitioning between ice melt/growth on the sides vs. on the top and bottom. The parameter $\Phi$ controls this and has differing values for ice melt and retreat. In principle, most of the ice growth is assumed to happen at the base of the ice while rather more of the melt happens on the sides of the ice due to warming of the water in the leads.

The heat fluxes through the ice are based on a simple one layer Semtner (1976) type model with snow on top. The temperature is assumed to be linear within the snow and within the ice. The ice contains brine pockets for a total ice salinity of 3.2 or the surface ocean salinity, whichever is less. The surface ocean temperature and salinity is half a $dz$ below the surface. The water right below the surface is assumed to be at the freezing temperature; a logarithmic boundary layer is computed having the temperature and salinity matched at freezing.

Similarly, there is a snow equation:

$${DAh_{s} \over Dt}=\left[A(W_{s}-W_{sm})\right]$$

where $W_{s}$ and $W_{sm}$ are the snowfall and snow melt rates, respectively, in units of equivalent water. Also in the model is the depth of the melt ponds in spring, $h_{mp}$ which can be up to 0.1 m, after which melting ice contributes to $W_{ro}$. The melt ponds are not part of the thermal conductivity computations. They could contribute to the albedo computation, but that has been largely commented out.

The ice model thermodynamic variables:

Name	Description
$\alpha _{w}=0.10$	shortwave albedo of water
$\alpha _{i}=0.60$	shortwave albedo of wet ice
$\alpha _{i}=0.65$	shortwave albedo of dry ice
$\alpha _{s}$	shortwave albedo of wet snow
$\alpha _{s}$	shortwave albedo of dry snow
$C_{k}$	snow correction factor
$C_{pi}$ = 2093 J kg$^{-1}$ K$^{-1}$	specific heat of ice
$C_{po}$ = 3990 J kg$^{-1}$ K$^{-1}$	specific heat of water
$\epsilon _{w}=0.97$	longwave emissivity of water
$\epsilon _{i}=0.97$	longwave emissivity of ice
$\epsilon _{s}=0.99$	longwave emissivity of snow
$E(T,r)$	enthalpy of the ice/brine system
$F_{T}\!\uparrow$	heat flux from the ocean into the ice
$H\!\downarrow$	sensible heat
$i_{w}$	fraction of the solar heating transmitted through a lead into the water below
$k_{i}$ = 2.04 W m$^{-1}$ K${^{-}1}$	thermal conductivity of ice
$k_{s}$ = 0.31 W m$^{-1}$ K${^{-}1}$	thermal conductivity of snow
$L_{i}$ = 302 MJ m$^{-3}$	latent heat of fusion of ice
$L_{s}$ = 110 MJ m$^{-3}$	latent heat of fusion of snow
$LE\!\downarrow$	latent heat
$LW\!\!\downarrow$	incoming longwave radiation
$m$ & $-0.054^{\circ }$C/PSU	coefficient in linear $T_{f}(S)=mS$ equation
$\Phi$	contribution to $A$ equation from freezing water
$Q_{ai}$	heat flux out of the snow/ice surface
$Q_{ao}$	heat flux out of the ocean surface
$Q_{i2}$	heat flux up out of the ice
$Q_{io}$	heat flux up into the ice
$Q_{s}$	heat flux up through the snow
$r$	brine fraction in ice
$\rho _{i}$ = 910 $m^{3}/kg$	density of ice
$S_{i}$ = 3.2 PSU	salinity of the ice
$SW\!\!\downarrow$	incoming shortwave radiation
$\sigma$ = $5.67\times 10^{-8}$ W m$^{-2}$ K$^{-4}$	Stefan-Boltzmann constant
$T_{0}$	temperature of the bottom of the ice
$T_{1}$	temperature of the interior of the ice
$T_{2}$	temperature at the upper surface of the ice
$T_{3}$	temperature at the upper surface of the snow
$T_{f}$	freezing temperature
$T_{{\rm {melt}}\_i}$ = $mS_{i}$	melting temperature of ice
$T_{{\rm {melt}}\_s}$ = 0$^{\circ }$ C	melting temperature of snow
$W_{ai}$	melt rate on the upper ice/snow surface
$W_{ao}$	freeze rate at the air/water interface
$W_{fr}$	rate of frazil ice growth
$W_{io}$	freeze rate at the ice/water interface
$W_{ro}$	rate of run-off of surface melt water
$W_{s}$	snowfall rate
$W_{sm}$	snow melt rate

The locations of the ice and snow temperatures and the heat fluxes:

The temperature profile is assumed to be linear between adjacent temperature points. The interior of the ice contains "brine pockets", leading to a prognostic equation for the temperature $T_{1}$.

The surface flux to the air is:

$$Q_{ai}=-H\!\downarrow -LE\!\downarrow -\epsilon _{s}LW\!\!\downarrow -(1-\alpha _{s})SW\!\!\downarrow +\epsilon _{s}\sigma (T_{3}+273)^{4}$$

The incoming shortwave and longwave radiations are assumed to come from an atmospheric model. The formulae for sensible heat, latent heat, and outgoing longwave radiation are the same as in Parkinson and Washington (1979) and are shown in Radiant Heat Fluxes. The sensible heat is a function of $T_{3}$, as is the heat flux through the snow $Q_{s}$. Setting $Q_{ai}=Q_{s}$, we can solve for $T_{3}$ by setting $T_{3}^{n+1}=T_{3}^{n}+\Delta T_{3}$ and linearizing in $\Delta T_{3}$. As in Parkinson and Washington, if $T_{3}$ is found to be above the melting temperature, it is set to $T_{\rm {melt}}$ and the extra energy goes into melting the snow or ice:

$${\begin{aligned}W_{ai}&={\frac {Q_{ai}-Q_{i2}}{\rho _{o}L_{3}}}\\L_{3}&\equiv \left[E(T_{3},1)-E(T_{1},R_{1})\right]\end{aligned}}$$

Note that $L_{3}=(1-r)L_{i}$ plus a small sensible heat correction. We can store water on the surface in melt pools to a fixed depth---everything extra melted at the surface is assumed to flow into the ocean as $W_{ro}$. The melt ponds however do not contribute to the heat flux computation.

Inside the ice there are brine pockets in which there is salt water at the in situ freezing temperature. It is assumed that the ice has a uniform overall salinity of $S_{i}$ and that the freezing temperature is a linear function of salinity. The brine fraction $r$ is given by

$$r={\frac {S_{i}m}{T_{1}}}$$

The enthalpy of the combined ice/brine system is given by

$$E(T,r)=r(L_{i}+C_{po}T)+(1-r)C_{pi}T$$

Substituting in for $r$ and differentiating gives:

$${\frac {\partial E}{\partial T}}=-{\frac {S_{i}mL_{i}}{T_{1}^{2}}}+C_{pi}$$

Inside the snow, we have

$$Q_{s}={\frac {k_{s}}{h_{s}}}(T_{2}-T_{3})$$

The heat conduction in the upper part of the ice layer is

$$Q_{I2}={\frac {2k_{i}}{h_{i}}}(T_{1}-T_{2})$$

These can be set equal to each other to solve for $T_{2}$

$$T_{2}={\frac {T_{3}+C_{k}T_{1}}{1+C_{k}}}$$

where

$$C_{k}\equiv {\frac {2k_{i}h_{s}}{h_{i}k_{s}}}\,.$$

Substituting into (\ref{qi2}), we get:

$$Q_{s}=Q_{I2}={\frac {2k_{i}}{h_{i}}}{\frac {(T_{1}-T_{3})}{(1+C_{k})}}$$

Note that in the absence of snow, $C_{k}$ becomes zero and we recover the formula for the no-snow case in which $T_{3}=T_{2}$.

At the bottom of the ice, we have

$$Q_{I0}={\frac {2k_{i}}{h_{i}}}(T_{0}-T_{1})$$

The difference between $Q_{I0}$ and $Q_{I2}$ goes into the enthalpy of the ice:

$$\rho _{i}h_{i}\left[{\frac {\partial E}{\partial t}}+{\vec {v}}\cdot \nabla E\right]=Q_{I0}-Q_{I2}$$

We can use the chain rule to obtain an equation for timestepping $T_{1}$:

$$\rho _{i}h_{i}{\frac {\partial E}{\partial T}}\left[{\frac {\partial T_{1}}{\partial t}}+{\vec {v}}\cdot \nabla T_{1}\right]=Q_{I0}-Q_{I2}$$

where

$${\begin{aligned}Q_{I0}-Q_{I2}&={\frac {2k_{i}}{h_{i}}}\left[(T_{0}-T_{1})-{\frac {(T_{1}-T_{3})}{1+C_{k}}}\right]\\&={\frac {2k_{i}}{h_{i}}}\left[(T_{0}+{\frac {T_{3}-(2+C_{k})T_{1}}{1+C_{k}}}\right]\end{aligned}}$$

Ocean surface boundary conditions

The ocean receives surface stresses from both the atmosphere and the ice, according to the ice concentration:

$${\begin{aligned}K_{m}{\frac {\partial u_{w}}{\partial z}}&={\frac {A}{\rho _{o}}}\tau _{io}^{x}+{\frac {1-A}{\rho _{o}}}\tau _{ao}^{x}\\K_{m}{\frac {\partial v_{w}}{\partial z}}&={\frac {A}{\rho _{o}}}\tau _{io}^{y}+{\frac {1-A}{\rho _{o}}}\tau _{ao}^{y}\end{aligned}}$$

where the relevant variables are in the following table:

Variable	Value	Definition
$b$	3.14	factor
${\dot {E}}$		evaporation
$k$	0.4	von Karman's constant
$K_{m}$		vertical viscosity of seawater
$\nu$	$1.8\times 10^{-6}m^{2}s^{-1}$	kinematic viscosity of seawater
${\dot {P}}$		precipitation
$Pr$	13.0	molecular Prandtl number
$P_{rt}$	0.85	turbulent Prandtl number
$S_{0}$		surface salinity
$S$		internal ocean salinity
$Sc$	2432	molecular Schmidt number
$\tau _{io}$		stress on the ocean from the ice
$\tau _{ao}$		stress on the ocean from the wind
$T_{0}$		surface ocean temperature
$T$		internal ocean temperature
$u_{\tau }$		friction velocity $|\tau _{io}|^{1/2}\rho _{o}^{-1/2}$
$z_{0}$		roughness parameter

surface salinity ($T_{0}=mS_{0}$) in the presense of ice. We also have $T$ and $S$ at the uppermost computed ocean point ${1 \over 2}dz$ below the surface. In order to solve for $T_{0}$ and $S_{0}$, we assume a Yaglom and Kader, 1974 logarithmic boundary layer. The upper ocean heat flux is:

$${\frac {F_{T}}{\rho _{o}C_{po}}}=-C_{T_{z}}(T_{0}-T)\;\;\;z\rightarrow 0$$

where

$${\begin{aligned}C_{T_{z}}&={\frac {u_{\tau }}{P_{rt}k^{-1}\ln(-z/z_{0})+B_{T}}}\\B_{T}&=b\left({\frac {z_{0}u_{\tau }}{\nu }}\right)^{1/2}Pr^{2/3}\end{aligned}}$$

Likewise, we have the following equation for the surface salt flux:

$$F_{S}=-C_{S_{z}}(S_{0}-S)\;\;\;z\rightarrow 0$$

where

$${\begin{aligned}C_{S_{z}}&={\frac {u_{\tau }}{P_{rt}k^{-1}\ln(-z/z_{0})+B_{S}}}\\B_{S}&=b\left({\frac {z_{0}u_{\tau }}{\nu }}\right)^{1/2}Sc^{2/3}\end{aligned}}$$

The ocean model receives the following heat and salt fluxes:

$${\begin{aligned}F_{T}&=AQ_{io}+(1-A)Q_{ao}-W_{o}L_{o}\\F_{S}&=(W_{o}-AW_{ro})(S_{i}-S_{0})+(1-A)S_{o}({\dot {P}}-{\dot {E}})\\W_{o}&\equiv AW_{io}+(1-A)W_{ao}\end{aligned}}$$

Mellor and Kantha describe solving simultaneously for the five unknowns $W_{o}$, $T_{0}$, $S_{0}$, $F_{T}$ and $F_{S}$. Instead, we use the old value of $T_{0}$ to find $W_{io}$ and therefore $W_{o}$. Using the new value of $W_{io}$, solve for a new value of $S_{0}$ and then find the new $T_{0}$ as the freezing temperature for that salinity:

$${\begin{aligned}W_{io}&={1 \over L_{o}}\left[{Q_{io} \over \rho _{o}}+C_{po}C_{T_{z}}(T_{o}-T)\right]\\S_{0}&={C_{S_{z}}S+(W_{ai}H(-W_{ai})-W_{io}))S_{i} \over C_{S_{z}}+W_{ai}H(-W_{ai})+W_{ro}-W_{io}}\end{aligned}}$$

Note that the term $W_{ai}H(-W_{ai})$ is the contribution of melting only, not refreezing of melt ponds.

Frazil ice formation

Following Steele et al. (1989), we check to see if any of the ocean temperatures are below freezing at the end of each timestep. If so, frazil ice is formed, changing the local temperature and salinity. The ice that forms is assumed to instantly float up to the surface and add to the ice layer there. We balance the mass, heat, and salt before and after the ice is formed:

$${\begin{aligned}m_{w_{1}}&=m_{w_{2}}+m_{i}\\m_{w_{1}}(C_{pw}T_{1}+L)&=m_{w_{2}}(C_{pw}T_{2}+L)+m_{i}C_{pi}T_{2}\\m_{w_{1}}S_{1}&=m_{w_{2}}S_{2}\,.\end{aligned}}$$

The variables are defined in this table:

Variable	Value	Definition
$C_{pi}$	1994 J kg$^{-1}$ K$^{-1}$	specific heat of ice
$C_{pw}$	3987 J kg$^{-1}$ K$^{-1}$	specific heat of water
$\gamma$	$m_{i}/m_{w_{2}}$	fraction of water that froze
$L$	3.16e5 J kg$^{-1}$	latent heat of fusion
$m_{i}$		mass of ice formed
$m_{w_{1}}$		mass of water before freezing
$m_{w_{2}}$		mass of water after freezing
$m$	$-0.0543$	constant in freezing equation
$n$	$7.59\times 10^{-4}$	constant in freezing equation
$S_{1}$		salinity before freezing
$S_{2}$		salinity after freezing
$T_{1}$		temperature before freezing
$T_{2}$		temperature after freezing

Defining $\gamma =m_{i}/m_{w_{2}}$ and dropping terms of order $\gamma ^{2}$ leads to:

$$T_{2}=T_{1}+\gamma \left[{\frac {L}{C_{pw}}}+T_{1}\left(1-{\frac {C_{pi}}{C_{pw}}}\right)\right]$$

$$S_{2}=S_{1}(1+\gamma )\,.$$

We also want the final temperature and salinity to be on the freezing line, which we approximate as:

$$T_{f}=mS+nz\,.$$

We can then solve for $\gamma$:

$$\gamma ={\frac {-T_{1}+mS_{1}+nz}{{L}{C_{pw}}+T_{1}\left(1-{\frac {C_{pi}}{C_{pw}}}\right)-mS_{1}}}\,.$$

The ocean is checked at each depth $k$ and at each timestep for supercooling. If the water is below freezing, the temperature and salinity are adjusted as in equations (\ref{t2eq}) and (\ref{s2eq}) and the ice above is thickened by the amount:

$$\Delta h=\gamma _{k}\Delta z_{k}{\frac {\rho _{w}}{\rho _{i}}}\,.$$