SEAICE.850777 2001/10/15 This document contains an overview of the thermodynamic and dynamic sea ice models used in GISS climate models, and also their programming details. --------------------------- Thermodynamic Sea Ice Model --------------------------- In the dynamic ocean models used at GISS, each horizontal grid box is either all ocean or all continent. Each horizontal ocean grid cell can be fractionally covered by sea ice which has uniform thickness within the grid cell. In each ocean grid cell with sea ice, one mass variable indicates the snow amount and a second mass variable indicates the ice thickness. There are four separate temperature layers; the first thermal layer is generally thin but its thickness depends weakly on the snow amount. Table 1 lists the thermodynamic sea ice variables for each ocean cell and Table 2 lists the different values of sea ice parameters for the different models used at GISS. The specified ocean models do not calculate sea ice cover, but rather use climatological or time-varying observed ocean surface temperature data sets including horizontal sea ice cover. Table 1: Prognostic thermodynamic sea ice variables --------------------------------------------------- RSI = horizontal ratio of sea ice cover to ocaen cover MSI1 (kg/m) = fixed upper ice thickness plus variable snow amount MSI2 (kg/m) = variable lower ice thickness HSIn (J/m) = heat content (including negative latent heat) Table 2: Thermodynamic sea ice parameters ----------------------------------------- Spec C Spec Q C Ocean model -ify Grid -ify Flux Grid GFDL MICOM ----------- ---- ---- ---- ---- ---- ---- ---- Atmosphere model C C B B B B B Max snow mass (kg/m) 92 92 100 100 92 Min temp layer 1 (kg/m) 46 46 46 46 91 Min ice thick (kg/m) 364 364 459 459 364 Rain compression factor 2 2 0 0 Density of ice (kg/m) 910 910 917 917 910 Density of snow (kg/m) 100 100 300 300 Therm conduc ice (W/Cm) 2.18 2.18 Therm conduc sno (W/Cm) .35 .35 FLEAD = lead parameter .06 .10 MELTSI frequency (1/hr) 1 The thermodynamic sea ice coding in the C grid models conserve mass and static energy precisely, by which is meant that any gain of sea ice mass or static energy causes an opposite compensating gain of mass or energy to either the atmosphere or liquid ocean. Except for the sea ice melting subroutine, MELTSI, the processes discussed below are applied each hour in the C-grid and Q-flux ocean models. The atmospheric model condensation subroutines provide the mass of precipitation, the energy of precipitation, and the age of snow to the sea ice subroutines. The age of snow, which influences the albedo, depends on both time and temperature. Snowfall on sea ice increases the snow mass. If the snow mass exceeds the maximum snow mass then some snow is compacted into ice bringing the snow mass down to .9 times the maximum snow mass. The rain compression factor determines how much snow is compressed into ice when rain falls on sea ice. The rain temperature equilibrates with the first thermal layer. Some or all rain may freeze or some snow or ice may melt. Unfrozen rain and melted ice runs into the ocean at 0C. This coding occurs in subroutine PRECSI. The radiation subroutines provide surface insolation and downward thermal radiation flux to the sea ice subroutines. The spectrally integrated sea ice albedo varies from .45 for ice with out snow to .95 for ice with deep fresh snow. Insolation for each grid cell is distributed into the different surface types depending on the integrated surface albedo of each type. Insolation does not penetrate, but is absorbed in the first thermal layer. Atmospheric surface fluxes of water vapor and energy are calculated for each surface type for each grid cell in subroutine SURFCE. Because the mass of sea ice thermal layer 1 may be thin, SURFCE uses an implicit time scheme to calculate the surface fluxes. SURFCE also calculates the thermal flux between layers 1 and 2 because it is needed for the implicit time scheme. The thermal conductivity between layers 1 and 2 depends on the amount of snow. Subroutine SEAICE receives the atmospheric surface fluxes and the flux between thermal layers 1 and 2 from SURFCE. SEAICE calculates the thermal fluxes between thermal layers 2 and 3, 3 and 4, and between layer 4 and the ocean below. SEAICE then applies these fluxes to the layers. If the sea ice temperature in layers 1 or 4 tends to rise above 0C, then the ice tempreature stays at 0C and snow or sea ice is melted which runs into the ocean below. When surface fluxes tend to cause the first layer of the open ocean to cool below the freezing point, the ocean stays at the freezing point and the required amount of sea ice is formed with minimum thickness. If sea ice tends to draw sufficient heat from the ocean below to cool it to below freezing, then the ocean stays at the freezing point and sea ice thickens. If the horizontal open ocean fraction becomes less than the minimum lead criteria, then sea ice is contracted horizontally, thus increasing its thickness. Sea ice contains no salt, so when sea ice forms or thickens, the salt concentration of the first ocean layer increases. This coding occurs in subroutine OSOURC. If the temperature of the first ocean layer rises above 0C, then sea ice is melted vertically and horizontally drawing the necessary energy from the ocean. If sea ice becomes thinner than the minimum sea ice thickness, then it is contracted horizontally so that it remains at the minimum sea ice thickness. This coding occurs in subroutine MELTSI. Many of the above processes require relayering of the mass and thermal layers so that the requirements indicated in Table 2 are continually obeyed. Relayering uses the (diffusive) upstream scheme for thermal advection. Relayering is kept to a minimum in each subroutine by calculating at each layer edge the net mass flux that transports the thermal flux as opposed to transporting the thermal fluxes by different processes whose mass fluxes may be of opposite sign. There are two reasons to break the sea ice coding into several subroutines. One reason is to modularize the code, but the more important reason is to calculate model diagnostics. Model diagnostics break the sea ice changes into individual terms and confirm the model's precise mass and energy conservation. Details of subroutine SURFCE (calculates the surface fluxes) ------------------------------------------------------------ The thermal conductivity between layers 1 and 2 divided by the distance between the layer centers is given by: dF1/dT (W/Cm) = 2. / [(MSI1-SNOW)/ii + SNOW/ss], where SNOW (kg/m) = mass of snow, i,s (kg/m) = densities of ice and snow, and i,s (W/Cm) = thermal conductivities of ice and snow. Details of subroutine SEAICE (applies surface fluxes to sea ice) ---------------------------------------------------------------- Dew falling on sea ice increases the ice amount but does not change the snow amount. Snow or ice that is melted from thermal layer 1, falls directly into the ocean without having a chance to refreeze in the lower layers, contrary to the coding for land ice. Details of subroutine OSOURC (forms sea ice by cold ocean) ---------------------------------------------------------- At the beginning of OSOURC, the first layer of the ocean is horizontally partitioned into an open ocean fraction and a sea ice covered fraction. The open ocean fraction receives the open ocean surface fluxes which may cause sea ice to form on the open ocean with the minimum sea ice thickness. The sea ice covered fraction receives the sea ice - ocean flux which may cause additional sea ice to form that increases the ice thickness. At the end of OSOURC, the new and old sea ice are merged into a single thickness conserving sea ice mass, sea ice energy, and the updated horizontal sea ice fraction. At the end of OSOURC, the horizontal fraction of open ocean, 1-RSI, in a grid cell must be at least FLEADi/(MSI1+MSI2). If it is less, then the sea ice is contracted horizontally to maintain this constraint. Details of subroutine MELTSI (melts sea ice by warm water) ---------------------------------------------------------- In the C grid ocean models, the ocean energy above 0C in the first ocean layer (12 m) is calculated. Each hour, only 1/72 of this energy is used to melt sea ice. The fraction of this energy used to reduce the vertical thickness of sea ice is RSI; the fraction to melt sea ice horizontally is 1-RSI. Possible improvements to sea ice thermodynamics ----------------------------------------------- Add prognostic salt mass for each thermal layer. Would need to specify amount of salt when sea ice first forms and rate at which salt penetrates downward. Should rates be temperature dependent? Possible effects may be to reduce ocean overturning and to determine sea ice temperature and melting temperature. Note that heat content is the prognostic variable. The sea ice temperature and melting temperature could be made dependent on the salinity via the specific heat capacity and the latent heat of melting. The specific heat capacity and the latent heat of melting at 0C are currently constant. Make the temperature of both sea ice and water at the interface be the freezing temperature of ocean water. The melting or freezing at the interface will depend on how much heat each medium can deliver to the interface. Compact snow into ice that would be under water due to the weight of snow avove it. Allow solar radiation to penetrate through the ice. Improve sea ice albedo. --------------------- Dynamic Sea Ice Model --------------------- The dynamic sea ice model (used presently only in the C grid ocean models) adds five prognostic variables that are listed in Table 3. The time rate of change of sea ice velocity components are given by: u/t = T/M - (u-U)CV/M - uBRM + v(f+utan/a) - P/ - F/ v/t = T/M - (v-V)CV/M - vBRM - u(f+utan/a) - P/ - F/ where the variables are defined in Table 4. Table 3: Prognostic dynamic sea ice variables --------------------------------------------- RSIX = east-west gradient of horizontal fraction of sea ice RSIY = north-south gradient of horizontal fraction of sea ice F (Pa) = internal sea ice pressure u,v (m/s) = sea ice velocity components Table 4: Variables used in dynamical equations for sea ice ---------------------------------------------------------- a (m) = radius of the Earth = 6375000 B (m/kgs) = island and coast line blocking factor C = sea ice - ocean drag coefficient f (1/s) = Coriolis parameter F,F (Pa/m) = gradient of internal sea ice pressure M (kg/m) = sea ice mass per unit area = MSI1 + MSI2 P,P (Pa/m) = pressure gradient from ocean and sea ice tilt and atmospheric pressure R = horizontal sea ice fraction = RSI u,v (m/s) = sea ice velocity components u,v (m/s) = sea ice velocity components at beginning of step U,V (m/s) = ocean velocity components of ocean layer 1 t (s) = time (an integration time step is one hour) T,T (N/m) = atmospheric momentum stress components V (m/s) = [(u-U) + (v-V)] = latitude (kg/m) = density of ice = 910 The equations for u and v are simultaneous linear differential equations: u/t = Au + Bv + C (C here is not the drag coefficient) v/t = Av - Bu + D. They can be solved in closed form: u = [u + (AC-BD)/(A+B)] exp(At) cos(Bt) + [v + (AD+BC)/(A+B)] exp(At) sin(Bt) - (AC-BD)/(A+B) v = [v + (AD+BC)/(A+B)] exp(At) cos(Bt) + [u + (AC-BD)/(A+B)] exp(At) sin(Bt) - (AD+BC)/(A+B) In the finite difference representation, the sea ice velocity components are defined on the C grid (i.e. they are perpendicular to the grid cell edges). This is ideal for sea ice advection and for calculating the pressure gradient forces, but the Coriolis term is calculated less accurately. The pair of equations are solved twice, once for u using averaged values of v, and again for v using averaged values of u. The two linear equations for sea ice velocity components are integrated in closed form. Given the sea ice velocity at the beginning of a time step, the time integral of sea ice velocity during the time step and the sea ice velocity at the end of the time step are calculated. These closed form integrations are iterated upon to allow the internal sea ice pressure to equilibrate with the current parameter values. A modified Linear Upstream Scheme is used to advect sea ice. This should produce smooth distributions of sea ice yet be considerably less diffusive than a standard upstream scheme. The modification is to insist that RSIRSIX and RSIRSIY be between 0 and 1 at each grid cell edge. This causes additional diffusion which is not present in the normal Linear Upstream Scheme. This modification could be eliminated. Implementation of each term on the right hand side of the sea ice velocity tendency equations follows. Note that the advective terms for sea ice velocity are ignored. Subroutine SURFCE determines the surface momentum stress between the atmosphere and sea ice. Inside SURFCE, the sea ice velocity is ignored, i.e. it is assumed to be zero. The momentum stress for each component is assumed to be constant during integration of sea ice velocity. Momentum stress is calculated on the primary grid and is weighted to the C-grid velocity components proportional to sea ice area of the two grid cells. For sea ice - ocean drag, the magnitude of sea ice minus ocean velocity is calculated once using the sea ice velocity at the beginning of the time step. This can lead to hourly oscillations. The C-grid model's minimum value of the drag coefficient, C, is .0015 for minimum .4 m thick ice. The coefficient increases linearly with ice thickness and is equal to .005 for 3 m thick ice. The time integrated drag also accelerates the ocean below the ice. The island and coast line blocking factor, B, for sea ice advection is used to decelerate the sea ice velocity. It is used sparingly. It helps prevent ice from leaving the Arctic Ocean via the Northwest Passage. The blocking effect increases proportional to RSI times the sea ice mass per unit area in the ocean grid cell. The Coriolis and metric terms are simple using a four point average of the other component. The form of this term and the lack of momentum advection prevent conservation of angular momentum. The components of the sea ice pressure force, P and P, are caused by the ocean and sea ice tilt and the atmospheric pressure. They are assumed to be constant during the integration step for sea ice velocity. The terms for internal sea ice pressure gradient, F and F, follow Flato and Hibler [1992] with a few modifications. The prognostic internal sea ice pressure is updated using an iteration scheme with at most 20 iterations. In each iteration, the internal sea ice pressure gradient term is used in solving the two simultaneous linear equations for sea ice velocity. The time integral of this sea ice velocity estimates the sea ice area crossing each grid cell edge. If the sea ice area is convergent in a cell and the sea ice pressure is below its maximum, then the ice pressure increment is calculated so that the ice area convergence would be zero were this increment the only change to the sea ice pressure field. If the sea ice area is divergent in a cell and the sea ice pressure is positive, then ice pressure increment is calculated so that the ice area divergence would be zero were this increment the only change to the sea ice pressure field. The increments are applied to all grid cells simultaneously before the next iteration. The formula for maximum sea ice pressure is: Fmax (Pa) = 27500 M exp[20(RSI-1)] / . The scheme used here for internal sea ice pressure causes sea ice area convergence to approach zero, and applies all pressure increments simultaneously at the end of an iteration. Grid cell order of calculations in unimportant. The largest increment for the 20-th iteration is usually less than 100 Pa. This scheme differs from that of Flato and Hibler [1992] which causes sea ice velocity convergence to approach zero, and applies the pressure increments as they are calculated for each grid cell. Their method converges more rapidly, but if convergence is not complete, then the order of calculations affects the results. Summary ------- The advection of sea ice by the modified Linear Upstream Scheme is certainly superior to the simple upstream scheme or second order differencing scheme used by other models. Sea ice dynamics and advection are performed every hour with prognostic variables for u, v and F. u and v are integrated with one hour time steps. Flato and Hibler [1992] often use daily time steps, build up F from scratch each step, and u and v are determined by setting the time derivatives to zero. In iteratively solving for F, we want sea ice area to converge to zero, whereas Flato and Hibler want sea ice velocity to be non-divergent. With our wide Northwest Passage, we will definitely need an island or coast line blocking term to prevent too much sea ice from escaping from the Arctic Ocean. Possible improvements to sea ice dynamics ----------------------------------------- Increase the sea ice drag parallel to a coast line. Let subroutine SURFCE recognize the sea ice and ocean surface velocities instead of letting it assume that the surface has 0 velocity everywhere. Modify internal sea ice pressure to cause ֎(u,v) to converge to 0. Let parallel shear stresses of sea ice affect u and v. Let the ocean surface velocity, that the sea ice sees, be calculated from the ocean first layer velocity using an Ekman spiral. Currently, the sea ice sees the ocean first layer mean velocity. ------------------------------ Land ice calving in Antarctica ------------------------------ The purposes of putting land ice calving into the atmosphere-ocean model are these: 1. In the prior model, sea level dropped about 6.5 mm/year. Water evaporated from the ocean and precipitated onto the ice sheets and ice caps. Only a fraction of it was returned to the atmosphere via evaporation or to the ocean via melting. 2. When sea ice advection was first implemented, the southern hemisphere sea ice cover was too small. To correct these deficiencies, land ice calving was implemented from Antarctica. For each grid cell in Antarctica, the net annual snow accumulation (precipitation minus evaporation minus liquid runoff) is calculated from a 5 year run of the atmospheric model. Assuming that Antarctica is in mass equilibrium, the net snow accumulation is converted into land ice transport. Using the steepest ice top topography gradients, land ice transport is directed downhill toward the ocean, usually. The total land ice transport which reaches the coast is 298610 kg/year (for run A073b), compared with the IPCC [1996] estimate of 201610 kg/year. 40% of the 298610 kg/year of land ice transport reaches a known ice shelf. There the transports of ice mass and ice energy are added to accumulating arrays. When the ice mass accumulating array of a grid cell exceeds 10 kg, the ice mass and energy are released as an iceberg into the ocean. An iceberg's movement in the ocean is a linear combination of the surface air velocity, the ocean velocity, and the sea ice velocity. An iceberg draws heat from the first five ocean layers (158 m) to warm up and melt. The rate of heating depends on the ocean temperature and the square root of the of the iceberg mass. The 60% of land ice transport that does not reach an ice shelf is continually dumped into the ocean as 2 m thick sea ice. Figures ------- Figure 1. February and August northern hemisphere ocean ice cover (10 m) for the observations [NASA/GSFC] and two atmosphere-ocean model control simulations with constant 1950 atmospheric composition. Figure 2. Same as Figure 1 but for southern hemisphere. Figure 3. Global vector plot of sea ice mass flux from a 10-year annual average of a control simulation. Sea ice mass flux is proportional to the area of the arrows. Figure 4. Observed sea ice velocity on a northern polar projection. Figure 5. A control simulation's annual sea ice velocity on a northern polar projection. Figure 6. Annual glacial ice flow from atmospheric model. Ice shelves near Antarctica are indicated by the letter S.