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

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

gravity | |

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.

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

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:

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

where

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

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 :

where

## 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 | |

evaporation | ||

0.4 | von Karman's constant | |

vertical viscosity of seawater | ||

kinematic viscosity of seawater | ||

precipitation | ||

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:

where

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

where

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: