2.5. 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 (1) by neglecting the first and third terms on the right-hand side. We use the remapping method of [38], 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 (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 (4) and (5) 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 (10) 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 (11) 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 (8) and (6) 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, (8) and (6) apply with \(\eta_R = H_R - H_C\) and \(\eta_n = h_n - H_C\).
Figure Linear approximation of thickness distribution function 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.
Linear approximation of thickness distribution function¶
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 (12) and (13) 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 (8) and (6). 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 [5], 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.