2. Science Guide¶
2.1. Coupling with other climate model components¶
The sea ice model exchanges information with the other model components via a flux coupler. CICE has been coupled into numerous climate models with a variety of coupling techniques. This document is oriented primarily toward the CESM Flux Coupler [41] from NCAR, the first major climate model to incorporate CICE. The flux coupler was originally intended to gather state variables from the component models, compute fluxes at the model interfaces, and return these fluxes to the component models for use in the next integration period, maintaining conservation of momentum, heat, and fresh water. However, several of these fluxes are now computed in the ice model itself and provided to the flux coupler for distribution to the other components, for two reasons. First, some of the fluxes depend strongly on the state of the ice, and vice versa, implying that an implicit, simultaneous determination of the ice state and the surface fluxes is necessary for consistency and stability. Second, given the various ice types in a single grid cell, it is more efficient for the ice model to determine the net ice characteristics of the grid cell and provide the resulting fluxes, rather than passing several values of the state variables for each cell. These considerations are explained in more detail below.
The fluxes and state variables passed between the sea ice model and the
CESM flux coupler are listed in Table 1. By convention,
directional fluxes are positive downward. In CESM, the sea ice model may
exchange coupling fluxes using a different grid than the computational
grid. This functionality is activated using the namelist variable
gridcpl_file. Another namelist variable highfreq, allows the
high-frequency coupling procedure implemented in the Regional Arctic
System Model (RASM). In particular, the relative atmosphere-ice velocity
(\(\vec{U}_a-\vec{u}\)) is used instead of the full atmospheric
velocity for computing turbulent fluxes in the atmospheric boundary
layer.
Table 1: Data exchanged between the CESM flux coupler and the sea ice model
| Variable | Description | Interaction with flux coupler |
|---|---|---|
| \(z_o\) | Atmosphere level height | From atmosphere model via flux coupler to sea ice model |
| \(\vec{U}_a\) | Wind velocity | From atmosphere model via flux coupler to sea ice model |
| \(Q_a\) | Specific humidity | From atmosphere model via flux coupler to sea ice model |
| \(\rho_a\) | Air density | From atmosphere model via flux coupler to sea ice model |
| \(\Theta_a\) | Air potential temperature | From atmosphere model via flux coupler to sea ice model |
| \(T_a\) | Air temperature | From atmosphere model via flux coupler to sea ice model |
| \(F_{sw\downarrow}\) | Incoming shortwave radiation (4 bands) | From atmosphere model via flux coupler to sea ice model |
| \(F_{L\downarrow}\) | Incoming longwave radiation | From atmosphere model via flux coupler to sea ice model |
| \(F_{rain}\) | Rainfall rate | From atmosphere model via flux coupler to sea ice model |
| \(F_{snow}\) | Snowfall rate | From atmosphere model via flux coupler to sea ice model |
| \(F_{frzmlt}\) | Freezing/melting potential | From ocean model via flux coupler to sea ice model |
| \(T_w\) | Sea surface temperature | From ocean model via flux coupler to sea ice model |
| \(S\) | Sea surface salinity | From ocean model via flux coupler to sea ice model |
| \(\nabla H_o\) | Sea surface slope | From ocean model via flux coupler to sea ice model |
| \(\vec{U}_w\) | Surface ocean currents | From ocean model via flux coupler to sea ice model |
| \(\vec{\tau}_a\) | Wind stress | From sea ice model via flux coupler to atmosphere model |
| \(F_s\) | Sensible heat flux | From sea ice model via flux coupler to atmosphere model |
| \(F_l\) | Latent heat flux | From sea ice model via flux coupler to atmosphere model |
| \(F_{L\uparrow}\) | Outgoing longwave radiation | From sea ice model via flux coupler to atmosphere model |
| \(F_{evap}\) | Evaporated water | From sea ice model via flux coupler to atmosphere model |
| \(\alpha\) | Surface albedo (4 bands) | From sea ice model via flux coupler to atmosphere model |
| \(T_{sfc}\) | Surface temperature | From sea ice model via flux coupler to atmosphere model |
| \(F_{sw\Downarrow}\) | Penetrating shortwave radiation | From sea ice model via flux coupler to ocean model |
| \(F_{water}\) | Fresh water flux | From sea ice model via flux coupler to ocean model |
| \(F_{hocn}\) | Net heat flux to ocean | From sea ice model via flux coupler to ocean model |
| \(F_{salt}\) | Salt flux | From sea ice model via flux coupler to ocean model |
| \(\vec{\tau}_w\) | Ice-ocean stress | From sea ice model via flux coupler to ocean model |
| \(F_{bio}\) | Biogeochemical fluxes | From sea ice model via flux coupler to ocean model |
| \(a_{i}\) | Ice fraction | From sea ice model via flux coupler to both ocean and atmosphere models |
| \(T^{ref}_{a}\) | 2m reference temperature (diagnostic) | From sea ice model via flux coupler to both ocean and atmosphere models |
| \(Q^{ref}_{a}\) | 2m reference humidity (diagnostic) | From sea ice model via flux coupler to both ocean and atmosphere models |
| \(F_{swabs}\) | Absorbed shortwave (diagnostic) | From sea ice model via flux coupler to both ocean and atmosphere models |
The ice fraction \(a_i\) (aice) is the total fractional ice coverage of a grid cell. That is, in each cell,
where \(a_{i}\) is the sum of fractional ice areas for each category of ice. The ice fraction is used by the flux coupler to merge fluxes from the ice model with fluxes from the other components. For example, the penetrating shortwave radiation flux, weighted by \(a_i\), is combined with the net shortwave radiation flux through ice-free leads, weighted by (\(1-a_i\)), to obtain the net shortwave flux into the ocean over the entire grid cell. The flux coupler requires the fluxes to be divided by the total ice area so that the ice and land models are treated identically (land also may occupy less than 100% of an atmospheric grid cell). These fluxes are “per unit ice area” rather than “per unit grid cell area.”
In some coupled climate models (for example, recent versions of the U.K.
Hadley Centre model) the surface air temperature and fluxes are computed
within the atmosphere model and are passed to CICE. In this case the
logical parameter calc_Tsfc in ice_therm_vertical is set to false.
The fields fsurfn (the net surface heat flux from the atmosphere), flatn
(the surface latent heat flux), and fcondtopn (the conductive flux at
the top surface) for each ice thickness category are copied or derived
from the input coupler fluxes and are passed to the thermodynamic driver
subroutine, thermo_vertical. At the end of the time step, the surface
temperature and effective conductivity (i.e., thermal conductivity
divided by thickness) of the top ice/snow layer in each category are
returned to the atmosphere model via the coupler. Since the ice surface
temperature is treated explicitly, the effective conductivity may need
to be limited to ensure stability. As a result, accuracy may be
significantly reduced, especially for thin ice or snow layers. A more
stable and accurate procedure would be to compute the temperature
profiles for both the atmosphere and ice, together with the surface
fluxes, in a single implicit calculation. This was judged impractical,
however, given that the atmosphere and sea ice models generally exist on
different grids and/or processor sets.
2.1.1. Atmosphere¶
The wind velocity, specific humidity, air density and potential temperature at the given level height \(z_\circ\) are used to compute transfer coefficients used in formulas for the surface wind stress and turbulent heat fluxes \(\vec\tau_a\), \(F_s\), and \(F_l\), as described below. Wind stress is arguably the primary forcing mechanism for the ice motion, although the ice–ocean stress, Coriolis force, and slope of the ocean surface are also important [69]. The sensible and latent heat fluxes, \(F_s\) and \(F_l\), along with shortwave and longwave radiation, \(F_{sw\downarrow}\), \(F_{L\downarrow}\) and \(F_{L\uparrow}\), are included in the flux balance that determines the ice or snow surface temperature when calc_Tsfc = true. As described in Section Thermodynamics, these fluxes depend nonlinearly on the ice surface temperature \(T_{sfc}\). The balance equation is iterated until convergence, and the resulting fluxes and \(T_{sfc}\) are then passed to the flux coupler.
The snowfall precipitation rate (provided as liquid water equivalent and converted by the ice model to snow depth) also contributes to the heat and water mass budgets of the ice layer. Melt ponds generally form on the ice surface in the Arctic and refreeze later in the fall, reducing the total amount of fresh water that reaches the ocean and altering the heat budget of the ice; this version includes two new melt pond parameterizations. Rain and all melted snow end up in the ocean.
Wind stress and transfer coefficients for the turbulent heat fluxes are computed in subroutine atmo_boundary_layer following [41]. For clarity, the equations are reproduced here in the present notation.
The wind stress and turbulent heat flux calculation accounts for both stable and unstable atmosphere–ice boundary layers. Define the “stability”
where \(\kappa\) is the von Karman constant, \(g\) is gravitational acceleration, and \(u^*\), \(\Theta^*\) and \(Q^*\) are turbulent scales for velocity, temperature, and humidity, respectively:
The wind speed has a minimum value of 1 m/s. We have ignored ice motion in \(u^*\), and \(T_{sfc}\) and \(Q_{sfc}\) are the surface temperature and specific humidity, respectively. The latter is calculated by assuming a saturated surface, as described in Section Thermodynamic surface forcing balance.
Neglecting form drag,the exchange coefficients \(c_u\), \(c_\theta\) and \(c_q\) are initialized as
and updated during a short iteration, as they depend upon the turbulent scales. The number of iterations is set by the namelist variable natmiter. (For the case with form drag, see section Variable exchange coefficients.) Here, \(z_{ref}\) is a reference height of 10m and \(z_{ice}\) is the roughness length scale for the given sea ice category. \(\Upsilon\) is constrained to have magnitude less than 10. Further, defining \(\chi = \left(1-16\Upsilon\right)^{0.25}\) and \(\chi \geq 1\), the “integrated flux profiles” for momentum and stability in the unstable (\(\Upsilon <0\)) case are given by
In a departure from the parameterization used in [41], we use profiles for the stable case following [40],
The coefficients are then updated as
where \(\lambda = \ln\left(z_\circ/z_{ref}\right)\). The first iteration ends with new turbulent scales from equations (2). After five iterations the latent and sensible heat flux coefficients are computed, along with the wind stress:
where \(L_{vap}\) and \(L_{ice}\) are latent heats of vaporization and fusion, \(\rho_a\) is the density of air and \(c_p\) is its specific heat. Again following [40], we have added a constant to the sensible heat flux coefficient in order to allow some heat to pass between the atmosphere and the ice surface in stable, calm conditions.
The atmospheric reference temperature \(T_a^{ref}\) is computed from \(T_a\) and \(T_{sfc}\) using the coefficients \(c_u\), \(c_\theta\) and \(c_q\). Although the sea ice model does not use this quantity, it is convenient for the ice model to perform this calculation. The atmospheric reference temperature is returned to the flux coupler as a climate diagnostic. The same is true for the reference humidity, \(Q_a^{ref}\).
Additional details about the latent and sensible heat fluxes and other quantities referred to here can be found in Section Thermodynamic surface forcing balance.
For CICE run in stand-alone mode (i.e., uncoupled), the AOMIP shortwave and longwave radiation formulas are available in ice_forcing.F90. In function longwave_rosati_miyakoda, downwelling longwave is computed as
where the atmospheric vapor pressure (mb) is \(e_a = 1000 Q_a/(0.622+0.378Q_a)\), \(\epsilon=0.97\) is the ocean emissivity, \(\sigma\) is the Stephan-Boltzman constant, \(f_{cld}\) is the cloud cover fraction, and \(T_a\) is the surface air temperature (K). The first term on the right is upwelling longwave due to the mean (merged) ice and ocean surface temperature, \(T_s\) (K), and the other terms on the right represent the net longwave radiation patterned after [60]. The downwelling longwave formula of [57] is also available in function longwave_parkinson_washington:
The value of \(F_{lw\uparrow}\) is different for each ice thickness category, while \(F_{lw\downarrow}\) depends on the mean value of surface temperature averaged over all of the thickness categories and open water.
The AOMIP shortwave forcing formula (in subroutine compute_shortwave) incorporates the cloud fraction and humidity through the atmospheric vapor pressure:
where \(\cos Z\) is the cosine of the solar zenith angle.
2.1.2. Ocean¶
New sea ice forms when the ocean temperature drops below its freezing temperature. In the Bitz and Lipscomb thermodynamics, [7] \(T_f=-\mu S\), where \(S\) is the seawater salinity and \(\mu=0.054 \ ^\circ\)/ppt is the ratio of the freezing temperature of brine to its salinity (linear liquidus approximation). For the mushy thermodynamics, \(T_f\) is given by a piecewise linear liquidus relation. The ocean model calculates the new ice formation; if the freezing/melting potential \(F_{frzmlt}\) is positive, its value represents a certain amount of frazil ice that has formed in one or more layers of the ocean and floated to the surface. (The ocean model assumes that the amount of new ice implied by the freezing potential actually forms.)
If \(F_{frzmlt}\) is negative, it is used to heat already existing ice from below. In particular, the sea surface temperature and salinity are used to compute an oceanic heat flux \(F_w\) (\(\left|F_w\right| \leq \left|F_{frzmlt}\right|\)) which is applied at the bottom of the ice. The portion of the melting potential actually used to melt ice is returned to the coupler in \(F_{hocn}\). The ocean model adjusts its own heat budget with this quantity, assuming that the rest of the flux remained in the ocean.
In addition to runoff from rain and melted snow, the fresh water flux \(F_{water}\) includes ice melt water from the top surface and water frozen (a negative flux) or melted at the bottom surface of the ice. This flux is computed as the net change of fresh water in the ice and snow volume over the coupling time step, excluding frazil ice formation and newly accumulated snow. Setting the namelist option update_ocn_f to true causes frazil ice to be included in the fresh water and salt fluxes.
There is a flux of salt into the ocean under melting conditions, and a (negative) flux when sea water is freezing. However, melting sea ice ultimately freshens the top ocean layer, since the ocean is much more saline than the ice. The ice model passes the net flux of salt \(F_{salt}\) to the flux coupler, based on the net change in salt for ice in all categories. In the present configuration, ice_ref_salinity is used for computing the salt flux, although the ice salinity used in the thermodynamic calculation has differing values in the ice layers.
A fraction of the incoming shortwave \(F_{sw\Downarrow}\) penetrates the snow and ice layers and passes into the ocean, as described in Section Thermodynamic surface forcing balance.
Many ice models compute the sea surface slope \(\nabla H_\circ\) from geostrophic ocean currents provided by an ocean model or other data source. In our case, the sea surface height \(H_\circ\) is a prognostic variable in POP—the flux coupler can provide the surface slope directly, rather than inferring it from the currents. (The option of computing it from the currents is provided in subroutine evp_prep.) The sea ice model uses the surface layer currents \(\vec{U}_w\) to determine the stress between the ocean and the ice, and subsequently the ice velocity \(\vec{u}\). This stress, relative to the ice,
is then passed to the flux coupler (relative to the ocean) for use by the ocean model. Here, \(\theta\) is the turning angle between geostrophic and surface currents, \(c_w\) is the ocean drag coefficient, \(\rho_w\) is the density of seawater, and \(\hat{k}\) is the vertical unit vector. The turning angle is necessary if the top ocean model layers are not able to resolve the Ekman spiral in the boundary layer. If the top layer is sufficiently thin compared to the typical depth of the Ekman spiral, then \(\theta=0\) is a good approximation. Here we assume that the top layer is thin enough.
For CICE run in stand-alone mode (i.e., uncoupled), a thermodynamic slab ocean mixed-layer parameterization is available in ice_ocean.F90. The turbulent fluxes are computed above the water surface using the same parameterizations as for sea ice, but with parameters appropriate for the ocean. The surface flux balance takes into account the turbulent fluxes, oceanic heat fluxes from below the mixed layer, and shortwave and longwave radiation, including that passing through the sea ice into the ocean. If the resulting sea surface temperature falls below the salinity-dependent freezing point, then new ice (frazil) forms. Otherwise, heat is made available for melting the ice.
2.1.3. Variable exchange coefficients¶
In the default CICE setup, atmospheric and oceanic neutral drag coefficients (\(c_u\) and \(c_w\)) are assumed constant in time and space. These constants are chosen to reflect friction associated with an effective sea ice surface roughness at the ice–atmosphere and ice–ocean interfaces. Sea ice (in both Arctic and Antarctic) contains pressure ridges as well as floe and melt pond edges that act as discrete obstructions to the flow of air or water past the ice, and are a source of form drag. Following [76] and based on recent theoretical developments [49][48], the neutral drag coefficients can now be estimated from properties of the ice cover such as ice concentration, vertical extent and area of the ridges, freeboard and floe draft, and size of floes and melt ponds. The new parameterization allows the drag coefficients to be coupled to the sea ice state and therefore to evolve spatially and temporally. This parameterization is contained in the subroutine neutral_drag_coeffs and is accessed by setting formdrag = true in the namelist.
Following [76], consider the general case of fluid flow obstructed by N randomly oriented obstacles of height \(H\) and transverse length \(L_y\), distributed on a domain surface area \(S_T\). Under the assumption of a logarithmic fluid velocity profile, the general formulation of the form drag coefficient can be expressed as
where \(z_0\) is a roughness length parameter at the top or bottom surface of the ice, \(\gamma\) is a geometric factor, \(c\) is the resistance coefficient of a single obstacle, and \(S_c\) is a sheltering function that takes into account the shielding effect of the obstacle,
with \(D\) the distance between two obstacles and \(s_l\) an attenuation parameter.
As in the original drag formulation in CICE (sections Atmosphere and Ocean), \(c_u\) and \(c_w\) along with the transfer coefficients for sensible heat, \(c_{\theta}\), and latent heat, \(c_{q}\), are initialized to a situation corresponding to neutral atmosphere–ice and ocean–ice boundary layers. The corresponding neutral exchange coefficients are then replaced by coefficients that explicitly account for form drag, expressed in terms of various contributions as
The contributions to form drag from ridges (and keels underneath the ice), floe edges and melt pond edges can be expressed using the general formulation of equation (12) (see [76] for details). Individual terms in equation (15) are fully described in [76]. Following [4] the skin drag coefficient is parametrized as
where \(m_s\) (\(m_k\)) is a sheltering parameter that depends on the average sail (keel) height, \(H_s\) (\(H_k\)), but is often assumed constant, \(D_s\) (\(D_k\)) is the average distance between sails (keels), and \(c_{ss}\) (\(c_{sk}\)) is the unobstructed atmospheric (oceanic) skin drag that would be attained in the absence of sails (keels) and with complete ice coverage, \(a_{ice}=1\).
Calculation of equations (12) – (16) requires that small-scale geometrical properties of the ice cover be related to average grid cell quantities already computed in the sea ice model. These intermediate quantities are briefly presented here and described in more detail in [76]. The sail height is given by
and the distance between sails
where \(0<\alpha<1\) and \(0<\beta<1\) are weight functions, \(\alpha_{s}\) and \(\alpha_{k}\) are the sail and keel slope, \(\phi_s\) and \(\phi_k\) are constant porosities for the sails and keels, and we assume constant ratios for the average keel depth and sail height (\(H_k/H_s=R_h\)) and for the average distances between keels and between sails (\(D_k/D_s=R_d\)). With the assumption of hydrostatic equilibrium, the effective ice plus snow freeboard is \(H_{f}=\bar{h_i}(1-\rho_i/\rho_w)+\bar{h_s}(1-\rho_s/\rho_w)\), where \(\rho_i\), \(\rho_w\) and \(\rho_s\) are respectively the densities of sea ice, water and snow, \(\bar{h_i}\) is the mean ice thickness and \(\bar{h_s}\) is the mean snow thickness (means taken over the ice covered regions). For the melt pond edge elevation we assume that the melt pond surface is at the same level as the ocean surface surrounding the floes [20][21][22] and use the simplification \(H_p = H_f\). Finally to estimate the typical floe size \(L_A\), distance between floes, \(D_F\), and melt pond size, \(L_P\) we use the parameterizations of [49] to relate these quantities to the ice and pond concentrations. All of these intermediate quantities are available as history output, along with Cdn_atm, Cdn_ocn and the ratio Cdn_atm_ratio_n between the total atmospheric drag and the atmospheric neutral drag coefficient.
We assume that the total neutral drag coefficients are thickness category independent, but through their dependance on the diagnostic variables described above, they vary both spatially and temporally. The total drag coefficients and heat transfer coefficients will also depend on the type of stratification of the atmosphere and the ocean, and we use the parameterization described in section Atmosphere that accounts for both stable and unstable atmosphere–ice boundary layers. In contrast to the neutral drag coefficients the stability effect of the atmospheric boundary layer is calculated separately for each ice thickness category.
The transfer coefficient for oceanic heat flux to the bottom of the ice may be varied based on form drag considerations by setting the namelist variable fbot_xfer_type to Cdn_ocn; this is recommended when using the form drag parameterization. Its default value of the transfer coefficient is 0.006 (fbot_xfer_type = ’constant’).
2.2. Model components¶
The Arctic and Antarctic sea ice packs are mixtures of open water, thin first-year ice, thicker multiyear ice, and thick pressure ridges. The thermodynamic and dynamic properties of the ice pack depend on how much ice lies in each thickness range. Thus the basic problem in sea ice modeling is to describe the evolution of the ice thickness distribution (ITD) in time and space.
The fundamental equation solved by CICE is [74]:
where \({\bf u}\) is the horizontal ice velocity, \(\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})\), \(f\) is the rate of thermodynamic ice growth, \(\psi\) is a ridging redistribution function, and \(g\) is the ice thickness distribution function. We define \(g({\bf x},h,t)\,dh\) as the fractional area covered by ice in the thickness range \((h,h+dh)\) at a given time and location.
Equation (19) is solved by partitioning the ice pack in each grid cell into discrete thickness categories. The number of categories can be set by the user, with a default value \(N_C = 5\). (Five categories, plus open water, are generally sufficient to simulate the annual cycles of ice thickness, ice strength, and surface fluxes [6][45].) Each category \(n\) has lower thickness bound \(H_{n-1}\) and upper bound \(H_n\). The lower bound of the thinnest ice category, \(H_0\), is set to zero. The other boundaries are chosen with greater resolution for small \(h\), since the properties of the ice pack are especially sensitive to the amount of thin ice [50]. The continuous function \(g(h)\) is replaced by the discrete variable \(a_{in}\), defined as the fractional area covered by ice in the open water by \(a_{i0}\), giving \(\sum_{n=0}^{N_C} a_{in} = 1\) by definition.
Category boundaries are computed in init_itd using one of several formulas, summarized in Table 2. Setting the namelist variable kcatbound equal to 0 or 1 gives lower thickness boundaries for any number of thickness categories \(N_C\). Table 2 shows the boundary values for \(N_C\) = 5 and linear remapping of the ice thickness distribution. A third option specifies the boundaries based on the World Meteorological Organization classification; the full WMO thickness distribution is used if \(N_C\) = 7; if \(N_C\) = 5 or 6, some of the thinner categories are combined. The original formula (kcatbound = 0) is the default. Category boundaries differ from those shown in Table 2 for the delta-function ITD. Users may substitute their own preferred boundaries in init_itd.
Table 2 : Data exchanged between the CESM flux coupler and the sea ice model Lower boundary values for thickness categories, in meters, for the three distribution options ( kcatbound ) and linear remapping ( kitd = 1 ). In the WMO case, the distribution used depends on the number of categories used.
| distribution | original | round | WMO | ||
|---|---|---|---|---|---|
| kcatbound | 0 | 1 | 2 | ||
| \(N_C\) | 5 | 5 | 5 | 6 | 7 |
| categories | lower bound (m) | ||||
| 1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| 2 | 0.64 | 0.60 | 0.30 | 0.15 | 0.10 |
| 3 | 1.39 | 1.40 | 0.70 | 0.30 | 0.15 |
| 4 | 2.47 | 2.40 | 1.20 | 0.70 | 0.30 |
| 5 | 4.57 | 3.60 | 2.00 | 1.20 | 0.70 |
| 6 | 2.00 | 1.20 | |||
| 7 | 2.00 | ||||
In addition to the fractional ice area, \(a_{in}\), we define the following state variables for each category \(n\). In a change from previous CICE versions, we no longer carry snow and ice energy as separate variables; instead they and sea ice salinity are carried as tracers on snow and ice volume.
- \(v_{in}\), the ice volume, equal to the product of \(a_{in}\) and the ice thickness \(h_{in}\).
- \(v_{sn}\), the snow volume, equal to the product of \(a_{in}\) and the snow thickness \(h_{sn}\).
- \(e_{ink}\), the internal ice energy in layer \(k\), equal to the product of the ice layer volume, \(v_{in}/N_i\), and the ice layer enthalpy, \(q_{ink}\). Here \(N_i\) is the total number of ice layers, with a default value \(N_i = 4\), and \(q_{ink}\) is the negative of the energy needed to melt a unit volume of ice and raise its temperature to; it is discussed in Section Thermodynamics. (NOTE: In the current code, \(e_i<0\) and \(q_i<0\) with \(e_i = v_iq_i\).)
- \(e_{snk}\), the internal snow energy in layer \(k\), equal to the product of the snow layer volume, \(v_{sn}/N_s\), and the snow layer enthalpy, \(q_{snk}\), where \(N_s\) is the number of snow layers. (Similarly, \(e_s<0\) in the code.) CICE allows multiple snow layers, but the default value is \(N_s=1\).
- \(S_i\), the bulk sea ice salt content in layer \(k\), equal to the product of the ice layer volume and the sea ice salinity tracer.
- \(T_sfn\), the surface temperature.
Since the fractional area is unitless, the volume variables have units of meters (i.e., m\(^3\) of ice or snow per m\(^2\) of grid cell area), and the energy variables have units of J/m\(^2\).
The three terms on the right-hand side of Equation (19) describe three kinds of sea ice transport: (1) horizontal transport in \((x,y)\) space; (2) transport in thickness space \(h\) due to thermodynamic growth and melting; and (3) transport in thickness space \(h\) due to ridging and other mechanical processes. We solve the equation by operator splitting in three stages, with two of the three terms on the right set to zero in each stage. We compute horizontal transport using the incremental remapping scheme of [12] as adapted for sea ice by [46]; this scheme is discussed in Section Horizontal transport. Ice is transported in thickness space using the remapping scheme of [45], as described in Section Transport in thickness space. The mechanical redistribution scheme, based on [74], [61], [27], [19], and [47] is outlined in Section Mechanical redistribution. To solve the horizontal transport and ridging equations, we need the ice velocity \({\bf u}\), and to compute transport in thickness space, we must know the the ice growth rate \(f\) in each thickness category. We use the elastic-viscous-plastic (EVP) ice dynamics scheme of [33], as modified by [10], [31], [34] and [35], or a new elastic-anisotropic-plastic model [83][81][77] to find the velocity, as described in Section Dynamics. Finally, we use a thermodynamic model to compute \(f\) (Section Thermodynamics). The order in which these computations are performed in the code itself was chosen so that quantities sent to the coupler are consistent with each other and as up-to-date as possible. The Delta-Eddington radiative scheme computes albedo and shortwave components simultaneously, and in order to have the most up-to-date values available for the coupler at the end of the timestep, the order of radiation calculations is shifted. Albedo and shortwave components are computed after the ice state has been modified by both thermodynamics and dynamics, so that they are consistent with the ice area and thickness at the end of the step when sent to the coupler. However, they are computed using the downwelling shortwave from the beginning of the timestep. Rather than recompute the albedo and shortwave components at the beginning of the next timestep using new values of the downwelling shortwave forcing, the shortwave components computed at the end of the last timestep are scaled for the new forcing.
2.2.1. Tracers¶
The basic conservation equations for ice area fraction \(a_{in}\), ice volume \(v_{in}\), and snow volume \(v_{sn}\) for each thickness category \(n\) are
The ice and snow volumes can be written equivalently in terms of tracers, ice thickness \(h_{in}\) and snow depth \(h_{sn}\):
Although we maintain ice and snow volume instead of the thicknesses as state variables in CICE, the tracer form is used for volume transport (section Horizontal transport). There are many other tracers available, whose values are contained in the trcrn array. Their transport equations typically have one of the following three forms
Equation (25) describes the transport of surface temperature, whereas Equation (26) and Equation (27) describe the transport of ice and snow enthalpy, salt, and passive tracers such as volume-weighted ice age and snow age. Each tracer field is given an integer index, trcr_depend, which has the value 0, 1, or 2 depending on whether the appropriate conservation equation is Equation (25), Equation (26), or Equation (27), respectively. The total number of tracers is \(N_{tr}\ge 1\). In the default configuration there are two tracers: surface temperature and volume-weighted ice age. Tracers for melt ponds (Sections Tracers that depend on other tracers (e.g., melt ponds) and Melt ponds), level ice area and level ice volume (used to diagnose ridged ice area and volume, Section Mechanical redistribution) are also available. Users may add any number of additional tracers that are transported conservatively provided that trcr_depend is defined appropriately. See Section Tracers for guidance on adding tracers.
2.2.1.1. Tracers that depend on other tracers (e.g., melt ponds)¶
Tracers may be defined that depend on other tracers. Melt pond tracers provide an example (these equations pertain to cesm and topo tracers; level-ice tracers are similar with an extra factor of \(a_{lvl}\), see Equations (116)–(119). Conservation equations for pond area fraction \(a_{pnd}a_i\) and pond volume \(h_{pnd}a_{pnd}a_i\), given the ice velocity \(\bf u\), are
(These equations represent quantities within one thickness category; all melt pond calculations are performed for each category, separately.) Equation (29) expresses conservation of melt pond volume, but in this form highlights that the quantity tracked in the code is the pond depth tracer \(h_{pnd}\), which depends on the pond area tracer \(a_{pnd}\). Likewise, \(a_{pnd}\) is a tracer on ice area (Equation (28)), which is a state variable, not a tracer.
For a generic quantity \(q\) that represents a mean value over the ice fraction, \(q a_i\) is the average value over the grid cell. Thus for cesm or topo melt ponds, \(h_{pnd}\) can be considered the actual pond depth, \(h_{pnd}a_{pnd}\) is the mean pond depth over the sea ice, and \(h_{pnd}a_{pnd}a_i\) is the mean pond depth over the grid cell. These quantities are illustrated in Figure 1.
Figure 1
Figure 1 : Melt pond tracer definitions. The graphic on the right illustrates the grid cell fraction of ponds or level ice as defined by the tracers. The chart on the left provides corresponding ice thickness and pond depth averages over the grid cell, sea ice and pond area fractions.
Tracers may need to be modified for physical reasons outside of the “core” module or subroutine describing their evolution. For example, when new ice forms in open water, the new ice does not yet have ponds on it. Likewise when sea ice deforms, we assume that pond water (and ice) on the portion of ice that ridges is lost to the ocean.
When new ice is added to a grid cell, the grid cell total area of melt ponds is preserved within each category gaining ice, \(a_{pnd}^{t+\Delta t}a_{i}^{t+\Delta t} = a_{pnd}^{t}a_{i}^{t}\), or
Similar calculations are performed for all tracer types, for example tracer-on-tracer dependencies such as \(h_{pnd}\), when needed:
In this case (adding new ice), \(h_{pnd}\) does not change because \(a_{pnd}^{t+\Delta t}a_{i}^{t+\Delta t} = a_{pnd}^{t}a_{i}^{t}\).
When ice is transferred between two thickness categories, we conserve the total pond area summed over categories \(n\),
Thus,
This is more complicated because of the \(\Delta\) term on the right-hand side, which is handled manually in ice_itd.F90. Such tracer calculations are scattered throughout the code, wherever there are changes to the ice thickness distribution.
Note that if a quantity such as \(a_{pnd}\) becomes zero in a grid cell’s thickness category, then all tracers that depend on it also become zero. If a tracer should be conserved (e.g., aerosols and the liquid water in topo ponds), additional code must be added to track changes in the conserved quantity.
More information about the melt pond schemes is in section Melt ponds.
2.2.1.2. Ice age¶
The age of the ice, \(\tau_{age}\), is treated as an ice-volume tracer (trcr_depend = 1). It is initialized at 0 when ice forms as frazil, and the ice ages the length of the timestep during each timestep. Freezing directly onto the bottom of the ice does not affect the age, nor does melting. Mechanical redistribution processes and advection alter the age of ice in any given grid cell in a conservative manner following changes in ice area. The sea ice age tracer is validated in [32].
Another age-related tracer, the area covered by first-year ice \(a_{FY}\), is an area tracer (trcr_depend = 0) that corresponds more closely to satellite-derived ice age data for first-year ice than does \(\tau_{age}\). It is re-initialized each year on 15 September (yday = 259) in the northern hemisphere and 15 March (yday = 75) in the southern hemisphere, in non-leap years. This tracer is increased when new ice forms in open water, in subroutine add_new_ice in ice_therm_itd.F90. The first-year area tracer is discussed in [2].
2.2.1.3. Sea ice biogeochemistry¶
Ice algal photosynthesis leads to carbon fixation and pigment buildup throughout much of the pack ice in springtime, including warm layers in contact with seawater and the central ice matrix. Turbulence moves seed organisms and trace elements upward across the ocean interface from the mixed layer [3], and convection allows them to penetrate deep into the brine channel network. Gravity drainage is strongly coupled to carbonate and alkalinity relationships; bio-active gases including molecular oxygen, dimethyl sulfide (DMS) and methane exhibit dynamic cycling within and around the ice matrix. All may ultimately be transferred across the upper interface, into snow cover and toward the atmosphere. Chlorophyll is often concentrated in thin layers, capable of absorbing downwelling radiation and redistributing the implied heat locally.
Although biogeochemical activity occurs throughout the sea ice column, in the present release we restrict ourselves to simulations of activity taking place in the lowest few vertical centimeters. A thin intermediate zone located in this vicinity can be considered distinct because it is relatively warm, extremely porous and necessarily in close contact with the seawater nutrient supply. The ecosystems supported are therefore especially intense, so that they constitute a logical starting point for code development. The resulting band of biological activity will be referred to most often in what follows as the ‘bottom zone.’ It constitutes the dominant habitat for ice algae across most of the noncoastal Arctic and among land fast ice systems of the Southern Hemisphere. In the literature and some portions of the code, however, the alternate term ‘skeletal layer’ is used. For many purposes, the two are interchangeable. But in fact, the latter is a reference to dendritic structures more typically observable in the wintertime; ‘bottom zone’ is the more general term.
Ecological succession in the bottom zone interacts intimately with chemistry of the upper ocean–atmosphere system. Photosynthesis is constrained initially by light limitation but soon becomes a matter of resource availability [3][11]. Radiation intensity dictates the timing of bloom inception. The direct inflow of salt laden water can then be filtered immediately for the dissolved nutrient forms, which are converted locally into biomass. In association with the carbon storage, considerable DMS and molecular oxygen are generated [17]. Melt and flush periods cause bottom ice organisms to be rejected from the matrix mechanically. The dominant primary producers among them are pennate diatoms, a class subject to rapid sinking [42]. Hence the amount of bottom layer biological activity dictates high latitude nutrient and trace gas levels in the early spring, then carbon drawdown moving into the period of breakup.
Light transmission passes into the bottom zone through overlying layers of snow and ice. Attenuation by the ice algae is integrated and averaged over the bottom three centimeters in which they reside. Our current CICE release treats inflow of the primary Arctic nutrients nitrate and silicate, along with the major pathway for nitrogen recycling via ammonia. Both fertilizers and byproducts are permitted to exchange back out of the porous bottom volume on a time scale of order hours [3][42]. Biomass multiplies exponentially as light and nutrient restrictions are lifted, following the rising polar sun. Uptake accelerates until there is a transition to flux limitation [17]. At variable times and locations, a kinetic balance is often reached in which dissolved nitrogen is sequestered faster than it can leave the porous volume. Computational terms involved in the ecodynamic simulation resemble Michaelis-Menten enzyme kinetics, and have been adapted here from a series of pioneering ice biogeochemistry models [3][39][42]. CICE biogeochemical applications include a series of Pan-Arctic simulations, conceptually extending from primary production and carbon cycling to the release of trace gases into the ice domain [11][17].
The bottom ice biogeochemistry module described here is designed for ready attachment of extra-nutrient cycling, including byproduct and detrital processing. Our own mechanism for the generation of DMS is included in the release as an example. For full details please see descriptions and schematics in [17]. Generally speaking, the primary nutrient flow as nitrogen simultaneously supports formation of silicate frustules by the ice diatoms, along with carbon biomass production and the implied chlorophyll buildup. Synthesis of organosulfur metabolites is handled simultaneously and in direct proportion. DMS spins off of this subsystem at much higher levels than those familiar from open water studies, since the ice algae are especially abundant and suffer strong cryological, osmotic and oxidant stress in an extreme environmental regime [70]. The sulfur transformations are governed by a typical pattern of routing kinetics within and below the ice [17], which reduces the net yield of volatile gas but still permits considerable buildup and release during the spring thaw.
Biogeochemical processing is modelled as a single layer of reactive tracers attached to the sea ice bottom. Optional settings are available via the zbgc_nml namelist in ice_in. The prefix tr_bgc signifies a biogeochemical tracer, while the postscript _sk indicates that a particular quantity is restricted to the band of bottom or skeletal ice material. History fields are controlled in the icefields_bgc_nml namelist. As with other CICE history fields, the suffix _ai indicates that the field is multiplied by ice area and is therefore a grid-cell-mean quantity.
Setting the flag skl_bgc to true turns on a reduced ‘NP’ (Nitrogen-Plankton) biogeochemistry consisting of nitrate and one algal tracer. Ammonium (tr_bgc_Am_sk), silicate (tr_bgc_Sil_sk), dimethyl sulfoniopropionate (DMSP) in particulate form (tr_bgc_DMSPp_sk), DMSP in dissolved form (tr_bgc_DMSPd_sk), and volatile dimethyl sulfide (tr_bgc_DMS_sk) may also be included by setting the respective flags to true in ice_in.
All biogeochemical tracers (\(T_b\)) are brine concentrations multiplied by the skeletal layer thickness (\(h_{sk}\)). Bulk tracer concentrations are written to diagnostic files and are found by the expression \(T_b \phi_{sk}/h_{sk}\) where \(\phi_{sk}\) is the skeletal layer porosity. Both \(h_{sk}\) and \(\phi_{sk}\) are defined in ice_zbgc_shared.F90.
Tracers \(T_b\) are area conserved and follow the horizontal transport equation ([eq:transport_aT]). The initial concentrations of tracers are specified in subroutine init_bgc in ice_zbgc.F90. Silicate and nitrate may be initialized and forced with climatology. Parameters sil_data_type and nit_data_type are set in ice_in and take the values ‘clim’ (climatology) and ‘default’ (constant). nit_data_type may also take the value ‘sss’ which sets nitrate to the ocean salinity value. For climatology, the data directory bgc_data_dir must also be specified. If restore_bgc is true, then values are linearly restored to climatology according to the restoring timescale trestore.
For each horizontal grid point, local biogeochemical tracer equations are solved in ice_algae.F90. There are two types of ice–ocean tracer flux formulations: ‘Jin2006’ modelled after the growth rate dependent piston velocity of [39], and ‘constant’ modelled after the constant piston velocity of [17]. The formulation is specified in ice_in by setting bgc_flux_type equal to ‘Jin2006’ or ‘constant’.
In addition to horizontal advection and transport among thickness categories, biogeochemical tracers (\(T_b\) where \(b = 1,\ldots, N_b\)) satisfy a set of local coupled equations of the form
where \(R_b\) represents the nonlinear biochemical reaction terms detailed in [17] and \(\Delta z\) is a length scale representing the molecular sublayer of the ice–ocean interface. Its value is absorbed in the piston velocity parameters. The piston velocity \(w_b\) depends on the particular tracer and the flux formulation.
For ‘Jin2006’, the piston velocity is a function of ice growth and melt rates. All tracers (algae included) flux with the same piston velocity during ice growth, \(dh/dt > 0\):
with parameters \(m_1\), \(m_2\), \(m_3\) and \(p_g\) defined in skl_biogeochemistry. For ice melt, \(dh/dt < 0\), all tracers with the exception of ice algae flux with
with \(p_m\) defined in skl_biogeochemistry. The ‘Jin2006’ formulation also requires that for both expressions, \(|w_b| \leq 0.9 h_{sk}/\Delta t\). The concentration difference at the ice–ocean boundary for each tracer, \(\Delta T_b\), depends on the sign of \(w_b\). For growing ice, \(w_b <0\), \(\Delta T_b = T_b/h_{sk} - T_{io}\), where \(T_{io}\) is the ocean concentration of tracer \(i\). For melting ice, \(w_b > 0\), \(\Delta T_b = T_b/h_{sk}\).
In ‘Jin2006’, the algal tracer (\(N_a\)) responds to ice melt in the same manner as the other tracers from Equation (36). However, this is not the case for ice growth. Unlike dissolved nutrients, algae are able to cling to the ice matrix and resist expulsion during desalination. For this reason, algal tracers do not flux between ice and ocean during ice growth unless the ice algal brine concentration is less than the ocean algal concentration (\(N_o\)). Then the ocean seeds the sea ice concentration according to
The ‘constant’ formulation uses a fixed piston velocity (PVc) for positive ice growth rates for all tracers except \(N_a\). As in ‘Jin2006’, congelation ice growth seeds the sea ice algal population according to Equation (37) when \(N_a < N_o h_{sk}/\phi_{sk}\). For bottom ice melt, all tracers follow the prescription
All of the necessary fluxes required for coupling to an ocean model component are contained in a single array, flux_bio.
This bottom-specific biogeochemistry module necessarily simplifies the full sea ice ecosystem. For the moment, zooplankton comprise an imaginary tracer of indefinite mass that make their presence felt mainly through the recycling of ammonia [17][39]. Consumer organisms are essentially siphoned off of the standard primary production pathway as a fixed fraction of overall growth [3]. Transfer velocities relating the upper water column to ice fluids and solutes are now parameterized as a function of off-line congelation rates, based on laboratory brine expulsion data [3][17][39]. Thus far we have tested removal from the active bottom zone as a single adjustable time constant, but results have not been formalized or adequately evaluated against data. Internally consistent connections with gravity drainage will be implemented in the vertical biogeochemistry module being developed for future release in CICE.
2.2.1.4. Aerosols¶
Aerosols may be deposited on the ice and gradually work their way through it until the ice melts and they are passed into the ocean. They are defined as ice and snow volume tracers (Equations (26) and (27)), with the snow and ice each having two tracers for each aerosol species, one in the surface scattering layer (delta-Eddington SSL) and one in the snow or ice interior below the SSL.
Rather than updating aerosols for each change to ice/snow thickness due to evaporation, melting, snow–ice formation, etc., during the thermodynamics calculation, these changes are deduced from the diagnostic variables (melts, meltb, snoice, etc) in ice_aerosol.F90. Three processes change the volume of ice or snow but do not change the total amount of aerosol, thus causing the aerosol concentration (the value of the tracer itself) to increase: evaporation, snow deposition and basal ice growth. Basal and lateral melting remove all aerosols in the melted portion. Surface ice and snow melt leave a significant fraction of the aerosols behind, but they do “scavenge” a fraction of them given by the parameter kscav = [0.03, 0.2, 0.02, 0.02, 0.01, 0.01] (only the first 3 are used in CESM, for their 3 aerosol species). Scavenging also applies to snow–ice formation. When sea ice ridges, a fraction of the snow on the ridging ice is thrown into the ocean, and any aerosols in that fraction are also lost to the ocean.
As upper SSL or interior layers disappear from the snow or ice, aerosols are transferred to the next lower layer, or into the ocean when no ice remains. The atmospheric flux faero_atm contains the rates of aerosol deposition for each species, while faero_ocn has the rate at which the aerosols are transferred to the ocean.
The aerosol tracer flag tr_aero must be set to true in ice_in, and the number of aerosol species is set in comp_ice; CESM uses 3. Results using the aerosol code in the context of CESM are documented in [30]. Global diagnostics are available when print_global is true, and history variables include the mass density for each layer (snow and ice SSL and interior), and atmospheric and oceanic fluxes, for each species.
2.2.1.5. Brine height¶
The brine height, \(h_b\), is the distance from the ice–ocean interface to the brine surface. When tr_brine is set true in ice_in and TRBRI is set equal to 1 in comp_ice, the brine surface can move relative to the ice surface. Physically, this occurs when the ice is permeable and there is a nonzero pressure head: the difference between the brine height and the equilibrium sea surface. Brine height motion is computed in ice_brine.F90 from thermodynamic variables (melts, meltb, meltt, etc) and the ice microstructural state deduced from internal bulk salinities and temperature. In the current release, this tracer is for diagnostic purposes only; it is driven by the prognostic salinity parameterization but is not used for computing salinity. In future releases it will be used for transporting biogeochemical tracers vertically through the ice.
Brine height is transported horizontally as the fraction \(f_{bri} = h_b/h_i\), a volume conserved tracer (Equation (26)). Note that unlike the sea ice porosity, brine height fraction may be greater than 1 when \(h_b > h_i\).
Vertical transport processes are, generally, a result of the brine motion. Therefore the vertical transport equations for biogeochemical tracers will be defined only where brine is present. This region, from the ice–ocean interface to the brine height, defines the domain of the vertical bio-grid, whose resolution is independent of the sea ice vertical grid and is set at compile time (see Section Bio-grid). The ice microstructural state, determined in ice_brine.F90, is computed from sea ice salinities and temperatures linearly interpolated to the bio-grid. When \(h_b > h_i\), the upper surface brine is assumed to have the same microstructural properties as the ice surface.
Changes to \(h_b\) occur from ice and snow melt, ice bottom boundary changes, and from pressure adjustments. The computation of \(h_b\) at \(t+\Delta t\) is a two step process. First, \(h_b\) is updated from changes in ice and snow thickness, ie.
Second, pressure driven adjustments arising from meltwater flushing and snow loading are applied to \(h'_b\). Brine flow due to pressure forces are governed by Darcy’s equation
The vertical component of the net permeability tensor \(\Pi^*\) is computed as
where the sea ice is composed of \(N\) vertical layers with \(i\)th layer thickness \(\Delta z_i\) and permeability \(\Pi_i\) ([eq:topo_permea]). The average sea ice density is \(\bar{\rho}\) specified in ice_zbgc_shared.F90. The hydraulic head is \(h_p = h_b - h_{sl}\) where \(h_{sl}\) is the sea level given by
Assuming constant \(h_i\) and \(h_s\) during Darcy flow, the rate of change of \(h_b\) is
where \(w_o = \Pi^* \bar{\rho} g/(h_i\mu\phi_{top})\) and \(\phi_{top}\) is the upper surface porosity. When the Darcy flow is downward into the ice (\(w_o < 0\)), then \(\phi_{top}\) equals the sea ice porosity in the uppermost layer. However, when the flow is upwards into the snow, then \(\phi_{top}\) equals the snow porosity phi_snow specified in ice_in. If a negative number is specified for phi_snow, then the default value is used: phi_snow \(=1 - \rho_s/\rho_w\).
Since \(h_{sl}\) remains relatively unchanged during Darcy flow, Equation (42) has the approximate solution
The contribution \(\Delta h_b|_{h_i,h_s}\) arises from snow and ice melt and bottom ice changes. Since the ice and brine bottom boundaries coincide, changes in the ice bottom from growth or melt, \((\Delta h_i)_{bot}\), equal the bottom brine boundary changes. The surface contribution from ice and snow melt, however, is opposite in sign. The ice contribution is as follows. If \(h_i > h_b\) and the ice surface is melting, ie. \((\Delta h_i)_{top} < 0\)), then melt water increases the brine height:
For snow melt (\(\Delta h_s < 0\)), it is assumed that all snow melt water contributes a source of surface brine. The total change from snow melt and ice thickness changes is
The above brine height calculation is used only when \(h_i\) and \(h_b\) exceed a minimum thickness, thinS, specified in ice_zbgc_shared. Otherwise
provided that \(|h_{sl}-h_b| \leq 0.001\). This formulation ensures small Darcy velocities when \(h_b\) first exceeds thinS.
Both the volume fraction \(f_{bri}\) and the area-weighted brine height \(h_b\) can be written to the history file. The history variable fbri contains the volume-averaged brine fraction tracer
while hbri is comparable to hi (\(h_i\))
where the sums are taken over thickness categories.
2.2.2. Horizontal transport¶
We wish to solve the continuity or transport equation (Equation (20)) for the fractional ice area in each thickness category \(n\). Equation (20) describes the conservation of ice area under horizontal transport. It is obtained from Equation (19) by discretizing \(g\) and neglecting the second and third terms on the right-hand side, which are treated separately (Sections Transport in thickness space and Mechanical redistribution.
There are similar conservation equations for ice volume (Equation (21)), snow volume (Equation (22)), ice energy and snow energy:
By default, ice and snow are assumed to have constant densities, so that volume conservation is equivalent to mass conservation. Variable-density ice and snow layers can be transported conservatively by defining tracers corresponding to ice and snow density, as explained in the introductory comments in ice_transport_remap.F90. Prognostic equations for ice and/or snow density may be included in future model versions but have not yet been implemented.
Two transport schemes are available: upwind and the incremental remapping scheme of [12] as modified for sea ice by [46]. The remapping scheme has several desirable features:
- It conserves the quantity being transported (area, volume, or energy).
- It is non-oscillatory; that is, it does not create spurious ripples in the transported fields.
- It preserves tracer monotonicity. That is, it does not create new extrema in the thickness and enthalpy fields; the values at time \(m+1\) are bounded by the values at time \(m\).
- It is second-order accurate in space and therefore is much less diffusive than first-order schemes (e.g., upwind). The accuracy may be reduced locally to first order to preserve monotonicity.
- It is efficient for large numbers of categories or tracers. Much of the work is geometrical and is performed only once per grid cell instead of being repeated for each quantity being transported.
The time step is limited by the requirement that trajectories projected backward from grid cell corners are confined to the four surrounding cells; this is what is meant by incremental remapping as opposed to general remapping. This requirement leads to a CFL-like condition,
For highly divergent velocity fields the maximum time step must be reduced by a factor of two to ensure that trajectories do not cross. However, ice velocity fields in climate models usually have small divergences per time step relative to the grid size.
The remapping algorithm can be summarized as follows:
- Given mean values of the ice area and tracer fields in each grid cell, construct linear approximations of these fields. Limit the field gradients to preserve monotonicity.
- Given ice velocities at grid cell corners, identify departure regions for the fluxes across each cell edge. Divide these departure regions into triangles and compute the coordinates of the triangle vertices.
- Integrate the area and tracer fields over the departure triangles to obtain the area, volume, and energy transported across each cell edge.
- Given these transports, update the state variables.
Since all scalar fields are transported by the same velocity field, step (2) is done only once per time step. The other three steps are repeated for each field in each thickness category. These steps are described below.
2.2.2.1. Reconstructing area and tracer fields¶
First, using the known values of the state variables, the ice area and tracer fields are reconstructed in each grid cell as linear functions of \(x\) and \(y\). For each field we compute the value at the cell center (i.e., at the origin of a 2D Cartesian coordinate system defined for that grid cell), along with gradients in the \(x\) and \(y\) directions. The gradients are limited to preserve monotonicity. When integrated over a grid cell, the reconstructed fields must have mean values equal to the known state variables, denoted by \(\bar{a}\) for fractional area, \(\tilde{h}\) for thickness, and \(\hat{q}\) for enthalpy. The mean values are not, in general, equal to the values at the cell center. For example, the mean ice area must equal the value at the centroid, which may not lie at the cell center.
Consider first the fractional ice area, the analog to fluid density \(\rho\) in [12]. For each thickness category we construct a field \(a({\bf r})\) whose mean is \(\bar{a}\), where \({\bf r} = (x,y)\) is the position vector relative to the cell center. That is, we require
where \(A=\int_A dA\) is the grid cell area. Equation (48) is satisfied if \(a({\bf r})\) has the form
where \(\left<\nabla a\right>\) is a centered estimate of the area gradient within the cell, \(\alpha_a\) is a limiting coefficient that enforces monotonicity, and \({\bf \bar{r}}\) is the cell centroid:
It follows from Equation (49) that the ice area at the cell center (\(\mathbf{r} = 0\)) is
where \(a_x = \alpha_a (\partial a / \partial x)\) and \(a_y = \alpha_a (\partial a / \partial y)\) are the limited gradients in the \(x\) and \(y\) directions, respectively, and the components of \({\bf \bar{r}}\), \(\overline{x} = \int_A x \, dA / A\) and \(\overline{y} = \int_A y \, dA / A\), are evaluated using the triangle integration formulas described in Section Integrating fields. These means, along with higher-order means such as \(\overline{x^2}\), \(\overline{xy}\), and \(\overline{y^2}\), are computed once and stored.
Next consider the ice and snow thickness and enthalpy fields. Thickness is analogous to the tracer concentration \(T\) in [12], but there is no analog in [12] to the enthalpy. The reconstructed ice or snow thickness \(h({\bf r})\) and enthalpy \(q(\mathbf{r})\) must satisfy
where \(\tilde{h}=h(\tilde{\bf r})\) is the thickness at the center of ice area, and \(\hat{q}=q(\hat{\bf r})\) is the enthalpy at the center of ice or snow volume. Equations ([eq:mean_thickness]) and ([eq:mean_enthalpy]) are satisfied when \(h({\bf r})\) and \(q({\bf r})\) are given by
where \(\alpha_h\) and \(\alpha_q\) are limiting coefficients. The center of ice area, \({\bf\tilde{r}}\), and the center of ice or snow volume, \({\bf \hat{r}}\), are given by
Evaluating the integrals, we find that the components of \({\bf \tilde{r}}\) are
and the components of \({\bf \hat{r}}\) are
where
From Equation (52) and Equation (53), the thickness and enthalpy at the cell center are given by
where \(h_x\), \(h_y\), \(q_x\) and \(q_y\) are the limited gradients of thickness and enthalpy. The surface temperature is treated the same way as ice or snow thickness, but it has no associated enthalpy. Tracers obeying conservation equations of the form Equation (26) and Equation (27) are treated in analogy to ice and snow enthalpy, respectively.
We preserve monotonicity by van Leer limiting. If \(\bar{\phi}(i,j)\) denotes the mean value of some field in grid cell \((i,j)\), we first compute centered gradients of \(\bar{\phi}\) in the \(x\) and \(y\) directions, then check whether these gradients give values of \(\phi\) within cell \((i,j)\) that lie outside the range of \(\bar{\phi}\) in the cell and its eight neighbors. Let \(\bar{\phi}_{\max}\) and \(\bar{\phi}_{\min}\) be the maximum and minimum values of \(\bar{\phi}\) over the cell and its neighbors, and let \(\phi_{\max}\) and \(\phi_{\min}\) be the maximum and minimum values of the reconstructed \(\phi\) within the cell. Since the reconstruction is linear, \(\phi_{\max}\) and \(\phi_{\min}\) are located at cell corners. If \(\phi_{\max} > \bar{\phi}_{\max}\) or :math:` phi_{min} < bar{phi}_{min}`, we multiply the unlimited gradient by \(\alpha = \min(\alpha_{\max}, \alpha_{\min})\), where
Otherwise the gradient need not be limited.
Earlier versions of CICE (through 3.14) computed gradients in physical space. Starting in version 4.0, gradients are computed in a scaled space in which each grid cell has sides of unit length. The origin is at the cell center, and the four vertices are located at (0.5, 0.5), (-0.5,0.5),(-0.5, -0.5) and (0.5, -0.5). In this coordinate system, several of the above grid-cell-mean quantities vanish (because they are odd functions of x and/or y), but they have been retained in the code for generality.
2.2.2.2. Locating departure triangles¶
The method for locating departure triangles is discussed in detail by [12]. The basic idea is illustrated in Figure 2, which shows a shaded quadrilateral departure region whose contents are transported to the target or home grid cell, labeled \(H\). The neighboring grid cells are labeled by compass directions: \(NW\), \(N\), \(NE\), \(W\), and \(E\). The four vectors point along the velocity field at the cell corners, and the departure region is formed by joining the starting points of these vectors. Instead of integrating over the entire departure region, it is convenient to compute fluxes across cell edges. We identify departure regions for the north and east edges of each cell, which are also the south and west edges of neighboring cells. Consider the north edge of the home cell, across which there are fluxes from the neighboring \(NW\) and \(N\) cells. The contributing region from the \(NW\) cell is a triangle with vertices \(abc\), and that from the \(N\) cell is a quadrilateral that can be divided into two triangles with vertices \(acd\) and \(ade\). Focusing on triangle \(abc\), we first determine the coordinates of vertices \(b\) and \(c\) relative to the cell corner (vertex \(a\)), using Euclidean geometry to find vertex \(c\). Then we translate the three vertices to a coordinate system centered in the \(NW\) cell. This translation is needed in order to integrate fields (Section Integrating fields) in the coordinate system where they have been reconstructed (Section Reconstructing area and tracer fields). Repeating this process for the north and east edges of each grid cell, we compute the vertices of all the departure triangles associated with each cell edge.
Figure 2
Figure 2 : In incremental remapping, conserved quantities are remapped from the shaded departure region, a quadrilateral formed by connecting the backward trajectories from the four cell corners, to the grid cell labeled \(H\). The region fluxed across the north edge of cell \(H\) consists of a triangle (\(abc\)) in the \(NW\) cell and a quadrilateral (two triangles, \(acd\) and \(ade\)) in the \(N\) cell.
Figure 3, reproduced from [12], shows all possible triangles that can contribute fluxes across the north edge of a grid cell. There are 20 triangles, which can be organized into five groups of four mutually exclusive triangles as shown in Table 3. In this table, \((x_1, y_1)\) and \((x_2,y_2)\) are the Cartesian coordinates of the departure points relative to the northwest and northeast cell corners, respectively. The departure points are joined by a straight line that intersects the west edge at \((0,y_a)\) relative to the northwest corner and intersects the east edge at \((0,y_b)\) relative to the northeast corner. The east cell triangles and selecting conditions are identical except for a rotation through 90 degrees.
Figure 3
Figure 3 : The 20 possible triangles that can contribute fluxes across the north edge of a grid cell.
Table 3 : Evaluation of contributions from the 20 triangles across the north cell edge. The coordinates \(x_1\), \(x_2\), \(y_1\), \(y_2\), \(y_a\), and \(y_b\) are defined in the text. We define \(\tilde{y}_1 = y_1\) if \(x_1>0\), else \(\tilde{y}_1 = y_a\). Similarly, \(\tilde{y}_2 = y_2\) if \(x_2<0\), else \(\tilde{y}_2 = y_b\).
| Triangle | Triangle | Selecting logical | |
| group | label | condition | |
| 1 | NW | \(y_a>0\) and \(y_1\geq0\) and \(x_1<0\) | |
| NW1 | \(y_a<0\) and \(y_1\geq0\) and \(x_1<0\) | ||
| W | \(y_a<0\) and \(y_1<0\) and \(x_1<0\) | ||
| W2 | \(y_a>0\) and \(y_1<0\) and \(x_1<0\) | ||
| 2 | NE | \(y_b>0\) and \(y_2\geq0\) and \(x_2>0\) | |
| NE1 | \(y_b<0\) and \(y_2\geq0\) and \(x_2>0\) | ||
| E | \(y_b<0\) and \(y_2<0\) and \(x_2>0\) | ||
| E2 | \(y_b>0\) and \(y_2<0\) and \(x_2>0\) | ||
| 3 | W1 | \(y_a<0\) and \(y_1\geq0\) and \(x_1<0\) | |
| NW2 | \(y_a>0\) and \(y_1<0\) and \(x_1<0\) | ||
| E1 | \(y_b<0\) and \(y_2\geq0\) and \(x_2>0\) | ||
| NE2 | \(y_b>0\) and \(y_2<0\) and \(x_2>0\) | ||
| 4 | H1a | \(y_a y_b\geq 0\) and \(y_a+y_b<0\) | |
| N1a | \(y_a y_b\geq 0\) and \(y_a+y_b>0\) | ||
| H1b | \(y_a y_b<0\) and \(\tilde{y}_1<0\) | ||
| N1b | \(y_a y_b<0\) and \(\tilde{y}_1>0\) | ||
| 5 | H2a | \(y_a y_b\geq 0\) and \(y_a+y_b<0\) | |
| N2a | \(y_a y_b\geq 0\) and \(y_a+y_b>0\) | ||
| H2b | \(y_a y_b<0\) and \(\tilde{y}_2<0\) | ||
| N2b | \(y_a y_b<0\) and \(\tilde{y}_2>0\) | ||
This scheme was originally designed for rectangular grids. Grid cells in CICE actually lie on the surface of a sphere and must be projected onto a plane. The projection used in CICE maps each grid cell to a square with sides of unit length. Departure triangles across a given cell edge are computed in a coordinate system whose origin lies at the midpoint of the edge and whose vertices are at (-0.5, 0) and (0.5, 0). Intersection points are computed assuming Cartesian geometry with cell edges meeting at right angles. Let CL and CR denote the left and right vertices, which are joined by line CLR. Similarly, let DL and DR denote the departure points, which are joined by line DLR. Also, let IL and IR denote the intersection points (0, \(y_a\)) and (0,:math:y_b) respectively, and let IC = (\(x_c\), 0) denote the intersection of CLR and DLR. It can be shown that \(y_a\), \(y_b\), and \(x_c\) are given by
Each departure triangle is defined by three of the seven points (CL, CR, DL, DR, IL, IR, IC).
Given a 2D velocity field u, the divergence \(\nabla\cdot{\bf u}\) in a given grid cell can be computed from the local velocities and written in terms of fluxes across each cell edge:
where \(L\) is an edge length and the indices \(N, S, E, W\) denote compass directions. Equation (54) is equivalent to the divergence computed in the EVP dynamics (Section Dynamics). In general, the fluxes in this expression are not equal to those implied by the above scheme for locating departure regions. For some applications it may be desirable to prescribe the divergence by prescribing the area of the departure region for each edge. This can be done in CICE 4.0 by setting l_fixed_area = true in ice_transport_driver.F90 and passing the prescribed departure areas (edgearea_e and edgearea_n) into the remapping routine. An extra triangle is then constructed for each departure region to ensure that the total area is equal to the prescribed value. This idea was suggested and first implemented by Mats Bentsen of the Nansen Environmental and Remote Sensing Center (Norway), who applied an earlier version of the CICE remapping scheme to an ocean model. The implementation in CICE 4.0 is somewhat more general, allowing for departure regions lying on both sides of a cell edge. The extra triangle is constrained to lie in one but not both of the grid cells that share the edge. Since this option has yet to be fully tested in CICE, the current default is l_fixed_area = false.
We made one other change in the scheme of [12] for locating triangles. In their paper, departure points are defined by projecting cell corner velocities directly backward. That is,
where \(\mathbf{x}_D\) is the location of the departure point relative to the cell corner and \(\mathbf{u}\) is the velocity at the corner. This approximation is only first-order accurate. Accuracy can be improved by estimating the velocity at the midpoint of the trajectory.
2.2.2.3. Integrating fields¶
Next, we integrate the reconstructed fields over the departure triangles to find the total area, volume, and energy transported across each cell edge. Area transports are easy to compute since the area is linear in \(x\) and \(y\). Given a triangle with vertices \(\mathbf{x_i} = (x_i,y_i)\), \(i\in\{1,2,3\}\), the triangle area is
The integral \(F_a\) of any linear function \(f(\mathbf{r})\) over a triangle is given by
where \(\mathbf{x}_0 = (x_0,y_0)\) is the triangle midpoint,
To compute the area transport, we evaluate the area at the midpoint,
and multiply by \(A_T\). By convention, northward and eastward transport is positive, while southward and westward transport is negative.
Equation (56) cannot be used for volume transport, because the reconstructed volumes are quadratic functions of position. (They are products of two linear functions, area and thickness.) The integral of a quadratic polynomial over a triangle requires function evaluations at three points,
where \(\mathbf{x}_i^\prime = (\mathbf{x}_0+\mathbf{x}_i)/2\) are points lying halfway between the midpoint and the three vertices. [12] use this formula to compute transports of the product \(\rho \, T\), which is analogous to ice volume. Equation (57) does not work for ice and snow energies, which are cubic functions—products of area, thickness, and enthalpy. Integrals of a cubic polynomial over a triangle can be evaluated using a four-point formula [71]:
where \(\mathbf{x_i}^{\prime\prime}=(3 \mathbf{x}_0 + 2 \mathbf{x}_i)/5\). To evaluate functions at specific points, we must compute many products of the form \(a({\bf x}) \, h({\bf x})\) and \(a({\bf x}) \, h({\bf x}) \, q({\bf x})\), where each term in the product is the sum of a cell-center value and two displacement terms. In the code, the computation is sped up by storing some sums that are used repeatedly.
2.2.2.4. Updating state variables¶
Finally, we compute new values of the state variables in each ice category and grid cell. The new fractional ice areas \(a_{in}^\prime(i,j)\) are given by
where \(F_{aE}(i,j)\) and \(F_{aN}(i,j)\) are the area transports across the east and north edges, respectively, of cell \((i,j)\), and \(A(i,j)\) is the grid cell area. All transports added to one cell are subtracted from a neighboring cell; thus Equation (59) conserves total ice area.
The new ice volumes and energies are computed analogously. New thicknesses are given by the ratio of volume to area, and enthalpies by the ratio of energy to volume. Tracer monotonicity is ensured because
where \(h^\prime\) and \(q^\prime\) are the new-time thickness and enthalpy, given by integrating the old-time ice area, volume, and energy over a Lagrangian departure region with area \(A\). That is, the new-time thickness and enthalpy are weighted averages over old-time values, with non-negative weights \(a\) and \(ah\). Thus the new-time values must lie between the maximum and minimum of the old-time values.
2.2.3. Dynamics¶
There are now different rheologies available in the CICE code. The elastic-viscous-plastic (EVP) model represents a modification of the standard viscous-plastic (VP) model for sea ice dynamics [26]. The elastic-anisotropic-plastic (EAP) model, on the other hand, explicitly accounts for the observed sub-continuum anisotropy of the sea ice cover [83][81]. If kdyn = 1 in the namelist then the EVP rheology is used (module ice_dyn_evp.F90), while kdyn = 2 is associated with the EAP rheology (ice_dyn_eap.F90). At times scales associated with the wind forcing, the EVP model reduces to the VP model while the EAP model reduces to the anisotropic rheology described in detail in [83][77]. At shorter time scales the adjustment process takes place in both models by a numerically more efficient elastic wave mechanism. While retaining the essential physics, this elastic wave modification leads to a fully explicit numerical scheme which greatly improves the model’s computational efficiency.
The EVP sea ice dynamics model is thoroughly documented in [33], [31], [34] and [35] and the EAP dynamics in [77]. Simulation results and performance of the EVP and EAP models have been compared with the VP model and with each other in realistic simulations of the Arctic respectively in [37] and [77]. Here we summarize the equations and direct the reader to the above references for details. The numerical implementation in this code release is that of [34] and [35], with revisions to the numerical solver as in [8]. The implementation of the EAP sea ice dynamics into CICE is described in detail in [77].
2.2.3.1. Momentum¶
The force balance per unit area in the ice pack is given by a two-dimensional momentum equation [26], obtained by integrating the 3D equation through the thickness of the ice in the vertical direction:
where \(m\) is the combined mass of ice and snow per unit area and \(\vec{\tau}_a\) and \(\vec{\tau}_w\) are wind and ocean stresses, respectively. The strength of the ice is represented by the internal stress tensor \(\sigma_{ij}\), and the other two terms on the right hand side are stresses due to Coriolis effects and the sea surface slope. The parameterization for the wind and ice–ocean stress terms must contain the ice concentration as a multiplicative factor to be consistent with the formal theory of free drift in low ice concentration regions. A careful explanation of the issue and its continuum solution is provided in [35] and :cite:CGHM04`.
The momentum equation is discretized in time as follows, for the classic EVP approach. First, for clarity, the two components of Equation (60) are
In the code, \({\tt vrel}=a_i c_w \rho_w\left|{\bf U}_w - {\bf u}^k\right|\), where \(k\) denotes the subcycling step. The following equations illustrate the time discretization and define some of the other variables used in the code.
and vrel\(\cdot\)waterx(y) = taux(y).
We solve this system of equations analytically for \(u^{k+1}\) and \(v^{k+1}\). Define
where \({\bf F} = \nabla\cdot\sigma^{k+1}\). Then
Solving simultaneously for \(u^{k+1}\) and \(v^{k+1}\),
where
When the subcycling is finished for each (thermodynamic) time step, the ice–ocean stress must be constructed from taux(y) and the terms containing vrel on the left hand side of the equations.
The Hibler-Bryan form for the ice-ocean stress [28] is included in ice_dyn_shared.F90 but is currently commented out, pending further testing.
2.2.3.2. Internal stress¶
For convenience we formulate the stress tensor \(\bf \sigma\) in terms of \(\sigma_1=\sigma_{11}+\sigma_{22}\), \(\sigma_2=\sigma_{11}-\sigma_{22}\), and introduce the divergence, \(D_D\), and the horizontal tension and shearing strain rates, \(D_T\) and \(D_S\) respectively.
Elastic-Viscous-Plastic
In the EVP model the internal stress tensor is determined from a regularized version of the VP constitutive law,
where
and \(P\) is a function of the ice thickness and concentration, described in Section Mechanical redistribution. The dynamics component employs a “replacement pressure” (see [24], for example), which serves to prevent residual ice motion due to spatial variations of \(P\) when the rates of strain are exactly zero.
Viscosities are updated during the subcycling, so that the entire dynamics component is subcycled within the time step, and the elastic parameter \(E\) is defined in terms of a damping timescale \(T\) for elastic waves, \(\Delta t_e < T < \Delta t\), as
where \(T=E_\circ\Delta t\) and \(E_\circ\) (eyc) is a tunable parameter less than one. The stress equations (67)–(69) become
All coefficients on the left-hand side are constant except for \(P\), which changes only on the longer time step \(\Delta t\). This modification compensates for the decreased efficiency of including the viscosity terms in the subcycling. (Note that the viscosities do not appear explicitly.) Choices of the parameters used to define \(E\), \(T\) and \(\Delta t_e\) are discussed in Sections Revised approach and Choosing an appropriate time step.
The bilinear discretization used for the stress terms \(\partial\sigma_{ij}/\partial x_j\) in the momentum equation is now used, which enabled the discrete equations to be derived from the continuous equations written in curvilinear coordinates. In this manner, metric terms associated with the curvature of the grid are incorporated into the discretization explicitly. Details pertaining to the spatial discretization are found in [34].
Elastic-Anisotropic-Plastic
In the EAP model the internal stress tensor is related to the geometrical properties and orientation of underlying virtual diamond shaped floes (see Figure 5). In contrast to the isotropic EVP rheology, the anisotropic plastic yield curve within the EAP rheology depends on the relative orientation of the diamond shaped floes (unit vector \(\mathbf r\) in Figure 5), with respect to the principal direction of the deformation rate (not shown). Local anisotropy of the sea ice cover is accounted for by an additional prognostic variable, the structure tensor \(\mathbf{A}\) defined by
where \(\mathbb{S}\) is a unit-radius circle; A is a unit trace, 2\(\times\)2 matrix. From now on we shall describe the orientational distribution of floes using the structure tensor. For simplicity we take the probability density function \(\vartheta(\mathbf r )\) to be Gaussian, \(\vartheta(z)=\omega_{1}\exp(-\omega_{2}z^{2})\), where \(z\) is the ice floe inclination with respect to the axis \(x_{1}\) of preferential alignment of ice floes (see Figure 5), \(\vartheta(z)\) is periodic with period \(\pi\), and the positive coefficients \(\omega_{1}\) and \(\omega_{2}\) are calculated to ensure normalization of \(\vartheta(z)\), i.e. \(\int_{0}^{2\pi}\vartheta(z)dz=1\). The ratio of the principal components of \(\mathbf{A}\), \(A_{1}/A_{2}\), are derived from the phenomenological evolution equation for the structure tensor \(\mathbf A\),
where \(t\) is the time, and \(D/Dt\) is the co-rotational time derivative accounting for advection and rigid body rotation (\(D\mathbf A/Dt = d\mathbf A/dt -\mathbf W \cdot \mathbf A -\mathbf A \cdot \mathbf W^{T}\)) with \(\mathbf W\) being the vorticity tensor. \(\mathbf F_{iso}\) is a function that accounts for a variety of processes (thermal cracking, melting, freezing together of floes) that contribute to a more isotropic nature to the ice cover. \(\mathbf F_{frac}\) is a function determining the ice floe re-orientation due to fracture, and explicitly depends upon sea ice stress (but not its magnitude). Following [83], based on laboratory experiments by [62] we consider four failure mechanisms for the Arctic sea ice cover. These are determined by the ratio of the principal values of the sea ice stress \(\sigma_{1}\) and \(\sigma_{2}\): (i) under biaxial tension, fractures form across the perpendicular principal axes and therefore counteract any apparent redistribution of the floe orientation; (ii) if only one of the principal stresses is compressive, failure occurs through axial splitting along the compression direction; (iii) under biaxial compression with a low confinement ratio, (\(\sigma_{1}/\sigma_{2}<R\)), sea ice fails Coulombically through formation of slip lines delineating new ice floes oriented along the largest compressive stress; and finally (iv) under biaxial compression with a large confinement ratio, (\(\sigma_{1}/\sigma_{2}\ge R\)), the ice is expected to fail along both principal directions so that the cumulative directional effect balances to zero.
Figure 5
Figure 5 : Geometry of interlocking diamond-shaped floes (taken from [83]). \(\phi\) is half of the acute angle of the diamonds. \(L\) is the edge length. \(\boldsymbol n_{1}\), \(\boldsymbol n_{2}\) and \(\boldsymbol\tau_{1}\), \(\boldsymbol\tau_{2}\) are respectively the normal and tangential unit vectors along the diamond edges. \(\mathbf v=L\boldsymbol\tau_{2}\cdot\dot{\boldsymbol\epsilon}\) is the relative velocity between the two floes connected by the vector \(L \boldsymbol \tau_{2}\). \(\mathbf r\) is the unit vector along the main diagonal of the diamond. Note that the diamonds illustrated here represent one possible realisation of all possible orientations. The angle \(z\) represents the rotation of the diamonds’ main axis relative to their preferential orientation along the axis \(x_1\).
The new anisotropic rheology requires solving the evolution Equation (70) for the structure tensor in addition to the momentum and stress equations. The evolution equation for \(\mathbf{A}\) is solved within the EVP subcycling loop, and consistently with the momentum and stress evolution equations, we neglect the advection term for the structure tensor. Equation (70) then reduces to the system of two equations:
where the first terms on the right hand side correspond to the isotropic contribution, \(F_{iso}\), and \(M_{11}\) and \(M_{12}\) are the components of the term \(F_{frac}\) in Equation (70) that are given in [83] and [77]. These evolution equations are discretized semi-implicitly in time. The degree of anisotropy is measured by the largest eigenvalue (\(A_{1}\)) of this tensor (\(A_{2}=1-A_{1}\)). \(A_{1}=1\) corresponds to perfectly aligned floes and \(A_{1}=0.5\) to a uniform distribution of floe orientation. Note that while we have specified the aspect ratio of the diamond floes, through prescribing \(\phi\), we make no assumption about the size of the diamonds so that formally the theory is scale invariant.
As described in greater detail in [83], the internal ice stress for a single orientation of the ice floes can be calculated explicitly and decomposed, for an average ice thickness \(h\), into its ridging (r) and sliding (s) contributions
where \(P_{r}\) and \(P_{s}\) are the ridging and sliding strengths and the ridging and sliding stresses are functions of the angle \(\theta= \arctan(\dot\epsilon_{II}/\dot\epsilon_{I})\), the angle \(y\) between the major principal axis of the strain rate tensor (not shown) and the structure tensor (\(x_1\) axis in Figure 5, and the angle \(z\) defined in Figure 5. In the stress expressions above the underlying floes are assumed parallel, but in a continuum-scale sea ice region the floes can possess different orientations in different places and we take the mean sea ice stress over a collection of floes to be given by the average
where we have introduced the friction parameter \(k=P_{s}/P_{r}\) and where we identify the ridging ice strength \(P_{r}(h)\) with the strength \(P\) described in section 1 and used within the EVP framework.
As is the case for the EVP rheology, elasticity is included in the EAP description not to describe any physical effect, but to make use of the efficient, explicit numerical algorithm used to solve the full sea ice momentum balance. We use the analogous EAP stress equations,
where the anisotropic stress \(\boldsymbol\sigma^{EAP}\) is defined in a look-up table for the current values of strain rate and structure tensor. The look-up table is constructed by computing the stress (normalized by the strength) from Equations (73)–(75) for discrete values of the largest eigenvalue of the structure tensor, \(\frac{1}{2}\le A_{1}\le 1\), the angle \(0\le\theta\le2\pi\), and the angle \(-\pi/2\le y\le\pi/2\) between the major principal axis of the strain rate tensor and the structure tensor [77]. The updated stress, after the elastic relaxation, is then passed to the momentum equation and the sea ice velocities are updated in the usual manner within the subcycling loop of the EVP rheology. The structure tensor evolution equations are solved implicitly at the same frequency, \(\Delta t_{e}\), as the ice velocities and internal stresses. Finally, to be coherent with our new rheology we compute the area loss rate due to ridging as \(\vert\dot{\boldsymbol\epsilon}\vert\alpha_{r}(\theta)\), with \(\alpha_r(\theta)\) and \(\alpha_s(\theta)\) given by [82],
Both ridging rate and sea ice strength are computed in the outer loop of the dynamics.
2.2.3.3. Revised approach¶
A modification of the standard elastic-viscous-plastic (EVP) approach for sea ice dynamics has been proposed by [8], that generalizes the EVP elastic modulus \(E\) and the time stepping approach for both momentum and stress to use an under-relaxation technique. In general terms, the momentum and stress equations become
where \(\breve{{\bf u}}\) and \(\breve{\sigma}\) represent the converged VP solution and \(\alpha, \beta < 1\).
Momentum
The momentum equations become
where
\({\bf u}^\circ\) is the initial value of velocity at the beginning of the subcycling (\(k=0\)), and we use \({\bf u}^{k+1}\) for the ice–ocean stress and Coriolis terms. Equations (76) and (77) differ from Equations (63) and (64) only in the last term.
Solving simultaneously for \({\bf u}^{k+1}\) as before, we have
where
and \(b\) is the same as in Equation (66).
Stress
In CICE’s classic approach, the update to \(\sigma_1\) at subcycle step \(k+1\) is
If we set
then Equation (79) becomes
This is equivalent to Eq. (23) in [8], but using \(\sigma\) at the current subcycle \(k+1\) in the last term on the right-hand side. Likewise, setting
produces equations equivalent to Eq. (23) in [8] for \(\sigma_2\) and \(\sigma_{12}\). Therefore the only change needed in the stress code is to use \(\alpha_1\) and \(\alpha_2\) instead of \(2T / \Delta t_e\) and \(2T /e^2 \Delta t_e\).
However, [8] introduce another change to the EVP stress equations by altering the form of Young’s modulus in the elastic term: the coefficient of \(\partial\sigma_1/\partial t\) is \(1/E\), but it is \(e^2/E\) in the \(\sigma_2\) and \(\sigma_{12}\) equations. This change does not affect the VP equations to which the EVP equations should converge, but it does affect the transient path taken during the subcycling. Since EVP subcycling is finite, the numerical solutions obtained using this method differ from the original EVP code.
To implement this second change, we need define only \(\alpha_1 = {2T/\Delta t_e}\) as above and incorporate the factor of \(e^2\) from \(\alpha_2\) into the equations for \(\sigma_2\) and \(\sigma_{12}\):
To minimize code changes and unify the two approaches, we define and apply \(1/\alpha_1\) and \(\beta\) in the classic EVP code, and modify the elastic stress term. These under-relaxation parameters control the rate at which the iteration converges. Thus for classic EVP we set
Then
The stress equations for stressp (\(\sigma_1\)) are unchanged; the modified equations for stressm (\(\sigma_2\)) and stress12 (\(\sigma_{12}\)) take the form
For classic EVP,
For revised EVP, arlx1i and brlx are defined separately from \(\Delta t\), \(\Delta t_e\), \(T\) and \(e\), and
\(\tilde{\bf u}\) must also be defined for revised EVP as in Equation (78). The extra terms in \(\tilde{a}\) and \(\tilde{\bf u}\) are multiplied by a flag (revp) that equals 1 for revised EVP and 0 for classic EVP. Revised EVP is activated by setting the namelist parameter revised_evp = true. Note that in the current implementation, only the modified version of the elastic term is available for either the classic (revised_evp = false) or the revised EVP method. A final difference is that the revised approach initializes the stresses to 0 at the beginning of each time step, while the classic EVP approach uses the previous time step value.
2.2.4. Thickness changes¶
2.2.4.1. Transport in thickness space¶
Next we solve the equation for ice transport in thickness space due to thermodynamic growth and melt,
which is obtained from Equation (19) by neglecting the first and third terms on the right-hand side. We use the remapping method of [45], in which thickness categories are represented as Lagrangian grid cells whose boundaries are projected forward in time. The thickness distribution function \(g\) is approximated as a linear function of \(h\) in each displaced category and is then remapped onto the original thickness categories. This method is numerically smooth and is not too diffusive. It can be viewed as a 1D simplification of the 2D incremental remapping scheme described above.
We first compute the displacement of category boundaries in thickness space. Assume that at time \(m\) the ice areas \(a_n^m\) and mean ice thicknesses \(h_n^m\) are known for each thickness category. (For now we omit the subscript \(i\) that distinguishes ice from snow.) We use a thermodynamic model (Section Thermodynamics) to compute the new mean thicknesses \(h_n^{m+1}\) at time \(m+1\). The time step must be small enough that trajectories do not cross; i.e., \(h_n^{m+1} < h_{n+1}^{m+1}\) for each pair of adjacent categories. The growth rate at \(h = h_n\) is given by \(f_n = (h_n^{m+1} - h_n^m) / \Delta t\). By linear interpolation we estimate the growth rate \(F_n\) at the upper category boundary \(H_n\):
If \(a_n\) or \(a_{n+1} = 0\), \(F_n\) is set to the growth rate in the nonzero category, and if \(a_n = a_{n+1} = 0\), we set \(F_n = 0\). The temporary displaced boundaries are given by
The boundaries must not be displaced by more than one category to the left or right; that is, we require \(H_{n-1} < H_n^* < H_{n+1}\). Without this requirement we would need to do a general remapping rather than an incremental remapping, at the cost of added complexity.
Next we construct \(g(h)\) in the displaced thickness categories. The ice areas in the displaced categories are \(a_n^{m+1} = a_n^m\), since area is conserved following the motion in thickness space (i.e., during vertical ice growth or melting). The new ice volumes are \(v_n^{m+1} = (a_n h_n)^{m+1} = a_n^m h_n^{m+1}\). For conciseness, define \(H_L = H_{n-1}^*\) and \(H_R = H_{n}^*\) and drop the time index \(m+1\). We wish to construct a continuous function \(g(h)\) within each category such that the total area and volume at time \(m+1\) are \(a_n\) and \(v_n\), respectively:
The simplest polynomial that can satisfy both equations is a line. It is convenient to change coordinates, writing \(g(\eta) = g_1 \eta + g_0\), where \(\eta = h - H_L\) and the coefficients \(g_0\) and \(g_1\) are to be determined. Then Equations (81) and (82) can be written as
where \(\eta_R = H_R - H_L\) and \(\eta_n = h_n - H_L\). These equations have the solution
Since \(g\) is linear, its maximum and minimum values lie at the boundaries, \(\eta = 0\) and \(\eta_R\):
Equation (85) implies that \(g(0) < 0\) when \(\eta_n > 2 \eta_R/3\), i.e., when \(h_n\) lies in the right third of the thickness range \((H_L, H_R)\). Similarly, Equation (86) implies that \(g(\eta_R) < 0\) when \(\eta_n < \eta_R/3\), i.e., when \(h_n\) is in the left third of the range. Since negative values of \(g\) are unphysical, a different solution is needed when \(h_n\) lies outside the central third of the thickness range. If \(h_n\) is in the left third of the range, we define a cutoff thickness, \(H_C = 3 h_n - 2 H_L\), and set \(g = 0\) between \(H_C\) and \(H_R\). Equations (83) and (84) are then valid with \(\eta_R\) redefined as \(H_C - H_L\). And if \(h_n\) is in the right third of the range, we define \(H_C = 3 h_n - 2 H_R\) and set \(g = 0\) between \(H_L\) and \(H_C\). In this case, (83) and (84) apply with \(\eta_R = H_R - H_C\) and \(\eta_n = h_n - H_C\).
Figure 4 illustrates the linear reconstruction of \(g\) for the simple cases \(H_L = 0\), \(H_R = 1\), \(a_n = 1\), and \(h_n =\) 0.2, 0.4, 0.6, and 0.8. Note that \(g\) slopes downward (\(g_1 < 0\)) when \(h_n\) is less than the midpoint thickness, \((H_L + H_R)/2 = 1/2\), and upward when \(h_n\) exceeds the midpoint thickness. For \(h_n = 0.2\) and 0.8, \(g = 0\) over part of the range.
Figure 4
Figure 4 : Linear approximation of the thickness distribution function \(g(h)\) for an ice category with left boundary \(H_L = 0\), right boundary \(H_R = 1\), fractional area \(a_n = 1\), and mean ice thickness \(h_n = 0.2, 0.4, 0.6,\) and \(0.8\).
Finally, we remap the thickness distribution to the original boundaries by transferring area and volume between categories. We compute the ice area \(\Delta a_n\) and volume \(\Delta v_n\) between each original boundary \(H_n\) and displaced boundary \(H_n^*\). If \(H_n^* > H_n\), ice moves from category \(n\) to \(n+1\). The area and volume transferred are
If \(H_n^* < H_N\), ice area and volume are transferred from category \(n+1\) to \(n\) using Equations (87) and (88) with the limits of integration reversed. To evaluate the integrals we change coordinates from \(h\) to \(\eta = h - H_L\), where \(H_L\) is the left limit of the range over which \(g > 0\), and write \(g(\eta)\) using Equations (83) and (84). In this way we obtain the new areas \(a_n\) and volumes \(v_n\) between the original boundaries \(H_{n-1}\) and \(H_n\) in each category. The new thicknesses, \(h_n = v_n/a_n\), are guaranteed to lie in the range \((H_{n-1}, H_n)\). If \(g = 0\) in the part of a category that is remapped to a neighboring category, no ice is transferred.
Other conserved quantities are transferred in proportion to the ice volume \(\Delta v_{in}\). For example, the transferred ice energy in layer \(k\) is \(\Delta e_{ink} = e_{ink} (\Delta v_{in} / v_{in})\).
The left and right boundaries of the domain require special treatment. If ice is growing in open water at a rate \(F_0\), the left boundary \(H_0\) is shifted to the right by \(F_0 \Delta t\) before \(g\) is constructed in category 1, then reset to zero after the remapping is complete. New ice is then added to the grid cell, conserving area, volume, and energy. If ice cannot grow in open water (because the ocean is too warm or the net surface energy flux is downward), \(H_0\) is fixed at zero, and the growth rate at the left boundary is estimated as \(F_0 = f_1\). If \(F_0 < 0\), all ice thinner than \(\Delta h_0 = -F_0 \Delta t\) is assumed to have melted, and the ice area in category 1 is reduced accordingly. The area of new open water is
The right boundary \(H_N\) is not fixed but varies with \(h_N\), the mean ice thickness in the thickest category. Given \(h_N\), we set \(H_N = 3 h_N - 2 H_{N-1}\), which ensures that \(g(h) > 0\) for \(H_{N-1} < h < H_N\) and \(g(h) = 0\) for \(h \geq H_N\). No ice crosses the right boundary. If the ice growth or melt rates in a given grid cell are too large, the thickness remapping scheme will not work. Instead, the thickness categories in that grid cell are treated as delta functions following [6], and categories outside their prescribed boundaries are merged with neighboring categories as needed. For time steps of less than a day and category thickness ranges of 10 cm or more, this simplification is needed rarely, if ever.
The linear remapping algorithm for thickness is not monotonic for tracers, although significant errors rarely occur. Usually they appear as snow temperatures (enthalpy) outside the physical range of values in very small snow volumes. In this case we transfer the snow and its heat and tracer contents to the ocean.
2.2.4.2. Mechanical redistribution¶
The last term on the right-hand side of Equation (19) is \(\psi\), which describes the redistribution of ice in thickness space due to ridging and other mechanical processes. The mechanical redistribution scheme in CICE is based on [74], [61], [27], [19], and [47]. This scheme converts thinner ice to thicker ice and is applied after horizontal transport. When the ice is converging, enough ice ridges to ensure that the ice area does not exceed the grid cell area.
First we specify the participation function: the thickness distribution \(a_P(h) = b(h) \, g(h)\) of the ice participating in ridging. (We use “ridging” as shorthand for all forms of mechanical redistribution, including rafting.) The weighting function \(b(h)\) favors ridging of thin ice and closing of open water in preference to ridging of thicker ice. There are two options for the form of \(b(h)\). If krdg_partic = 0 in the namelist, we follow [74] and set
where \(G(h)\) is the fractional area covered by ice thinner than \(h\), and \(G^*\) is an empirical constant. Integrating \(a_P(h)\) between category boundaries \(H_{n-1}\) and \(H_n\), we obtain the mean value of \(a_P\) in category \(n\):
where \(a_{Pn}\) is the ratio of the ice area ridging (or open water area closing) in category \(n\) to the total area ridging and closing, and \(G_n\) is the total fractional ice area in categories 0 to \(n\). Equation (90) applies to categories with \(G_n < G^*\). If \(G_{n-1} < G^* < G_n\), then Equation (90) is valid with \(G^*\) replacing \(G_n\), and if \(G_{n-1} > G^*\), then \(a_{Pn} = 0\). If the open water fraction \(a_0 > G^*\), no ice can ridge, because “ridging” simply reduces the area of open water. As in [74] we set \(G^* = 0.15\).
If the spatial resolution is too fine for a given time step \(\Delta t\), the weighting function Equation (89) can promote numerical instability. For \(\Delta t = \mbox{1 hour}\), resolutions finer than \(\Delta x \sim \mbox{10 km}\) are typically unstable. The instability results from feedback between the ridging scheme and the dynamics via the ice strength. If the strength changes significantly on time scales less than \(\Delta t\), the viscous-plastic solution of the momentum equation is inaccurate and sometimes oscillatory. As a result, the fields of ice area, thickness, velocity, strength, divergence, and shear can become noisy and unphysical.
A more stable weighting function was suggested by [47]:
When integrated between category boundaries, Equation (91) implies
This weighting function is used if krdg_partic = 1 in the namelist. From Equation (91), the mean value of \(G\) for ice participating in ridging is \(a^*\), as compared to \(G^*/3\) for Equation (89). For typical ice thickness distributions, setting \(a^* = 0.05\) with krdg_partic = 1 gives participation fractions similar to those given by \(G^* = 0.15\) with krdg_partic = 0. See [47] for a detailed comparison of these two participation functions.
Thin ice is converted to thick, ridged ice in a way that reduces the total ice area while conserving ice volume and internal energy. There are two namelist options for redistributing ice among thickness categories. If krdg_redist = 0, ridging ice of thickness \(h_n\) forms ridges whose area is distributed uniformly between \(H_{\min} = 2 h_n\) and \(H_{\max} = 2 \sqrt{H^* h_n}\), as in [27]. The default value of \(H^*\) is 25 m, as in earlier versions of CICE. Observations suggest that \(H^* = 50\) m gives a better fit to first-year ridges [1], although the lower value may be appropriate for multiyear ridges [19]. The ratio of the mean ridge thickness to the thickness of ridging ice is \(k_n = (H_{\min} + H_{\max}) / (2 h_n)\). If the area of category \(n\) is reduced by ridging at the rate \(r_n\), the area of thicker categories grows simultaneously at the rate \(r_n/k_n\). Thus the net rate of area loss due to ridging of ice in category \(n\) is \(r_n(1-1/k_n)\).
The ridged ice area and volume are apportioned among categories in the thickness range \((H_{\min}, H_{\max})\). The fraction of the new ridge area in category \(m\) is
where \(H_L = \max(H_{m-1},H_{\min})\) and \(H_R= \min(H_m,H_{\max})\). The fraction of the ridge volume going to category \(m\) is
This uniform redistribution function tends to produce too little ice in the 3–5 m range and too much ice thicker than 10 m [1]. Observations show that the ITD of ridges is better approximated by a negative exponential. Setting krdg_redist = 1 gives ridges with an exponential ITD [47]:
for \(h \ge H_{\min}\), with \(g_R(h) = 0\) for \(h < H_{\min}\). Here, \(\lambda\) is an empirical e-folding scale and \(H_{\min}=2h_n\) (where \(h_n\) is the thickness of ridging ice). We assume that \(\lambda = \mu h_n^{1/2}\), where \(\mu\) (mu_rdg) is a tunable parameter with units . Thus the mean ridge thickness increases in proportion to \(h_n^{1/2}\), as in [27]. The value \(\mu = 4.0\) gives \(\lambda\) in the range 1–4 m for most ridged ice. Ice strengths with \(\mu = 4.0\) and krdg_redist = 1 are roughly comparable to the strengths with \(H^* = 50\) m and krdg_redist = 0.
From Equation (95) it can be shown that the fractional area going to category \(m\) as a result of ridging is
The fractional volume going to category \(m\) is
Equations (96) and (97) replace Equations (93) and (94) when krdg_redist = 1.
Internal ice energy is transferred between categories in proportion to ice volume. Snow volume and internal energy are transferred in the same way, except that a fraction of the snow may be deposited in the ocean instead of added to the new ridge.
The net area removed by ridging and closing is a function of the strain rates. Let \(R_{\mathrm{net}}\) be the net rate of area loss for the ice pack (i.e., the rate of open water area closing, plus the net rate of ice area loss due to ridging). Following [19], \(R_{\mathrm{net}}\) is given by
where \(C_s\) is the fraction of shear dissipation energy that contributes to ridge-building, \(D_D\) is the divergence, and \(\Delta\) is a function of the divergence and shear. These strain rates are computed by the dynamics scheme. The default value of \(C_s\) is 0.25.
Next, define \(R_{\mathrm{tot}} = \sum_{n=0}^N r_n\). This rate is related to \(R_{\mathrm{net}}\) by
Given \(R_{\mathrm{net}}\) from Equation (98), we use Equation (99) to compute \(R_{\mathrm{tot}}\). Then the area ridged in category \(n\) is given by \(a_{rn} = r_n \Delta t\), where \(r_n = a_{Pn} R_{\mathrm{tot}}\). The area of new ridges is \(a_{rn} / k_n\), and the volume of new ridges is \(a_{rn} h_n\) (since volume is conserved during ridging). We remove the ridging ice from category \(n\) and use Equations (93) and (94): (or (96) and (97)) to redistribute the ice among thicker categories.
Occasionally the ridging rate in thickness category \(n\) may be large enough to ridge the entire area \(a_n\) during a time interval less than \(\Delta t\). In this case \(R_{\mathrm{tot}}\) is reduced to the value that exactly ridges an area \(a_n\) during \(\Delta t\). After each ridging iteration, the total fractional ice area \(a_i\) is computed. If \(a_i > 1\), the ridging is repeated with a value of \(R_{\mathrm{net}}\) sufficient to yield \(a_i = 1\).
Two tracers for tracking the ridged ice area and volume are available. The actual tracers are for level (undeformed) ice area (alvl) and volume (vlvl), which are easier to implement for a couple of reasons: (1) ice ridged in a given thickness category is spread out among the rest of the categories, making it more difficult (and expensive) to track than the level ice remaining behind in the original category; (2) previously ridged ice may ridge again, so that simply adding a volume of freshly ridged ice to the volume of previously ridged ice in a grid cell may be inappropriate. Although the code currently only tracks level ice internally, both level ice and ridged ice are offered as history output. They are simply related:
Level ice area fraction and volume increase with new ice formation and decrease steadily via ridging processes. Without the formation of new ice, level ice asymptotes to zero because we assume that both level ice and ridged ice ridge, in proportion to their fractional areas in a grid cell (in the spirit of the ridging calculation itself which does not prefer level ice over previously ridged ice).
The ice strength \(P\) may be computed in either of two ways. If the namelist parameter kstrength = 0, we use the strength formula from [26]:
where \(P^* = 27,500 \, \mathrm {N/m}\) and \(C = 20\) are empirical constants, and \(h\) is the mean ice thickness. Alternatively, setting kstrength = 1 gives an ice strength closely related to the ridging scheme. Following [61], the strength is assumed proportional to the change in ice potential energy \(\Delta E_P\) per unit area of compressive deformation. Given uniform ridge ITDs (krdg_redist = 0), we have
where \(C_P = (g/2)(\rho_i/\rho_w)(\rho_w-\rho_i)\), \(\beta =R_{\mathrm{tot}}/R_{\mathrm{net}} > 1\) from Equation (99), and \(C_f\) is an empirical parameter that accounts for frictional energy dissipation. Following [19], we set \(C_f = 17\). The first term in the summation is the potential energy of ridging ice, and the second, larger term is the potential energy of the resulting ridges. The factor of \(\beta\) is included because \(a_{Pn}\) is normalized with respect to the total area of ice ridging, not the net area removed. Recall that more than one unit area of ice must be ridged to reduce the net ice area by one unit. For exponential ridge ITDs (krdg_redist = 1), the ridge potential energy is modified:
The energy-based ice strength given by Equations (101) or (102) is more physically realistic than the strength given by Equation (100). However, use of Equation (100) is less likely to allow numerical instability at a given resolution and time step. See [47] for more details.
2.2.4.3. Thermodynamics¶
The current CICE version includes three thermodynamics options, the “zero-layer” thermodynamics of [64] (ktherm = 0), the Bitz and Lipscomb model [7] (ktherm = 1) that assumes a fixed salinity profile, and a new “mushy” formulation (ktherm = 2) in which salinity evolves [78]. For each thickness category, CICE computes changes in the ice and snow thickness and vertical temperature profile resulting from radiative, turbulent, and conductive heat fluxes. The ice has a temperature-dependent specific heam to simulate the effect of brine pocket melting and freezing, for ktherm = 1 and 2.
Each thickness category \(n\) in each grid cell is treated as a horizontally uniform column with ice thickness \(h_{in} = v_{in}/a_{in}\) and snow thickness \(h_{sn} = v_{sn}/a_{in}\). (Henceforth we omit the category index \(n\).) Each column is divided into \(N_i\) ice layers of thickness \(\Delta h_i = h_i/N_i\) and \(N_s\) snow layers of thickness \(\Delta h_s = h_s/N_s\). The surface temperature (i.e., the temperature of ice or snow at the interface with the atmosphere) is \(T_{sf}\), which cannot exceed . The temperature at the midpoint of the snow layer is \(T_s\), and the midpoint ice layer temperatures are \(T_{ik}\), where \(k\) ranges from 1 to \(N_i\). The temperature at the bottom of the ice is held at \(T_f\), the freezing temperature of the ocean mixed layer. All temperatures are in degrees Celsius unless stated otherwise.
Each ice layer has an enthalpy \(q_{ik}\), defined as the negative of the energy required to melt a unit volume of ice and raise its temperature to . Because of internal melting and freezing in brine pockets, the ice enthalpy depends on the brine pocket volume and is a function of temperature and salinity. We can also define a snow enthalpy \(q_s\), which depends on temperature alone.
Given surface forcing at the atmosphere–ice and ice–ocean interfaces along with the ice and snow thicknesses and temperatures/enthalpies at time \(m\), the thermodynamic model advances these quantities to time \(m+1\) (ktherm=2 also advances salinity). The calculation proceeds in two steps. First we solve a set of equations for the new temperatures, as discussed in Section New temperatures. Then we compute the melting, if any, of ice or snow at the top surface, and the growth or melting of ice at the bottom surface, as described in Section Growth and melting. We begin by describing the surface forcing parameterizations, which are closely related to the ice and snow surface temperatures.
2.2.4.3.1. Melt ponds¶
Three explicit melt pond parameterizations are available in CICE, and all must use the delta-Eddington radiation scheme, described below. The default (ccsm3) shortwave parameterization incorporates melt ponds implicitly by adjusting the albedo based on surface conditions.
For each of the three explicit parameterizations, a volume \(\Delta V_{melt}\) of melt water produced on a given category may be added to the melt pond liquid volume:
where
is the fraction of the total melt water available that is added to the ponds, \(\rho_i\) and \(\rho_s\) are ice and snow densities, \(\Delta h_i\) and \(\Delta h_s\) are the thicknesses of ice and snow that melted, and \(F_{rain}\) is the rainfall rate. Namelist parameters are set for the level-ice (tr_pond_lvl) parameterization; in the cesm and topo pond schemes the standard values of \(r_{max}\) and \(r_{min}\) are 0.7 and 0.15, respectively.
Radiatively, the surface of an ice category is divided into fractions of snow, pond and bare ice. In these melt pond schemes, the actual pond area and depth are maintained throughout the simulation according to the physical processes acting on it. However, snow on the sea ice and pond ice may shield the pond and ice below from solar radiation. These processes do not alter the actual pond volume; instead they are used to define an “effective pond fraction” (and likewise, effective pond depth, snow fraction and snow depth) used only for the shortwave radiation calculation.
In addition to the physical processes discussed below, tracer equations and definitions for melt ponds are also described in Section Tracers and Figure 1.
CESM formulation (tr_pond_cesm = true)
Melt pond area and thickness tracers are carried on each ice thickness category as in Section Tracers. Defined simply as the product of pond area, \(a_p\), and depth, \(h_p\), the melt pond volume, \(V_{p}\), grows through addition of ice or snow melt water or rain water, and shrinks when the ice surface temperature becomes cold,
where \(dh_{i}\) and \(dh_{s}\) represent ice and snow melt at the top surface of each thickness category and \(r_2=0.01\). Here, \(T_p\) is a reference temperature equal to -2 \(^\circ\)C. Pond depth is assumed to be a linear function of the pond fraction (\(h_p=\delta_p a_p\)) and is limited by the category ice thickness (\(h_p \le 0.9 h_i\)). The pond shape (pndaspect) \(\delta_p = 0.8\) in the standard CESM pond configuration. The area and thickness are computed according to the assumed pond shape, and the pond area is then reduced in the presence of snow for the radiation calculation. Ponds are allowed only on ice at least 1 cm thick. This formulation differs slightly from that documented in [30].
Topographic formulation (tr_pond_topo = true)
The principle concept of this scheme is that melt water runs downhill under the influence of gravity and collects on sea ice with increasing surface height starting at the lowest height [20][21][22]. Thus, the topography of the ice cover plays a crucial role in determining the melt pond cover. However, CICE does not explicitly represent the topography of sea ice. Therefore, we split the existing ice thickness distribution function into a surface height and basal depth distribution assuming that each sea ice thickness category is in hydrostatic equilibrium at the beginning of the melt season. We then calculate the position of sea level assuming that the ice in the whole grid cell is rigid and in hydrostatic equilibrium.
Figure 6
Figure 6 : (a) Schematic illustration of the relationship between the height of the pond surface \(h_{pnd,tot}\), the volume of water \(V_{Pk}\) required to completely fill up to category \(k\), the volume of water \(V_{P} - V_{Pk}\), and the depth to which this fills up category \(k + 1\). Ice and snow areas \(a_i\) and \(a_s\) are also depicted. The volume calculation takes account of the presence of snow, which may be partially or completely saturated. (b) Schematic illustration indicating pond surface height \(h_{pnd,tot}\) and sea level \(h_{sl}\) measured with respect to the thinnest surface height category \(h_{i1}\), the submerged portion of the floe \(h_{sub}\), and hydraulic head \(\Delta H\) . A positive hydraulic head (pond surface above sea level) will flush melt water through the sea ice into the ocean; a negative hydraulic head can drive percolation of sea water up onto the ice surface. Here, \(\alpha=0.6\) and \(\beta=0.4\) are the surface height and basal depth distribution fractions. The height of the steps is the height of the ice above the reference level, and the width of the steps is the area of ice of that height. The illustration does not imply a particular assumed topography, rather it is assumed that all thickness categories are present at the sub-grid scale so that water will always flow to the lowest surface height class.
Once a volume of water is produced from ice and snow melting, we calculate the number of ice categories covered by water. At each time step, we construct a list of volumes of water \(\{V_{P1}, V_{P2}, . . . V_{P,k-1}, V_{Pk},\) \(V_{P,k+1}, . . . \}\), where \(V_{Pk}\) is the volume of water required to completely cover the ice and snow in the surface height categories from \(i = 1\) up to \(i = k\). The volume \(V_{Pk}\) is defined so that if the volume of water \(V_{P}\) is such that \(V_{Pk} < V_{P} < V_{P,k+1}\) then the snow and ice in categories \(n = 1\) up to \(n = k + 1\) are covered in melt water. Figure 6 (a) depicts the areas covered in melt water and saturated snow on the surface height (rather than thickness) categories \(h_{top,k}\). Note in the CICE code, we assume that \(h_{top,n}/h_{in} = 0.6\) (an arbitrary choice). The fractional area of the \(n\)th category covered in snow is \(a_{sn}\). The volume \(V_{P1}\), which is the region with vertical hatching, is the volume of water required to completely fill up the first thickness category, so that any extra melt water must occupy the second thickness category, and it is given by the expression
where \(V_{sw}\) is the fraction of the snow volume that can be occupied by water, and \(h_{s1}\) is the snow depth on ice height class 1. In a similar way, the volume required to fill up the first and second surface categories, \(V_{P2}\), is given by
The general expression for volume \(V_{Pk}\) is given by
(Note that we have implicitly assumed that \(h_{si} < h_{top,k+1} - h_{top,k}\) for all \(k\).) No melt water can be stored on the thickest ice thickness category. If the melt water volume exceeds the volume calculated above, the remaining melt water is released to the ocean.
At each time step, the pond height above the level of the thinnest surface height class, that is, the maximum pond depth, is diagnosed from the list of volumes \(V_{Pk}\). In particular, if the total volume of melt water \(V_{P}\) is such that \(V_{Pk} < V_{P} < V_{P,k+1}\) then the pond height \(h_{pnd,tot}\) is
where \(h_{par}\) is the height of the pond above the level of the ice in class \(k\) and partially fills the volume between \(V_{P,k}\) and \(V_{P,k+1}\). From Figure 6 (a) we see that \(h_{top,k} - h_{top,1}\) is the height of the melt water, which has volume \(V_{Pk}\), which completely fills the surface categories up to category \(k\). The remaining volume, \(V_{P} - V_{Pk}\), partially fills category \(k + 1\) up to the height \(h_{par}\) and there are two cases to consider: either the snow cover on category \(k + 1\), with height \(h_{s,k+1}\), is completely covered in melt water (i.e., \(h_{par} > h_{s,k+1}\)), or it is not (i.e., \(h_{par} \le h_{s,k+1}\)). From conservation of volume, we see from Figure 6 (a) that for an incompletely to completely saturated snow cover on surface ice class \(k + 1\),
and for a saturated snow cover with water on top of the snow on surface ice class \(k + 1\),
As the melting season progresses, not only does melt water accumulate upon the upper surface of the sea ice, but the sea ice beneath the melt water becomes more porous owing to a reduction in solid fraction [16]. The hydraulic head of melt water on sea ice (i.e., its height above sea level) drives flushing of melt water through the porous sea ice and into the underlying ocean. The mushy thermodynamics scheme (ktherm \(= 2\)) handles flushing. For ktherm \(\ne 2\) we model the vertical flushing rate using Darcy’s law for flow through a porous medium
where \(w\) is the vertical mass flux per unit perpendicular cross-sectional area (i.e., the vertical component of the Darcy velocity), \(\Pi_v\) is the vertical component of the permeability tensor (assumed to be isotropic in the horizontal), \(\mu\) is the viscosity of water, \(\rho_o\) is the ocean density, \(g\) is gravitational acceleration, \(\Delta H\) is the the hydraulic head, and \(h_i\) is the thickness of the ice through which the pond flushes. As proposed by [25] the vertical permeability of sea ice can be calculated from the liquid fraction \(\phi\):
Since the solid fraction varies throughout the depth of the sea ice, so does the permeability. The rate of vertical drainage is determined by the lowest (least permeable) layer, corresponding to the highest solid fraction. From the equations describing sea ice as a mushy layer [18], the solid fraction is determined by:
where \(S\) is the bulk salinity of the ice, \(S_{br}(T)\) is the concentration of salt in the brine at temperature \(T\) and \(c_i\) is the concentration of salt in the ice crystals (set to zero).
The hydraulic head is given by the difference in height between the upper surface of the melt pond \(h_{pnd,tot}\) and the sea level \(h_{sl}\). The value of the sea level \(h_{sl}\) is calculated from
where \(0.4 \sum^{N}_{n=1} a_{in} h_{i,n}\) is the mean thickness of the basal depth classes, and \(h_{sub}\) is the depth of the submerged portion of the floe. Figure 6 (b) depicts the relationship between the hydraulic head and the depths and heights that appear in Equation (112). The depth of the submerged portion of the floe is determined from hydrostatic equilibrium to be
where \(\rho_m\) is the density of melt water, \(V_P\) is the total pond volume, \(V_s\) is the total snow volume, and \(V_i\) is the total ice volume.
When the surface energy balance is negative, a layer of ice is formed at the upper surface of the ponds. The rate of growth of the ice lid is given by the Stefan energy budget at the lid-pond interface
where \(L_0\) is the latent heat of fusion of pure ice per unit volume, \(T_i\) and \(T_p\) are the ice surface and pond temperatures, and \(k_i\) and \(k_p\) are the thermal conductivity of the ice lid and pond respectively. The second term on the right hand-side is close to zero since the pond is almost uniformly at the freezing temperature [73]. Approximating the temperature gradient in the ice lid as linear, the Stefan condition yields the classic Stefan solution for ice lid depth
where \(\Delta T\) is the temperature difference between the top and the bottom of the lid. Depending on the surface flux conditions the ice lid can grow, partially melt, or melt completely. Provided that the ice lid is thinner than a critical lid depth (1 cm is suggested) then the pond is regarded as effective, i.e. the pond affects the optical properties of the ice cover. Effective pond area and pond depth for each thickness category are passed to the radiation scheme for calculating albedo. Note that once the ice lid has exceeded the critical thickness, snow may accumulate on the lid causing a substantial increase in albedo. In the current CICE model, melt ponds only affect the thermodynamics of the ice through the albedo. To conserve energy, the ice lid is dismissed once the pond is completely refrozen.
As the sea ice area shrinks due to melting and ridging, the pond volume over the lost area is released to the ocean immediately. In [21], the pond volume was carried as an ice area tracer, but in [22] and here, pond area and thickness are carried as separate tracers, as in Section Tracers.
Unlike the cesm and level-ice melt pond schemes, the liquid pond water in the topo parameterization is not necessarily virtual; it can be withheld from being passed to the ocean model until the ponds drain by setting the namelist variable l_mpond_fresh = .true. The refrozen pond lids are still virtual. Extra code needed to track and enforce conservation of water has been added to ice_itd.F90 (subroutine zap_small_areas), ice_mechred.F90 (subroutine ridge_shift), ice_therm_itd.F90 (subroutines linear_itd and lateral_melt), and ice_therm_vertical.F90 (subroutine thermo_vertical), along with global diagnostics in ice_diagnostics.F90.
Level-ice formulation (tr_pond_lvl = true)
This meltpond parameterization represents a combination of ideas from the empirical CESM melt pond scheme and the topo approach, and is documented in [36]. The ponds evolve according to physically based process descriptions, assuming a thickness-area ratio for changes in pond volume. A novel aspect of the new scheme is that the ponds are carried as tracers on the level (undeformed) ice area of each thickness category, thus limiting their spatial extent based on the simulated sea ice topography. This limiting is meant to approximate the horizontal drainage of melt water into depressions in ice floes. (The primary difference between the level-ice and topo meltpond parameterizations lies in how sea ice topography is taken into account when determining the areal coverage of ponds.) Infiltration of the snow by melt water postpones the appearance of ponds and the subsequent acceleration of melting through albedo feedback, while snow on top of refrozen pond ice also reduces the ponds’ effect on the radiation budget.
Melt pond processes, described in more detail below, include addition of liquid water from rain, melting snow and melting surface ice, drainage of pond water when its weight pushes the ice surface below sea level or when the ice interior becomes permeable, and refreezing of the pond water. If snow falls after a layer of ice has formed on the ponds, the snow may block sunlight from reaching the ponds below. When melt water forms with snow still on the ice, the water is assumed to infiltrate the snow. If there is enough water to fill the air spaces within the snowpack, then the pond becomes visible above the snow, thus decreasing the albedo and ultimately causing the snow to melt faster. The albedo also decreases as snow depth decreases, and thus a thin layer of snow remaining above a pond-saturated layer of snow will have a lower albedo than if the melt water were not present.
The level-ice formulation assumes a thickness-area ratio for changes in pond volume, while the CESM scheme assumes this ratio for the total pond volume. Pond volume changes are distributed as changes to the area and to the depth of the ponds using an assumed aspect ratio, or shape, given by the parameter \(\delta_p\) (pndaspect), \(\delta_p = {\Delta h_p / \Delta a_{p}}\) and \(\Delta V = \Delta h_p \Delta a_{p} = \delta_p\Delta a_p^2 = \Delta h_{p}^2/\delta_p\). Here, \(a_{p} = a_{pnd} a_{lvl}\), the mean pond area over the ice.
Given the ice velocity \(\bf u\), conservation equations for level ice fraction \(a_{lvl}a_i\), pond area fraction \(a_{pnd}a_{lvl}a_i\), pond volume \(h_{pnd}a_{pnd}a_{lvl}a_i\) and pond ice volume \(h_{ipnd}a_{pnd}a_{lvl}a_i\) are
(We have dropped the category subscript here, for clarity.) Equations (118) and (119) express conservation of melt pond volume and pond ice volume, but in this form highlight that the quantities tracked in the code are the tracers \(h_{pnd}\) and \(h_{ipnd}\), pond depth and pond ice thickness. Likewise, the level ice fraction \(a_{lvl}\) is a tracer on ice area fraction (Equation (116)), and pond fraction \(a_{pnd}\) is a tracer on level ice (Equation (117)).
Pond ice. The ponds are assumed to be well mixed fresh water, and therefore their temperature is 0\(^\circ\)C. If the air temperature is cold enough, a layer of clear ice may form on top of the ponds. There are currently three options in the code for refreezing the pond ice. Only option A tracks the thickness of the lid ice using the tracer \(h_{ipnd}\) and includes the radiative effect of snow on top of the lid.
A. The frzpnd = ‘hlid’ option uses a Stefan approximation for growth of fresh ice and is invoked only when \(\Delta V_{melt}=0\).
The basic thermodynamic equation governing ice growth is
assuming a linear temperature profile through the ice thickness \(h_i\). In discrete form, the solution is
where
When \(\Delta V_{melt}>0\), any existing pond ice may also melt. In this case,
where \(F_\circ\) is the net downward surface flux.
In either case, the change in pond volume associated with growth or melt of pond ice is
where \(\rho_0\) is the density of fresh water.
B. The frzpnd = ‘cesm’ option uses the same empirical function as in the CESM melt pond parameterization.
Radiative effects. Freshwater ice that has formed on top of a melt pond is assumed to be perfectly clear. Snow may accumulate on top of the pond ice, however, shading the pond and ice below. The depth of the snow on the pond ice is initialized as \(h_{ps}^0 = F_{snow}\Delta t\) at the first snowfall after the pond ice forms. From that time until either the pond ice or the pond snow disappears, the pond snow depth tracks the depth of snow on sea ice (\(h_s\)) using a constant difference \(\Delta\). As \(h_s\) melts, \(h_{ps}=h_s-\Delta\) will be reduced to zero eventually, at which time the pond ice is fully uncovered and shortwave radiation passes through.
To prevent a sudden change in the shortwave reaching the sea ice (which can prevent the thermodynamics from converging), thin layers of snow on pond ice are assumed to be patchy, thus allowing the shortwave flux to increase gradually as the layer thins. This is done using the same parameterization for patchy snow as is used elsewhere in CICE, but with its own parameter \(h_{s1}\):
If any of the pond ice melts, the radiative flux allowed to pass through the ice is reduced by the (roughly) equivalent flux required to melt that ice. This is accomplished (approximately) with \(a_{pnd}^{eff} = (1-f_{frac})a_{pnd}a_{lvl}\), where (see Equation (122))
Snow infiltration by pond water. If there is snow on top of the sea ice, melt water may infiltrate the snow. It is a “virtual process” that affects the model’s thermodynamics through the input parameters of the radiation scheme; it does not melt the snow or affect the snow heat content.
A snow pack is considered saturated when its percentage of liquid water content is greater or equal to 15% (Sturm and others, 2009). We assume that if the volume fraction of retained melt water to total liquid content
then effectively there are no meltponds present, that is, \(a_{pnd}^{eff}=h_{pnd}^{eff}=0\). Otherwise, we assume that the snowpack is saturated with liquid water.
We assume that all of the liquid water accumulates at the base of the snow pack and would eventually melt the surrounding snow. Two configurations are therefore possible, (1) the top of the liquid lies below the snow surface and (2) the liquid water volume overtops the snow, and all of the snow is assumed to have melted into the pond. The volume of void space within the snow that can be filled with liquid melt water is
and we compare \(V_p\) with \(V_{mx}\).
Case 1: For \(V_p < V_{mx}\), we define \(V_p^{eff}\) to be the volume of void space filled by the volume \(V_p\) of melt water: \(\rho_0 V_p = (\rho_0-\rho_s) V_p^{eff},\) or in terms of depths,
The liquid water under the snow layer is not visible and therefore the ponds themselves have no direct impact on the radiation (\(a_{pnd}^{eff}=h_{pnd}^{eff}=0\)), but the effective snow thickness used for the radiation scheme is reduced to
Here, the factor \(a_p=a_{pnd}a_{lvl}\) averages the reduced snow depth over the ponds with the full snow depth over the remainder of the ice; that is, \(h_s^{eff} = h_s(1-a_p) + (h_s -h_p^{eff})a_p.\)
Case 2: Similarly, for \(V_p \ge V_{mx}\), the total mass in the liquid is \(\rho_0 V_p + \rho_s V_s = \rho_0 V_p^{eff},\) or
Thus the effective depth of the pond is the depth of the whole slush layer \(h_p^{eff}\). In this case, \(a_{pnd}^{eff}=a_{pnd}a_{lvl}\).
Drainage. A portion \(1-r\) of the available melt water drains immediately into the ocean. Once the volume changes described above have been applied and the resulting pond area and depth calculated, the pond depth may be further reduced if the top surface of the ice would be below sea level or if the sea ice becomes permeable.
We require that the sea ice surface remain at or above sea level. If the weight of the pond water would push the mean ice–snow interface of a thickness category below sea level, some or all of the pond water is removed to bring the interface back to sea level via Archimedes’ Principle written in terms of the draft \(d\),
There is a separate freeboard calculation in the thermodynamics which considers only the ice and snow and converts flooded snow to sea ice. Because the current melt ponds are “virtual” in the sense that they only have a radiative influence, we do not allow the pond mass to change the sea ice and snow masses at this time, although this issue may need to be reconsidered in the future, especially for the Antarctic.
The mushy thermodynamics scheme (ktherm \(= 2\)) handles flushing. For ktherm \(\ne 2\), the permeability of the sea ice is calculated using the internal ice temperatures \(T_i\) (computed from the enthalpies as in the sea ice thermodynamics). The brine salinity and liquid fraction are given by [55] [eq 3.6] \(S_{br} = {1/ (10^{-3} - 0.054/T_i)}\) and \(\phi = S/S_{br}\), where \(S\) is the bulk salinity of the combined ice and brine. The ice is considered permeable if \(\phi \ge 0.05\) with a permeability of \(p=3\times 10^{-8}\min(\phi^3)\) (the minimum being taken over all of the ice layers). A hydraulic pressure head is computed as \(P=g\rho_w\Delta h\) where \(\Delta h\) is the height of the pond and sea ice above sea level. Then the volume of water drained is given by
where \(d_p\) is a scaling factor (dpscale), and \(\mu=1.79\times 10^{-3}\) kg m \(^{-1}\) s \(^{-1}\) is the dynamic viscosity.
Conservation elsewhere. When ice ridges and when new ice forms in open water, the level ice area changes and ponds must be handled appropriately. For example, when sea ice deforms, some of the level ice is transformed into ridged ice. We assume that pond water (and ice) on the portion of level ice that ridges is lost to the ocean. All of the tracer volumes are altered at this point in the code, even though \(h_{pnd}\) and \(h_{ipnd}\) should not change; compensating factors in the tracer volumes cancel out (subroutine ridge_shift in ice_mechred.F90).
When new ice forms in open water, level ice is added to the existing sea ice, but the new level ice does not yet have ponds on top of it. Therefore the fractional coverage of ponds on level ice decreases (thicknesses are unchanged). This is accomplished in ice_therm_itd.F90 (subroutine add_new_ice) by maintaining the same mean pond area in a grid cell after the addition of new ice,
and solving for the new pond area tracer \(a_{pnd}^\prime\) given the newly formed ice area \(\Delta a_i = \Delta a_{lvl}\).
2.2.4.3.2. Thermodynamic surface forcing balance¶
The net surface energy flux from the atmosphere to the ice (with all fluxes defined as positive downward) is
where \(F_s\) is the sensible heat flux, \(F_l\) is the latent heat flux, \(F_{L\downarrow}\) is the incoming longwave flux, \(F_{L\uparrow}\) is the outgoing longwave flux, \(F_{sw}\) is the incoming shortwave flux, \(\alpha\) is the shortwave albedo, and \(i_0\) is the fraction of absorbed shortwave flux that penetrates into the ice. The albedo may be altered by the presence of melt ponds. Each of the explicit melt pond parameterizations (CESM, topo and level-ice ponds) should be used in conjunction with the Delta-Eddington shortwave scheme, described below.
Shortwave radiation: Delta-Eddington
Two methods for computing albedo and shortwave fluxes are available, the default (“ccsm3”) method, described below, and a multiple scattering radiative transfer scheme that uses a Delta-Eddington approach. “Inherent” optical properties (IOPs) for snow and sea ice, such as extinction coefficient and single scattering albedo, are prescribed based on physical measurements; reflected, absorbed and transmitted shortwave radiation (“apparent” optical properties) are then computed for each snow and ice layer in a self-consistent manner. Absorptive effects of inclusions in the ice/snow matrix such as dust and algae can also be included, along with radiative treatment of melt ponds and other changes in physical properties, for example granularization associated with snow aging. The Delta-Eddington formulation is described in detail in [9]. Since publication of this technical paper, a number of improvements have been made to the Delta-Eddington scheme, including a surface scattering layer and internal shortwave absorption for snow, generalization for multiple snow layers and more than four layers of ice, and updated IOP values.
The namelist parameters R_ice and R_pnd adjust the albedo of bare or ponded ice by the product of the namelist value and one standard deviation. For example, if R_ice = 0.1, the albedo increases by \(0.1\sigma\). Similarly, setting R_snw = 0.1 decreases the snow grain radius by \(0.1\sigma\) (thus increasing the albedo). Two additional tuning parameters are available for this scheme, dT_mlt and rsnw_mlt. dT_mlt is the temperature change needed for a change in snow grain radius from non-melting to melting, and rsnw_mlt is the maximum snow grain radius when melting. An absorption coefficient for algae (kalg) may also be set. See [9] for details; the CESM melt pond and Delta-Eddington parameterizations are further explained and validated in [30].
Shortwave radiation: CCSM3
In the parameterization used in the previous version of the Community Climate System Model (CCSM3), the albedo depends on the temperature and thickness of ice and snow and on the spectral distribution of the incoming solar radiation. Albedo parameters have been chosen to fit observations from the SHEBA field experiment. For \(T_{sf} < -1^{circ}C\) and \(h_i > ` `ahmax\), the bare ice albedo is 0.78 for visible wavelengths (\(<700\)nm) and 0.36 for near IR wavelengths (\(>700\)nm). As \(h_i\) decreases from ahmax to zero, the ice albedo decreases smoothly (using an arctangent function) to the ocean albedo, 0.06. The ice albedo in both spectral bands decreases by 0.075 as \(T_{sf}\) rises from \(-1^{circ}C\) to . The albedo of cold snow (\(T_{sf} < -1^{circ}C\)) is 0.98 for visible wavelengths and 0.70 for near IR wavelengths. The visible snow albedo decreases by 0.10 and the near IR albedo by 0.15 as \(T_{sf}\) increases from \(-1^{circ}C\) to . The total albedo is an area-weighted average of the ice and snow albedos, where the fractional snow-covered area is
and \(h_{snowpatch} = 0.02 \ {\mathrm m}\). The envelope of albedo values is shown in Figure 7. This albedo formulation incorporates the effects of melt ponds implicitly; the explicit melt pond parameterization is not used in this case.
Figure 7
Figure 7 : Albedo as a function of ice thickness and temperature, for the two extrema in snow depth, for the default (CCSM3) shortwave option. Maximum snow depth is computed based on Archimedes’ Principle for the given ice thickness. These curves represent the envelope of possible albedo values.
The net absorbed shortwave flux is \(F_{swabs} = \sum (1-\alpha_j) F_{sw\downarrow}\), where the summation is over four radiative categories (direct and diffuse visible, direct and diffuse near infrared). The flux penetrating into the ice is \(I_0 = i_0 \, F_{swabs}\), where \(i_0 = 0.70 \, (1-f_{snow})\) for visible radiation and \(i_0 = 0\) for near IR. Radiation penetrating into the ice is attenuated according to Beer’s Law:
where \(I(z)\) is the shortwave flux that reaches depth \(z\) beneath the surface without being absorbed, and \(\kappa_i\) is the bulk extinction coefficient for solar radiation in ice, set to \(1.4 \ {\mathrm m^{-1}}\) for visible wavelengths [15]. A fraction \(\exp(-\kappa_i h_i)\) of the penetrating solar radiation passes through the ice to the ocean (\(F_{sw\Downarrow}\)).
Longwave radiation, turbulent fluxes
While incoming shortwave and longwave radiation are obtained from the atmosphere, outgoing longwave radiation and the turbulent heat fluxes are derived quantities. Outgoing longwave takes the standard blackbody form, \(F_{L\uparrow}=\epsilon\sigma \left(T_{sf}^{K}\right)^4\), where \(\epsilon=0.95\) is the emissivity of snow or ice, \(\sigma\) is the Stefan-Boltzmann constant and \(T_{sf}^{K}\) is the surface temperature in Kelvin. (The longwave fluxes are partitioned such that \(\epsilon F_{L\downarrow}\) is absorbed at the surface, the remaining \(\left(1-\epsilon\right)F_{L\downarrow}\) being returned to the atmosphere via \(F_{L\uparrow}\).) The sensible heat flux is proportional to the difference between air potential temperature and the surface temperature of the snow or snow-free ice,
\(C_s\) and \(C_l\) (below) are nonlinear turbulent heat transfer coefficients described in Section Atmosphere. Similarly, the latent heat flux is proportional to the difference between \(Q_a\) and the surface saturation specific humidity \(Q_{sf}\):
where \(q_1 = 1.16378 \times 10^7 \, \mathrm{kg/m^3}\), \(q_2 = 5897.8 \, \mathrm{K}\), \(T_{sf}^K\) is the surface temperature in Kelvin, and \(\rho_a\) is the surface air density.
The net downward heat flux from the ice to the ocean is given by [51]:
where \(\rho_w\) is the density of seawater, \(c_w\) is the specific heat of seawater, \(c_h = 0.006\) is a heat transfer coefficient, \(u_*=\sqrt{\left|\vec{\tau}_w\right|/\rho_w}\) is the friction velocity, and \(T_w\) is the sea surface temperature. A minimum value of \(u_*\) is available; we recommend \(u_{*\min} = 5\times 10^{-4}\) m/s, but the optimal value may depend on the ocean forcing used and can be as low as 0.
\(F_{bot}\) is limited by the total amount of heat available from the ocean, \(F_{frzmlt}\). Additional heat, \(F_{side}\), is used to melt the ice laterally following [52] and [68]. \(F_{bot}\) and the fraction of ice melting laterally are scaled so that \(F_{bot} + F_{side} \ge F_{frzmlt}\) in the case that :math:` F_{frzmlt}<0` (melting; see Section Growth and melting).
2.2.4.3.3. New temperatures¶
Zero-layer thermodynamics (ktherm = 0) An option for zero-layer thermodynamics [64] is available in this version of CICE by setting the namelist parameter ktherm to 0 and changing the number of ice layers, nilyr, in ice_domain_size.F90 to 1. In the zero-layer case, the ice is fresh and the thermodynamic calculations are much simpler than in the other configurations, which we describe here.
Bitz and Lipscomb thermodynamics (ktherm = 1)
The “BL99” thermodynamic sea ice model is based on [53] and [7], and is described more fully in [44]. The vertical salinity profile is prescribed and is unchanging in time. The snow is assumed to be fresh, and the midpoint salinity \(S_{ik}\) in each ice layer is given by
where \(z \equiv (k-1/2)/N_i\), \(S_{\max} = 3.2\) ppt, and \(a=0.407\) and \(b=0.573\) are determined from a least-squares fit to the salinity profile observed in multiyear sea ice by [63]. This profile varies from \(S=0\) at the top surface (\(z = 0\)) to \(S=S_{\max}\) at the bottom surface (\(z=1\)) and is similar to that used by [53]. Equation (125) is fairly accurate for ice that has drained at the top surface due to summer melting. It is not a good approximation for cold first-year ice, which has a more vertically uniform salinity because it has not yet drained. However, the effects of salinity on heat capacity are small for temperatures well below freezing, so the salinity error does not lead to significant temperature errors.
Temperature updates.
Given the temperatures \(T_{sf}^m\), \(T_s^m\), and \(T_{ik}^m\) at time \(m\), we solve a set of finite-difference equations to obtain the new temperatures at time \(m+1\). Each temperature is coupled to the temperatures of the layers immediately above and below by heat conduction terms that are treated implicitly. For example, the rate of change of \(T_{ik}\) depends on the new temperatures in layers \(k-1\), \(k\), and \(k+1\). Thus we have a set of equations of the form
where \({\bf A}\) is a tridiagonal matrix, \({\bf x}\) is a column vector whose components are the unknown new temperatures, and \({\bf b}\) is another column vector. Given \({\bf A}\) and \({\bf b}\), we can compute \({\bf x}\) with a standard tridiagonal solver.
There are four general cases: (1) \(T_{sf} < 0^{circ}C\), snow present; (2) \(T_{sf} = 0^{circ}C\), snow present; (3) \(T_{sf} < 0^{circ}C\), snow absent; and (4) \(T_{sf} = 0^{circ}C\), snow absent. For case 1 we have one equation (the top row of the matrix) for the new surface temperature, \(N_s\) equations for the new snow temperatures, and \(N_i\) equations for the new ice temperatures. For cases 2 and 4 we omit the equation for the surface temperature, which is held at , and for cases 3 and 4 we omit the snow temperature equations. Snow is considered absent if the snow depth is less than a user-specified minimum value, hs_min. (Very thin snow layers are still transported conservatively by the transport modules; they are simply ignored by the thermodynamics.)
The rate of temperature change in the ice interior is given by [53]:
where \(\rho_i = 917 \ \mathrm {kg/m^{3}}\) is the sea ice density (assumed to be uniform), \(c_i(T,S)\) is the specific heat of sea ice, \(K_i(T,S)\) is the thermal conductivity of sea ice, \(I_{pen}\) is the flux of penetrating solar radiation at depth \(z\), and \(z\) is the vertical coordinate, defined to be positive downward with \(z = 0\) at the top surface. If shortwave = ‘default’, the penetrating radiation is given by Beer’s Law:
where \(I_0\) is the penetrating solar flux at the top ice surface and \(\kappa_i\) is an extinction coefficient. If shortwave = ‘dEdd’, then solar absorption is computed by the Delta-Eddington scheme.
The specific heat of sea ice is given to an excellent approximation by [56]
where \(c_0 = 2106\) J/kg/deg is the specific heat of fresh ice at , \(L_0 = 3.34 \times 10^5\) J/kg is the latent heat of fusion of fresh ice at , and \(\mu = 0.054\) deg/ppt is the (liquidus) ratio between the freezing temperature and salinity of brine.
Following [79] and [53], the standard thermal conductivity (conduct = ‘MU71’) is given by
where \(K_0 = 2.03\) W/m/deg is the conductivity of fresh ice and \(\beta = 0.13\) W/m/ppt is an empirical constant. Experimental results [75] suggest that Equation (129) may not be a good description of the thermal conductivity of sea ice. In particular, the measured conductivity does not markedly decrease as \(T\) approaches , but does decrease near the top surface (regardless of temperature).
An alternative parameterization based on the “bubbly brine” model of [58] for conductivity is available (conduct = ‘bubbly’):
where \(\rho_i\) and \(\rho_0=917\) kg/m \(^3\) are densities of sea ice and pure ice. Whereas the parameterization in Equation (129) asymptotes to a constant conductivity of 2.03 W m\(^{-1}\) K \(^{-1}\) with decreasing \(T\), \(K_i\) in Equation (130) continues to increase with colder temperatures.
The equation for temperature changes in snow is analogous to Equation (127), with \(\rho_s = 330\) kg/m \(^3\), \(c_s = c_0\), and \(K_s = 0.30\) W/m/deg replacing the corresponding ice values. If shortwave = ‘default’, then the penetrating solar radiation is equal to zero for snow-covered ice, since most of the incoming sunlight is absorbed near the top surface. If shortwave = ‘dEdd’, however, then \(I_{pen}\) is nonzero in snow layers.
It is possible that more shortwave penetrates into an ice layer than is needed to completely melt the layer, or else it causes the computed temperature to be greater than the melting temperature, which until now has caused the vertical thermodynamics code to abort. A parameter frac = 0.9 sets the fraction of the ice layer than can be melted through. A minimum temperature difference for absorption of radiation is also set, currently dTemp = 0.02 (K). The limiting occurs in ice_therm_vertical.F90, for both the default and delta Eddington radiation schemes. If the available energy would melt through a layer, then penetrating shortwave is first reduced, possibly to zero, and if that is insufficient then the local conductivity is also reduced to bring the layer temperature just to the melting point.
We now convert Equation (127) to finite-difference form. The resulting equations are second-order accurate in space, except possibly at material boundaries, and first-order accurate in time. Before writing the equations in full we give finite-difference expressions for some of the terms.
First consider the terms on the left-hand side of Equation (127). We write the time derivatives as
where \(T\) is the temperature of either ice or snow and \(m\) is a time index. The specific heat of ice layer \(k\) is approximated as
which ensures that energy is conserved during a change in temperature. This can be shown by using Equation (128) to integrate \(c_i \, dT\) from \(T_{ik}^m\) to \(T_{ik}^{m+1}\); the result is \(c_{ik}(T_{ik}^{m+1} - T_{ik}^m)\), where \(c_{ik}\) is given by Equation (131). The specific heat is a nonlinear function of \(T_{ik}^{m+1}\), the unknown new temperature. We can retain a set of linear equations, however, by initially guessing \(T_{ik}^{m+1} = T_{ik}^m\) and then iterating the solution, updating \(T_{ik}^{m+1}\) in Equation (131) with each iteration until the solution converges.
Next consider the first term on the right-hand side of Equation (127). The first term describes heat diffusion and is discretized for a given ice or snow layer \(k\) as
where \(\Delta h\) is the layer thickness and \(K_{k}\) is the effective conductivity at the upper boundary of layer \(k\). This discretization is centered and second-order accurate in space, except at the boundaries. The flux terms on the right-hand side (RHS) are treated implicitly; i.e., they depend on the temperatures at the new time \(m+1\). The resulting scheme is first-order accurate in time and unconditionally stable. The effective conductivity \(K^*\) at the interface of layers \(k-1\) and \(k\) is defined as
which reduces to the appropriate values in the limits \(K_k \gg K_{k-1}\) (or vice versa) and \(h_k \gg h_{k-1}\) (or vice versa). The effective conductivity at the top (bottom) interface of the ice-snow column is given by \(K^*=2K/\Delta h\), where \(K\) and \(\Delta h\) are the thermal conductivity and thickness of the top (bottom) layer. The second term on the RHS of Equation (127) is discretized as
where \(\tau_k\) is the fraction of the penetrating solar radiation \(I_0\) that is transmitted through layer \(k\) without being absorbed.
We now construct a system of equations for the new temperatures. For \(T_{sf}<0^{circ}C\) we require
where \(F_{ct}\) is the conductive flux from the top surface to the ice interior, and both fluxes are evaluated at time \(m+1\). Although \(F_0\) is a nonlinear function of \(T_{sf}\), we can make the linear approximation
where \(T_{sf}^*\) is the surface temperature from the most recent iteration, and \(F_0^*\) and \((dF_0/dT_{sf})^*\) are functions of \(T_{sf}^*\). We initialize \(T_{sf}^* = T_{sf}^m\) and update it with each iteration. Thus we can rewrite Equation (133) as
Rearranging terms, we obtain
the first equation in the set of equations (126). The temperature change in ice/snow layer \(k\) is
where \(T_0 = T_{sf}\) in the equation for layer 1. In tridiagonal matrix form, Equation (135) becomes
where \(\eta_k = \Delta t/(\rho_k c_k \Delta h_k)\). In the equation for the bottom ice layer, the temperature at the ice–ocean interface is held fixed at \(T_f\), the freezing temperature of the mixed layer; thus the last term on the LHS is known and is moved to the RHS. If \(T_{sf} = 0^{circ}C\) , then there is no surface flux equation. In this case the first equation in Equation (126) is similar to Equation (136), but with the first term on the LHS moved to the RHS.
These equations are modified if \(T_{sf}\) and \(F_{ct}\) are computed within the atmospheric model and passed to CICE (calc_Tsfc = false; see Section Atmosphere). In this case there is no surface flux equation. The top layer temperature is computed by an equation similar to Equation (136) but with the first term on the LHS replaced by \(\eta_1 F_{ct}\) and moved to the RHS. The main drawback of treating the surface temperature and fluxes explicitly is that the solution scheme is no longer unconditionally stable. Instead, the effective conductivity in the top layer must satisfy a diffusive CFL condition:
For thin layers and typical coupling intervals (\(\sim 1\) hr), \(K^*\) may need to be limited before being passed to the atmosphere via the coupler. Otherwise, the fluxes that are returned to CICE may result in oscillating, highly inaccurate temperatures. The effect of limiting is to treat the ice as a poor heat conductor. As a result, winter growth rates are reduced, and the ice is likely to be too thin (other things being equal). The values of hs_min and \(\Delta t\) must therefore be chosen with care. If hs_min is too small, frequent limiting is required, but if hs_min is too large, snow will be ignored when its thermodynamic effects are significant. Likewise, infrequent coupling requires more limiting, whereas frequent coupling is computationally expensive.
This completes the specification of the matrix equations for the four cases. We compute the new temperatures using a tridiagonal solver. After each iteration we check to see whether the following conditions hold:
- \(T_{sf} \leq 0^{circ}C\).
- The change in \(T_{sf}\) since the previous iteration is less than a prescribed limit, \(\Delta T_{\max}\).
- \(F_0 \geq F_{ct}\). (If \(F_0 < F_{ct}\), ice would be growing at the top surface, which is not allowed.)
- The rate at which energy is added to the ice by the external fluxes equals the rate at which the internal ice energy is changing, to within a prescribed limit \(\Delta F_{\max}\).
We also check the convergence rate of \(T_{sf}\). If \(T_{sf}\) is oscillating and failing to converge, we average temperatures from successive iterations to improve convergence. When all these conditions are satisfied—usually within two to four iterations for \(\Delta T_{\max} \approx 0.01^{circ}C\) and \(\Delta F_{max} \approx 0.01 \ \mathrm{W/m^2}\)—the calculation is complete.
To compute growth and melt rates (Section Growth and melting, we derive expressions for the enthalpy \(q\). The enthalpy of snow (or fresh ice) is given by
Sea ice enthalpy is more complex, because of brine pockets whose salinity varies inversely with temperature. Since the salinity is prescribed, there is a one-to-one relationship between temperature and enthalpy. The specific heat of sea ice, given by Equation (128), includes not only the energy needed to warm or cool ice, but also the energy used to freeze or melt ice adjacent to brine pockets. Equation (128) can be integrated to give the energy \(\delta_e\) required to raise the temperature of a unit mass of sea ice of salinity \(S\) from \(T\) to \(T^\prime\):
If we let \(T^\prime = T_{m} \equiv -\mu S\), the temperature at which the ice is completely melted, we have
Multiplying by \(\rho_i\) to change the units from \(\mathrm {J/kg}\) to \(\mathrm {J/m^{3}}\) and adding a term for the energy needed to raise the meltwater temperature to , we obtain the sea ice enthalpy:
Note that Equation (137) is a quadratic equation in \(T\). Given the layer enthalpies we can compute the temperatures using the quadratic formula:
where
The other root is unphysical.
Mushy thermodynamics (ktherm = 2)
The “mushy” thermodynamics option treats the sea ice as a mushy layer [18] in which the ice is assumed to be composed of microscopic brine inclusions surrounded by a matrix of pure water ice. Both enthalpy and salinity are prognostic variables. The size of the brine inclusions is assumed to be much smaller than the size of the ice layers, allowing a continuum approximation: a bulk sea-ice quantity is taken to be the liquid-fraction-weighted average of that quantity in the ice and in the brine.
Enthalpy and mushy physics.
The mush enthalpy, \(q\), is related to the temperature, \(T\), and the brine volume, \(\phi\), by
where \(q_{br}\) is the brine enthalpy, \(q_i\) is the pure ice enthalpy, \(\rho_i\) and \(c_i\) are density and heat capacity of the ice, \(\rho_{w}\) and \(c_{w}\) are density and heat capacity of the brine and \(L_0\) is the latent heat of melting of pure ice. We assume that the specific heats of the ice and brine are fixed at the values of cp_ice and cp_ocn, respectively. The enthalpy is the energy required to raise the temperature of the sea ice to , including both sensible and latent heat changes. Since the sea ice contains salt, it usually will be fully melted at a temperature below \(0^{circ}C\). Equations (137) and (138) are equivalent except for the density used in the term representing the energy required to bring the melt water temperature to (\(\rho_i\) and \(\rho_w\) in equations (137) and (138), respectively).
The liquid fraction, \(\phi\), of sea ice is given by
where the brine salinity, \(S_{br}\), is given by the liquidus relation using the ice temperature.
Within the parameterizations of brine drainage the brine density is a function of brine salinity [55]:
Outside the parameterizations of brine drainage the densities of brine and ice are fixed at the values of \(\rho_w\) and \(\rho_i\), respectively.
The permeability of ice is computed from the liquid fraction as in [25]:
where \(\phi_\Pi\) is 0.05.
The liquidus relation used in the mushy layer module is based on observations of [5]. A piecewise linear relation can be fitted to observations of Z (the ratio of mass of salt (in g) to mass of pure water (in kg) in brine) to the melting temperature: \(Z = aT + b\). Salinity is the mass of salt (in g) per mass of brine (in kg) so is related to Z by
The data is well fitted with two linear regions,
where
\(T_0\) is the temperature at which the two linear regions meet. Fitting to the data, \(T_0=-7.636^\circ\)C, \(a_1=-18.48 \;\mathrm{g} \;\mathrm{kg}^{-1} \;\mathrm{K}^{-1}\), \(a_2=-10.3085\;\mathrm{g} \;\mathrm{kg}^{-1} \;\mathrm{K}^{-1}\), \(b_1=0\) and \(b_2=62.4 \;\mathrm{g}\;\mathrm{kg}^{-1}\).
Two stage outer iteration.
As for the BL99 thermodynamics [7] there are two qualitatively different situations that must be considered when solving for the vertical thermodynamics: the surface can be melting and at the melting temperature, or the surface can be colder than the melting temperature and not melting. In the BL99 thermodynamics these two situations were treated within the same iterative loop, but here they are dealt with separately. If at the beginning of the time step the ice surface is cold and not melting, we solve the ice temperatures assuming that this is also true at the end of the time step. Once we have solved for the new temperatures we test to see if the answer is consistent with this assumption. If the surface temperature is below the melting temperature then we have found the appropriate consistent solution. If the surface is above the melting temperature at the end of the initial solution attempt, we recalculate the new temperatures assuming the surface temperature is fixed at the melting temperature. Alternatively if the surface is at the melting temperature at the start of a time step, we assume initially that this is also the case at the end of the time step, solve for the new temperatures and then check that the surface conductive heat flux is less than the surface atmospheric heat flux as is required for a melting surface. If this is not the case, the temperatures are recalculated assuming the surface is colder than melting. We have found that solutions of the temperature equations that only treat one of the two qualitatively different solutions at a time are more numerically robust than if both are solved together. The surface state rarely changes qualitatively during the solution so the method is also numerically efficient.
Temperature updates.
During the calculation of the new temperatures and salinities, the liquid fraction is held fixed at the value from the previous time step. Updating the liquid fraction during the Picard iteration described below was found to be numerically unstable. Keeping the liquid fraction fixed drastically improves the numerical stability of the method without significantly changing the solution.
Temperatures are calculated in a similar way to BL99 with an outer Picard iteration of an inner tridiagonal matrix solve. The conservation equation for the internal ice temperatures is
where \(q\) is the sea ice enthalpy, \(K\) is the bulk thermal conductivity of the ice, \(w\) is the vertical Darcy velocity of the brine, \(q_{br}\) is the brine enthalpy and \(F\) is the internally absorbed shortwave radiation. The first term on the right represents heat conduction and the second term represents the vertical advection of heat by gravity drainage and flushing.
The conductivity of the mush is given by
where \(K_i = 2.3\) Wm:math:^{-1}K\(^{-1}\) is the conductivity of pure ice and \(K_{br}=0.5375\) Wm:math:^{-1}K\(^{-1}\) is the conductivity of the brine. The thermal conductivity of brine is a function of temperature and salinity, but here we take it as a constant value for the middle of the temperature range experienced by sea ice, \(-10^\circ\)C [65], assuming the brine liquidus salinity at \(-10^\circ\)C.
We discretize the terms that include temperature in the heat conservation equation as
where the superscript signifies whether the quantity is evaluated at the start (\(t_0\)) or the end (\(t\)) of the time step and the subscript indicates the vertical layer. Writing out the temperature dependence of the enthalpy term we have
The mush thermal conductivities are fixed at the start of the timestep. For the lowest ice layer \(T_{k+1}\) is replaced with \(T_{bot}\), the temperature of the ice base. \(\Delta h\) is the layer thickness and \(z^\prime_k\) is the distance between the \(k-1\) and \(k\) layer centers.
Similarly, for the snow layer temperatures we have the following discretized equation:
For the upper-most layer (either ice layer or snow layer if it present) \(T_{k-1}\) is replaced with \(T_{sf}\), the temperature of the surface.
If the surface is colder than the melting temperature then we also have to solve for the surface temperature, \(T_{sf}\). Here we follow the methodology of BL99 described above.
These discretized temperature equations form a tridiagional matrix for the new temperatures and are solved with a standard tridiagonal solver. A Picard iteration is used to incorporate nonlinearity in the equations. The surface heat flux is a function of surface temperature and with each iteration, the surface heat flux is calculated with the new surface temperature until convergence is achieved. Convergence normally occurs after a few iterations once the temperature changes during an iteration fall below \(5\times10^{-4}\;^\circ\mathrm{C}\) and the energy conservation error falls below 0.9ferrmax.
Salinity updates.
Several physical processes alter the sea ice bulk salinity. New ice forms with the salinity of the sea water from which it formed. Gravity drainage reduces the bulk salinity of newly formed sea ice, while flushing of melt water through the ice also alters the salinity profile.
The salinity equation takes the form
where \(w\) is a vertical Darcy velocity and \(G\) is a source term. The right-hand side depends indirectly on the bulk salinity through the liquid fraction (\(S = \phi S_{br}\)). Since \(\phi\) is fixed for the time step, we solve the salinity equation explicitly after the temperature equation is solved.
A. Gravity drainage. Sea ice initially retains all the salt present in the sea water from which it formed. Cold temperatures near the top surface of forming sea ice result in higher brine salinities there, because the brine is always at its melting temperature. This colder, saltier brine is denser than the underlying sea water and the brine undergoes convective overturning with the ocean. As the dense, cold brine drains out of the ice, it is replaced by fresher seawater, lowering the bulk salinity of the ice. Following [78], gravity drainage is assumed to occur as two simultaneously operating modes: a rapid mode operating principally near the ice base and a slow mode occurring everywhere.
Rapid drainage takes the form of a vertically varying upward Darcy flow. The contribution to the bulk salinity equation for the rapid mode is
where \(S\) is the bulk salinity and \(B_{br}\) is the brine salinity, specified by the liquidus relation with ice temperature. This equation is discretized using an upwind advection scheme,
The upward advective flow also carries heat, contributing a term to the heat conservation Equation (139),
where \(q_{br}\) is the brine enthalpy. This term is discretized as
where the maximum is taken over all the ice layers between layer \(k\) and the ice base. \(\tilde{w}_j\) is given by
where \(Ra_c\) is a critical Rayleigh number and \(Ra(z)\) is the local Rayleigh number at a particular level,
where \(\Delta \rho\) is the difference in density between the brine at \(z\) and the ocean, \(\Pi\) is the minimum permeability between \(z\) and the ocean, \(h\) is the ice thickness, \(\kappa\) is the brine thermal diffusivity and \(\eta\) is the brine dynamic viscosity. Equation ([eq:mushyvel]) reduces the flow rate for Rayleigh numbers below the critical Rayleigh number.
The unmodified flow rate, \(w\), is determined from a hydraulic pressure balance argument for upward flow through the mush and returning downward flow through ice free channels:
where
There are three tunable parameters in the above parameterization, \(a\), the diameter of the channel, \(\Delta x\), the horizontal size of the mush draining through each channel, and \(Ra_c\), the critical Rayleigh number. \(\rho_p\) is the density of brine in the channel which we take to be the density of brine in the mush at the level that the brine is draining from. \(l\) is the thickness of mush from the ice base to the top of the layer in question. We assume that \(\Delta x\) is proportional to \(l\) so that \(\Delta x = 2 \beta l\). \(a\) (a_rapid_mode), \(\beta\) (aspect_rapid_mode) and \(Ra_c\) (Ra_c_rapid_mode) are all namelist parameters with default values of \(0.5\;\mathrm{mm}\), 1 and 10, respectively. The value \(\beta=1\) gives a square aspect ratio for the convective flow in the mush.
The slow drainage mode takes the form of a simple relaxation of bulk salinity:
The decay constant, \(\lambda\), is modeled as
where \(S^\ast\) is a tuning parameter for the drainage strength, \(T_{bot}\) is the basal ice temperature, \(T_{sf}\) is the upper surface temperature and \(h\) is the ice thickness. The bulk salinity relaxes to a value, \(S_c(z)\), given by
where \(S_{br}(z)\) is the brine salinity at depth \(z\) and \(\phi_c\) is a critical liquid fraction. Both \(S^\ast\) and \(\phi_c\) are namelist parameters, dSdt_slow_mode \(=1.5\times10^{-7}\;\mathrm{m}\;\mathrm{s}^{-1}\;\mathrm{K}^{-1}\) and phi_c_slow_mode \(=0.05\).
B. Downwards flushing. Melt pond water drains through sea ice and flushes out brine, reducing the bulk salinity of the sea ice. This is modeled with the mushy physics option as a vertical Darcy flow through the ice that affects both the enthalpy and bulk salinity of the sea ice:
These equations are discretized with an upwind advection scheme. The flushing Darcy flow, \(w_f\), is given by
where \(\overline{\Pi}\) is the harmonic mean of the ice layer permeabilities and \(\Delta h\) is the hydraulic head driving melt water through the sea ice. It is the difference in height between the top of the melt pond and sea level.
Basal boundary condition.
In traditional Stefan problems the ice growth rate is calculated by determining the difference in heat flux on either side of the ice/ocean interface and equating this energy difference to the latent heat of new ice formed. Thus,
where \((1-\phi_i)\) is the solid fraction of new ice formed and the right hand is the difference in heat flux at the ice–ocean interface between the ice side and the ocean side of the interface. However, with mushy layers there is usually no discontinuity in solid fraction across the interface, so \(\phi_i=1\) and Equation (141) cannot be used explicitly. To circumvent this problem we set the interface solid fraction to be 0.15, a value that reproduces observations. \(\phi_i\) is a namelist parameter (phi_i_mushy = 0.85). The basal ice temperature is set to the liquidus temperature \(T_f\) of the ocean surface salinity.
Tracer consistency.
In order to ensure conservation of energy and salt content, the advection routines will occasionally limit changes to either enthalpy or bulk salinity. The mushy thermodynamics routine determines temperature from both enthalpy and bulk salinity. Since the limiting changes performed in the advection routine are not applied consistently (from a mushy physics point of view) to both enthalpy and bulk salinity, the resulting temperature may be changed to be greater than the limit allowed in the thermodynamics routines. If this situation is detected, the code corrects the enthalpy so the temperature is below the limiting value. Conservation of energy is ensured by placing the excess energy in the ocean, and the code writes a warning that this has occurred to the diagnostics file. This situation only occurs with the mushy thermodynamics, and it should only occur very infrequently and have a minimal effect on results. The addition of the heat to the ocean may reduce ice formation by a small amount afterwards.
2.2.4.3.4. Growth and melting¶
Melting at the top surface is given by
where \(q\) is the enthalpy of the surface ice or snow layer [1] (recall that \(q < 0\)) and \(\delta h\) is the change in thickness. If the layer melts completely, the remaining flux is used to melt the layers beneath. Any energy left over when the ice and snow are gone is added to the ocean mixed layer. Ice cannot grow at the top surface due to conductive fluxes; however, snow–ice can form. New snowfall is added at the end of the thermodynamic time step.
Growth and melting at the bottom ice surface are governed by
where \(F_{bot}\) is given by Equation (124) and \(F_{cb}\) is the conductive heat flux at the bottom surface:
If ice is melting at the bottom surface, \(q\) in Equation (143) is the enthalpy of the bottom ice layer. If ice is growing, \(q\) is the enthalpy of new ice with temperature \(T_f\) and salinity \(S_{max}\) (ktherm = 1) or ocean surface salinity (ktherm = 2). This ice is added to the bottom layer.
In general, frazil ice formed in the ocean is added to the thinnest ice category. The new ice is grown in the open water area of the grid cell to a specified minimum thickness; if the open water area is nearly zero or if there is more new ice than will fit into the thinnest ice category, then the new ice is spread over the entire cell.
If the latent heat flux is negative (i.e., latent heat is transferred from the ice to the atmosphere), snow or snow-free ice sublimates at the top surface. If the latent heat flux is positive, vapor from the atmosphere is deposited at the surface as snow or ice. The thickness change of the surface layer is given by
where \(\rho\) is the density of the surface material (snow or ice), and \(L_v = 2.501 \times 10^6 \ \mathrm{J/kg}\) is the latent heat of vaporization of liquid water at . Note that \(\rho L_v\) is nearly an order of magnitude larger than typical values of \(q\). For positive latent heat fluxes, the deposited snow or ice is assumed to have the same enthalpy as the existing surface layer.
After growth and melting, the various ice layers no longer have equal thicknesses. We therefore adjust the layer interfaces, conserving energy, so as to restore layers of equal thickness \(\Delta h_i = h_i / N_i\). This is done by computing the overlap \(\eta_{km}\) of each new layer \(k\) with each old layer \(m\):
where \(z_m\) and \(z_k\) are the vertical coordinates of the old and new layers, respectively. The enthalpies of the new layers are
Lateral melting is accomplished by multiplying the state variables by \(1-r_{side}\), where \(r_{side}\) is the fraction of ice melted laterally [52][68], and adjusting the ice energy and fluxes as appropriate. We assume a floe diameter of 300 m.
Snow ice formation.
At the end of the time step we check whether the snow is deep enough to lie partially below the surface of the ocean (freeboard). From Archimedes’ principle, the base of the snow is at sea level when
Thus the snow base lies below sea level when
In this case, for ktherm = 1 (BL99) we raise the snow base to sea level by converting some snow to ice:
In rare cases this process can increase the ice thickness substantially. For this reason snow–ice conversions are postponed until after the remapping in thickness space (Section Transport in thickness space), which assumes that ice growth during a single time step is fairly small.
For ktherm = 2 (mushy), we model the snow–ice formation process as follows: If the ice surface is below sea level then we replace some snow with the same thickness of sea ice. The thickness change chosen is that which brings the ice surface to sea level. The new ice has a porosity of the snow, which is calculated as
where \(\rho_s\) is the density of snow and \(\rho_i\) is the density of fresh ice. The salinity of the brine occupying the above porosity within the new ice is taken as the sea surface salinity. Once the new ice is formed, the vertical ice and snow layers are regridded into equal thicknesses while conserving energy and salt.
| [1] | The mushy thermodynamics option does not include the enthalpy associated with raising the meltwater temperature to in these calculations, unlike BL99, which does include it. This extra heat is returned to the ocean (or the atmosphere, in the case of evaporation) with the melt water. |