Sea-Ice Model

From WikiROMS
Jump to: navigation, search
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.


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:

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:

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

where the nonlinear viscosities are given by

Hibler's ice strength is given by:

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

in which 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 (the Young's modulus) which can be chosen to make the elastic term small compared to the other terms. We rearrange the VP rheology:

then add the elastic term:

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 (, and ), described below:

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
ice concentration,
nonlinear air drag coefficient
air drag coefficient
water drag coefficient
diffusion terms
Kronecker delta function
Young's modulus
eccentricity of the elliptical yield curve
strain rate tensor
nonlinear shear viscosity
Coriolis parameter
forces due to internal ice stress
Heaviside function
ice thickness
ice cutoff thickness
snow thickness on ice-covered fraction
ice mass,
ice strength
ice strength parameters
thermodynamic terms
internal ice stress tensor
air stress
water stress
the () components of ice velocity
10 meter air and surface water velocities
nonlinear bulk viscosity
surface elevation of the underlying water

Note that Hibler's variable is equivalent to our combination - his is the average thickness over the whole gridbox while our is the average thickness over the ice-covered fraction of the gridbox.


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:

 Ice diag1.png

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

The term 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 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 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:

where and 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, which can be up to 0.1 m, after which melting ice contributes to . 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
shortwave albedo of water
shortwave albedo of wet ice
shortwave albedo of dry ice
shortwave albedo of wet snow
shortwave albedo of dry snow
snow correction factor
= 2093 J kg K specific heat of ice
= 3990 J kg K specific heat of water
longwave emissivity of water
longwave emissivity of ice
longwave emissivity of snow
enthalpy of the ice/brine system
heat flux from the ocean into the ice
sensible heat
fraction of the solar heating transmitted through a lead into the water below
= 2.04 W m K thermal conductivity of ice
= 0.31 W m K thermal conductivity of snow
= 302 MJ m latent heat of fusion of ice
= 110 MJ m latent heat of fusion of snow
latent heat
incoming longwave radiation
& C/PSU coefficient in linear equation
contribution to equation from freezing water
heat flux out of the snow/ice surface
heat flux out of the ocean surface
heat flux up out of the ice
heat flux up into the ice
heat flux up through the snow
brine fraction in ice
= 910 density of ice
= 3.2 PSU salinity of the ice
incoming shortwave radiation
= W m K Stefan-Boltzmann constant
temperature of the bottom of the ice
temperature of the interior of the ice
temperature at the upper surface of the ice
temperature at the upper surface of the snow
freezing temperature
= melting temperature of ice
= 0 C melting temperature of snow
melt rate on the upper ice/snow surface
freeze rate at the air/water interface
rate of frazil ice growth
freeze rate at the ice/water interface
rate of run-off of surface melt water
snowfall rate
snow melt rate

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

Ice diag2.png

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 .

The surface flux to the air is:

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 , as is the heat flux through the snow . Setting , we can solve for by setting and linearizing in . As in Parkinson and Washington, if is found to be above the melting temperature, it is set to and the extra energy goes into melting the snow or ice:

Note that 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 . 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 and that the freezing temperature is a linear function of salinity. The brine fraction is given by

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

Substituting in for and differentiating gives:

Inside the snow, we have

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

These can be set equal to each other to solve for


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

Note that in the absence of snow, becomes zero and we recover the formula for the no-snow case in which .

At the bottom of the ice, we have

The difference between and goes into the enthalpy of the ice:

We can use the chain rule to obtain an equation for timestepping :


Ocean surface boundary conditions

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

where the relevant variables are in the following table:

Variable Value Definition
3.14 factor
0.4 von Karman's constant
vertical viscosity of seawater
kinematic viscosity of seawater
13.0 molecular Prandtl number
0.85 turbulent Prandtl number
surface salinity
internal ocean salinity
2432 molecular Schmidt number
stress on the ocean from the ice
stress on the ocean from the wind
surface ocean temperature
internal ocean temperature
friction velocity
roughness parameter

surface salinity () in the presense of ice. We also have and at the uppermost computed ocean point below the surface. In order to solve for and , we assume a Yaglom and Kader, 1974 logarithmic boundary layer. The upper ocean heat flux is:


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


The ocean model receives the following heat and salt fluxes:

Mellor and Kantha describe solving simultaneously for the five unknowns , , , and . Instead, we use the old value of to find and therefore . Using the new value of , solve for a new value of and then find the new as the freezing temperature for that salinity:

Note that the term 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:

The variables are defined in this table:

Variable Value Definition
1994 J kg K specific heat of ice
3987 J kg K specific heat of water
fraction of water that froze
3.16e5 J kg latent heat of fusion
mass of ice formed
mass of water before freezing
mass of water after freezing
constant in freezing equation
constant in freezing equation
salinity before freezing
salinity after freezing
temperature before freezing
temperature after freezing

Defining and dropping terms of order leads to:

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

We can then solve for :

The ocean is checked at each depth 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: