2.6. Mechanical redistribution

The last term on the right-hand side of Equation (1) is \(\psi\), which describes the redistribution of ice in thickness space due to ridging and other mechanical processes. The mechanical redistribution scheme in Icepack is based on [68], [57], [22], [15], and [39]. 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 [68] and set

(1)\[\begin{split}b(h) = \left\{\begin{array}{ll} \frac{2}{G^*}(1-\frac{G(h)}{G^*}) & \mbox{if $G(h)<G^*$} \\ 0 & \mbox{otherwise} \end{array} \right.\end{split}\]

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\):

(2)\[a_{Pn} = \frac{2}{G^*} (G_n - G_{n-1}) \left( 1 - \frac{G_{n-1}+G_n}{2 G^*} \right),\]

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 (2) applies to categories with \(G_n < G^*\). If \(G_{n-1} < G^* < G_n\), then Equation (2) 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 [68] we set \(G^* = 0.15\).

If the spatial resolution is too fine for a given time step \(\Delta t\), the weighting function Equation (1) 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 [39]:

(3)\[b(h) = \frac{\exp[-G(h)/a^*]} {a^*[1-\exp(-1/a^*)]}\]

When integrated between category boundaries, Equation (3) implies

(4)\[a_{Pn} = \frac {\exp(-G_{n-1}/a^*) - \exp(-G_{n}/a^*)} {1 - \exp(-1/a^*)}\]

This weighting function is used if krdg_partic = 1 in the namelist. From Equation (3), the mean value of \(G\) for ice participating in ridging is \(a^*\), as compared to \(G^*/3\) for Equation (1). 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 [39] 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 [22]. 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 [15]. 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

(5)\[f_m^{\mathrm{area}} = \frac{H_R - H_L} {H_{\max} - H_{\min}},\]

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

(6)\[f_m^{\mathrm{vol}} = \frac{(H_R)^2 - (H_L)^2} {(H_{\max})^2 - (H_{\min})^2}.\]

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 [39]:

(7)\[g_R(h) \propto \exp[-(h - H_{\min})/\lambda]\]

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 [22]. 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 (7) it can be shown that the fractional area going to category \(m\) as a result of ridging is

(8)\[f_m^{\mathrm{area}} = \exp[-(H_{m-1} - H_{\min}) / \lambda] - \exp[-(H_m - H_{\min}) / \lambda].\]

The fractional volume going to category \(m\) is

(9)\[f_m^{\mathrm{vol}} = \frac{(H_{m-1}+\lambda) \exp[-(H_{m-1}-H_{\min})/\lambda] - (H_m + \lambda) \exp[-(H_m - H_{\min}) / \lambda]} {H_{min} + \lambda}.\]

Equations (8) and (9) replace Equations (5) and (6) 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 [15], \(R_{\mathrm{net}}\) is given by

(10)\[R_{\mathrm{net}} = \frac{C_s}{2} (\Delta - |D_D|) - \min(D_D,0),\]

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

(11)\[R_{\mathrm{net}} = \left[ a_{P0} + \sum_{n=1}^N a_{Pn}\left(1-{1\over k_n}\right)\right] R_{\mathrm{tot}}.\]

Given \(R_{\mathrm{net}}\) from Equation (10), we use Equation (11) 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 (5) and (6) (or (8) and (9)) 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 available for output. They are simply related:

(12)\[\begin{split}\begin{aligned} a_{lvl} + a_{rdg} &=& a_i, \\ v_{lvl} + v_{rdg} &=& v_i.\end{aligned}\end{split}\]

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 [21]:

(13)\[P = P^* h \exp[-C(1-a_i)],\]

where \(P^* = 27,500 \, \mathrm {N/m^2}\) 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 [57], 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

(14)\[P = C_f \, C_p \, \beta \sum_{n=1}^{N_C} \left[ -a_{Pn} \, h_n^2 + \frac{a_{Pn}}{k_n} \left( \frac{(H_n^{\max})^3 - (H_n^{\min})^3} {3(H_n^{\max}-H_n^{\min})} \right) \right],\]

where \(C_P = (g/2)(\rho_i/\rho_w)(\rho_w-\rho_i)\), \(\beta =R_{\mathrm{tot}}/R_{\mathrm{net}} > 1\) from Equation (11), and \(C_f\) is an empirical parameter that accounts for frictional energy dissipation. Following [15], 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:

(15)\[P = C_f \, C_p \, \beta \sum_{n=1}^{N_C} \left[ -a_{Pn} \, h_n^2 + \frac{a_{Pn}}{k_n} \left( H_{\min}^2 + 2H_{\min}\lambda + 2 \lambda^2 \right) \right]\]

The energy-based ice strength given by Equations (14) or (15) is more physically realistic than the strength given by Equation (13). However, use of Equation (13) is less likely to allow numerical instability at a given resolution and time step. See [39] for more details.