Sea-Ice Model

Sea-Ice Model

<wikitex> 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:

$$M {du \over dt} = M f v - M g {\partial \zeta_w \over \partial x} + \tau_a^x + \tau_w^x + {\cal F}_x$$ $$M {dv \over dt} = - M f u - M g {\partial \zeta_w \over \partial y} + \tau_a^y + \tau_w^y + {\cal F}_y$$

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

 $$\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}) .$$ 

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):

$$\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)$$

$$P = P^* A h_i e^{-C(1-A)}$$ where the nonlinear viscosities are given by $$\zeta = { P \over 2 \left[ (\epsilon^2_{11} +

  \epsilon^2_{22} ) ( 1 + 1/e^2 ) + 4 e^{-2} \epsilon^2_{12}
  + 2 \epsilon_{11} \epsilon_{22} ( 1 - 1/e^2 ) \right] ^{1/2} }$$

$$\eta = { \zeta \over e^2 }.$$

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$), which will be described below:

$${\partial A h_i \over \partial t} =

 - {\partial (u A h_i) \over \partial x} -
 {\partial (v A h_i) \over \partial y}
 + S_h + {\cal D}_h$$

$${\partial A h_s \over \partial t} =

 - {\partial (u A h_s) \over \partial x} -
 {\partial (v A h_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 .$$

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. An optional ridging term can be added ( Gray and Killworth, 1995): $${\partial A \over \partial t} =

 - {\partial (uA) \over \partial x} - {\partial (vA) \over \partial y}
 - A \alpha(A) \, \nabla \cdot \vec{v} \, H(-\nabla \cdot \vec{v})
 + S_A + {\cal D}_A \qquad \qquad 0 \leq A \leq 1 .$$

where $\alpha(A)$ is an arbitrary function such that $\alpha(0) = 0$, $\alpha(1) = 1$, and $0 \leq \alpha(A) \leq 1$. The ridging term leads to an increase in $h_i$ under convergent flow as would be produced by ridging. The function $\alpha(A)$ should be chosen so that it is near zero until the ice concentration is large enough that ridging is expected to occur, then should increase smoothly to one.

The ice model variables:

Name	Description
$A(x,y,t)$	ice concentration, $0.0 \leq A \leq 1.0$
$\alpha(A)$	ridging function
$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$	Coriolis parameter
${\cal F}_x, {\cal F}_y$	forces due to internal ice stress
$g$	gravity
$H$	Heaviside function
$h_i$	ice thickness
$h_o$	ice cutoff thickness
$h_s$	snow thickness on ice-covered fraction
$M$	ice mass, $\rho_i A h_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 $A h_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:

 [Figure with W_xx]

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

$${D A h_i \over D t} = {\rho_o \over \rho_i} \left[ A (W_{io} - W_{ai}) + (1-A) W_{ao} + W_{fr} \right]$$ $${D A \over D t} = {\rho_o A \over \rho_i h_i} \left[ \Phi (1-A) w_{ao} + (1-A) W_{fr} \right] \qquad \qquad 0 \leq A \leq 1$$

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 5. 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.

[Note: I have gobs more here...]

Frazil Ice

</wikitex>