2.2. Fundamental Variables

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.

In addition to an ice thickness distribution, CICE includes an optional capability for a floe size distribution.

Ice floe horizontal size may change through vertical and lateral growth and melting of existing floes, freezing of new ice, wave breaking, and welding of floes in freezing conditions. The floe size distribution (FSD) is a probability function that characterizes this variability. The scheme is based on the theoretical framework described in [15] for a joint floe size and thickness distribution (FSTD), and was implemented by [38] and [Roach19]. The joint floe size distribution is carried as an area-weighted tracer, defined as the fraction of ice belonging to a given thickness category with lateral floe size belong to a given floe size class. This development includes interactions between sea ice and ocean surface waves. Input data on ocean surface wave spectra at a single time is provided for testing, but as with the other CICE datasets, it should not be used for production runs or publications.

Additional information about the ITD and joint FSTD for CICE can be found in the Icepack documentation.

The fundamental equation solved by CICE is [47]:

(1)\[\frac{\partial g}{\partial t} = -\nabla \cdot (g {\bf u}) - \frac{\partial}{\partial h} (f g) + \psi,\]

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.

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 \(0\ ^{\circ}\)C. (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 (1) 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 [7] as adapted for sea ice by [30]; this scheme is discussed in Section Horizontal Transport. Ice is transported in thickness space using the remapping scheme of [29]. The mechanical redistribution scheme, based on [47], [42], [12], [8], and [31] is outlined in the Icepack Documentation. 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 [17], as modified by [5], [16], [18] and [19], or a new elastic-anisotropic-plastic model [53][51][48] to find the velocity, as described in Section Dynamics. Finally, we use a thermodynamic model to compute \(f\). 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.