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

The "pressure gradient" term is also modeled as a term in the internal ice stress. This term represents the resistance which ice has to being compressed (ice strength) and is a function of ice thickness and concentration: $$P = P^* A h_i \exp[ -C (1-A)] H(-\nabla \cdot \vec{v})$$ The Heaviside function guarantees that the ice has no strength when the flow is divergent Gray and Killworth.

$$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 ice model variables:

Name	Description
$u,v$	Horizontal ice velocity components
$M$	Ice mass, $\rho_i A h_i$
$A$	Ice concentration, $0.0 \leq A \leq 1.0$
$h_i$	Ice thickness
$f$	Coriolis parameter
$g$	Gravity
$\zeta_w$	Surface elevation of the underlying water
$\tau$	Air and water stresses
${\cal F}_x, {\cal F}_y$	Forces due to internal ice stress
$\sigma_{ij}$	Internal ice stress tensor
$\epsilon_{ij}$	Strainrate tensor
$\delta_{ij}$	Kronecker delta function
$\zeta, \eta$	Nonlinear ice viscosities
$P$	Ice strength
$P^*, C$	Ice strength parameters
$E$	Young's modulus

The second major component of the model consists of continuity equations describing the evolution of the ice thickness characteristics. Three parameters are calculated: the ice thickness $h_i$, the snow thickness $h_s$, and the compactness, $A$, which is defined as the fraction of area covered by thick ice. 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. 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: $${\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.

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>