Improving numerical accuracy for the viscous-plastic formulation of sea ice

One goal of this investigation is to address the numerical simulation of the nonlinear VP sea ice model, specifically concerning the accuracy, stability, and efficiency. While many efforts have been made to improve the computational efficiency of sea ice model solvers, e.g. [21, 22, 9, 14], analysis of the corresponding convergence properties is lacking, and indeed many of these methods fail to converge[23]. We note that [23], which proposes and implements an iterated IMplicit-EXplicit (IMEX) time integration to solve the coupled sea ice model monolithically, does provide an analysis of the temporal convergence of the numerical solution. There it is demonstrated that a combination of the second-order Runge-Kutta method for the explicit time integration and a second-order backward difference method for the implicit integration of the momentum equation yields an overall second-order accuracy in time of the numerical solution when compared to a reference solution obtained using a tiny time step ($1$ second). Spatial convergence is investigated in [38], where it is shown that in the VP sea ice model, the simulated velocity field depends on the spatial resolution of the model and approaches the analytical solution as the spatial resolution is increased. However, the study is not quantified in terms of convergence rate. The method in [33] adopts the Crank Nicholson time discretization and centered difference (CD) spatial discretization together with the JFNK solver. Second-order convergence is obtained using a synthetic model by simultaneously refining the spatial resolution and time step. It is not observed everywhere during the simulation test, however.

To the best of our knowledge, convergence with respect to spatial resolution has not been well studied. In particular, no clear conclusion has been drawn in terms of spatial convergence. We adopt this as a starting point in our investigation to explore the related numerical properties. As we focus on spatial convergence, we choose the time step to be sufficiently small so that the time discretization error does not affect the convergence rate. This allows us to test the VP formulation using an explicit time-stepping scheme and therefore avoid the error caused by the nonlinear solver needed for the implicit time-stepping scheme. Based on a constructed analytical solution with appropriately added forcing term to the governing equations, we test for convergence on the VP sea ice model using both a second-order CD spatial discretization scheme as well as a third-order total variation diminishing (TVD) Runge-Kutta time integration scheme.

Nowadays, with increasing mesh resolutions of up to $10$ km [23], the general performance of the existing sea ice solvers is degrading, leading to a significant increase in numerical cost [10]. Higher-order spatial discretization methods may be able to offset this problem. Due to the natural discontinuity feature of ice concentration and ice thickness, traditional higher-order finite difference schemes typically have spurious oscillations near discontinuities (the Gibbs phenomenon), which may pollute smooth regions and even lead to instability, causing blowups of the schemes [34]. The weighted essentially non-oscillatory (WENO) method [27] is designed to achieve higher-order accuracy in smooth regions while sharply resolving discontinuities in an essentially non-oscillatory fashion. This study verifies that these desirable properties hold when implementing WENO for the sea ice model.

In most sea ice models, including the VP model, two idealized thickness levels, namely thick and thin, are often adopted to approximately characterize ice thickness in a relatively simple form. The two variables used to keep track of these levels are ice thickness, which is equivalent to the mass of ice in any grid cell, and ice concentration, which is defined as the fraction of the grid cell area covered by thick ice.

One issue is how the constraints on these two variables are imposed. In particular, assuming continuity of the ice thickness and the ice velocity, the ice thickness should remain non-negative (see [16, Theorem 3.10]). The non-negativity of ice concentration is similarly guaranteed. As will be demonstrated, preserving the non-negativity in both parameters is an important consideration for choosing a numerical solver. Moreover, although not explicitly providing a method to guarantee the upper bound of ice concentration to be $1$, such a constraint is described in Hibler [7] to be equivalent to adding a mechanical sink term in the model, which is turned on when the ice concentration reaches $1$ to prevent its further growth.

As far as we know, the numerical methods used to solve the model or impose the constraints have not been rigorously analyzed or compared. As already mentioned, the original work in [7] does not explicitly discuss how these constraints may be incorporated into the numerical implementation of the model. Since then some investigations have explicitly provided approaches both for numerical simulation of the model as well as for imposing these constraints. For example, Mehlmann [29] uses a finite element framework to solve the sea ice model and imposes restrictions on the trial spaces of ice thickness and ice concentration through a projection of the solution. Lipscomb and Hunke [24] adopt an incremental remapping scheme for sea ice transport. This is a Lagrangian approach that preserves the monotonicity by Van Leer limiting. That is, the gradients are reduced when necessary to ensure that the values in the reconstructed fields stay inside the range of the mean values in the cell and its neighbors. All of these mentioned approaches impose the model constraints through a post-processing procedure, which may, unfortunately, introduce discontinuity into the numerical solution, which affects the accuracy and might ultimately impact stability so that the solution does not converge.

Therefore, the other goal of this investigation is to develop a numerical approach that more intuitively imposes the model constraints without any post-processing procedure. To accomplish this task, we propose using the potential function method motivated by the analogous approach of using a double-well potential function in what is commonly referred to as the phase field method [5, 18], which is designed to solve interface problems by treating the interface as an object with finite thickness. Our proposed method, described in Section 4, offers a simple but elegant way to incorporate additional restrictions into the model and could be generalized in various settings. It further yields a modified transport model with the extra forcing terms coming from the potential energy function, which is consistent with Hibler’s statement on how to include the mechanical sink term in the model.

The rest of the paper is organized as follows. In Section 2 we describe the sea ice model. Section 3 provides a brief overview of standard numerical solvers for the sea ice model, focusing on the JFNK solver and EVP solver. We describe both the WENO scheme and the potential function method in Section 4. In Section 5 we conduct some numerical experiments and compare the performances of the WENO method with the more typically employed CD scheme. We also provide some numerical illustrations depicting the use of the potential function method. We make some concluding remarks in Section 6.

We illustrate the procedure of the JFNK implementation with a backward Euler time integration scheme for the momentum equation (2.2a). The time-discretized one-dimensional momentum equation is written as

where the superscript $n$ denotes the current time level. The numerical solution at the previous time level $n-1$ for (2.2) is known. Let ${\bf u}^{n}=\{u_{j}\}_{j=1}^{N}$ denote the approximation of $u^{n}$ obtained by some finite difference spatial discretization technique at each grid point $x_{j}$, $j=1,\dots N$.²²2To avoid cumbersome notation, for the rest of this paper we will use ${\bf u}$ to depict the vectorized solution of (2.2) (and not the continuous solution (2.1)). The spatial discretization scheme will be specified on both non-staggered and staggered grids in Section 4. For now we generically define the solution on $N$ grid points. At current time level $n$, we therefore seek a solution to

where $\mathbf{F}({\bf u}^{n})$ is the difference between the right- and left-hand sides of (3.1) following spatial discretization.

Since we are focusing on a single time step, we can simplify the notation by dropping the superscript $n$, so that we seek the solution ${\bf u}={\bf u}^{n}$. Using the velocity solution at the previous time level as the initial value ${\bf u}^{(0)}$, we iteratively solve a sequence of linearized systems to consecutively obtain ${\bf u}^{(1)},{\bf u}^{(2)},\cdots,{\bf u}^{(k)},\cdots$ until some stopping criterion is satisfied. Algorithm 1 summarizes the iterative technique. More details can be found in [19, 21, 3].

Observe that Algorithm 1 is an inexact Newton’s method as it approximates the Jacobian $\mathbf{J}$ by a first-order Taylor series expansion. The linear system of equations is, in general, solved by the preconditioned FGMRES method [31], which is a Krylov subspace method. This method is introduced as a matrix-free approach because forming and storing the Jacobian matrix is prohibitively expensive in CPU time and storage [21]. In the one-dimensional case, we are able to obtain the matrix representation of $\mathbf{J}$ by applying it to the basis vectors. We can then solve the linear system using a direct solver as the computational cost and efficiency in the one-dimensional case are not causes for concern.

In the EVP solver, the velocity at time level $n$ is obtained by explicitly solving the momentum equation from the previous time level $n-1$. In particular, the constitutive law was rewritten by Hunke and Dukowicz [9] to include a time dependence on the stress tensor. The velocity is then solved together with stress during subcycling. In the one-dimensional case [39], the stress-strain relationship

is equivalent to

By adding an artificial elastic strain with an elastic parameter $E$, we obtain

In the original version of the EVP model [9], the viscosities $\eta$ and $\zeta$ were held fixed throughout the subcycling procedure. However, because the viscosities were not regularly updated, such linearization of the internal stress term caused the computed principal stress states to lie outside the elliptical yield curve [8]. To address this issue, Hunke [8] proposed to include the viscosities within the subcycling, while simultaneously changing the definition of the elastic parameter $E$ to maintain the computational efficiency. Specifically, with $E$ defined in terms of a damping timescale for elastic waves and $T$ according to the equation $\displaystyle E=\frac{\zeta}{T}$, (3.4) can be rewritten as

The subcycling solution is advanced iteratively with subcycling time step $\Delta t_{e}$. This approach yields the time evolution of stress as a function of the velocity from the previous iterate according to

where the subcycling iterate is denoted with the superscript $s$. With the newly calculated stress in (3.6), the velocity is subcycled according to

Observe that in the EVP solver, the damping timescale $T$ is a tuning parameter satisfying $\Delta t_{e}<T<\Delta t$, which is in general set to be $T=0.36\Delta t$ following the documentation of the CICE model [11]. In addition, we denote the number of subcycles by $N_{sub}$, satisfying $N_{sub}\times\Delta t_{e}=\Delta t$.

A main goal of this investigation is to demonstrate the advantages of using higher-order methods to solve the sea ice model. We use the WENO method [27] as a prototype for two reasons. First, WENO is designed to have a higher-order convergence rate for smooth solutions than standard three point stencil CD schemes (which are second order), and second, WENO is able to maintain stable, non-oscillatory, and sharp discontinuity transitions so that it is suitable for sea ice with natural discontinuity feature of thickness and concentration. We verify that both of these properties hold in our numerical simulation of sea ice cover with and without sharp features.

We use the method of lines for time integration in each numerical test. To ensure stability and maintain the accuracy obtained in the spatial derivative approximation, we use the third-order TVD Runge-Kutta (TVRK3) [35] when employing both WENO and CD schemes. We then compare the numerical convergence properties of the results obtained using both methods. We implement the same time integration scheme to ensure that we are only comparing the spatial discretization performances for each scheme and not evaluating the time integration methods. Indeed, using an implicit time integration (backward Euler) for the momentum equation and explicit time integration (forward Euler) for the transport equations is common [22, 23, 3]. Our numerical example in Section 5.3 provides a case study for the mixed time integration approach. It is also, of course, possible to use a different spatial discretization for each equation in the system. As we did not observe any advantage in this approach, to reduce complexity, we did not do this.

In the WENO case we construct a non-staggered grid so that all variables are defined at the center of each grid cell. That is, we seek the solution at the $j=1,\dots,N$ midpoints of each cell, $\displaystyle x_{j-\frac{1}{2}}=x_{j}-\frac{1}{2}\Delta x$, with $x_{j}=j\Delta x$ and $\displaystyle\Delta x=\frac{L}{N}$ where $L$ is the domain length. Following [26] for nonlinear degenerate parabolic equations, we use higher-order finite differencing for (2.2a). In particular, $\displaystyle\frac{\partial u}{\partial x}$ and $\displaystyle\frac{\partial\sigma}{\partial x}$ are discretized using the fifth order finite difference WENO method for conservation laws [13] based on the left-biased stencil and right-biased stencil, respectively. The fifth order WENO method for conservation laws [13] is also used to solve the transport equations (2.2b) and (2.2c).

On the other hand, in the CD case, we construct a one-dimensional version of the staggered Arakawa C-grid [1], where the velocity $u$ is defined on vertices, $u_{0},\cdots,u_{N}$, and the traces $h$ and $A$ are defined at the center of each grid cell, $h_{\frac{1}{2}},\cdots,h_{N-\frac{1}{2}}$ and $A_{\frac{1}{2}},\cdots,A_{N-\frac{1}{2}}$ respectively. Correspondingly, the stress $\sigma$, the viscosities $\eta$ and $\zeta$ and the ice strength $P$ are also defined at the center of each grid cell. To solve for the velocity $u$ in the momentum equation (2.2a), we take $\displaystyle h_{i}=\frac{1}{2}(h_{i+\frac{1}{2}}+h_{i-\frac{1}{2}})$ for $i=1,\cdots,N-1$. We then approximate $\displaystyle\frac{\partial u}{\partial x}$ at each cell center as

so that $\displaystyle\sigma=(\eta+\zeta)\frac{\partial u}{\partial x}-\frac{P}{2}$ is defined at $x_{i+\frac{1}{2}}$. This leads to the approximation $\displaystyle\frac{\partial\sigma}{\partial x}$ at each vertex given by

Similarly, for the transport equation (2.2b) of $h$, we approximate $\displaystyle\frac{\partial}{\partial x}(uh)$ at each cell center as

We similarly solve $A$ using the same spatial discretization for the transport equation (2.2c). Finally, we note that although it is not a realistic assumption, we assume periodic boundary conditions to avoid errors introduced by boundary approximation. That said, the complexities introduced at the boundaries are not fundamentally different from other PDE models.

In terms of time integration, we consider the ordinary differential equation

where $L(u)$ is a discretization of the spatial operator. The TVRK3 method advances the current solution $u^{n}$ to the next time level $u^{n+1}$ according to Algorithm 2.

Due to its physical interpretation, the variable ice thickness $h$ in the sea ice model must remain non-negative. Similarly, the ice concentration value $A$ must be between 0 and 1. It is crucial for the numerical methods to preserve both of these properties to ensure that a meaningful solution is obtained. Motivated by the double-well potential function in the phase field method [5, 18], where the resulting equation is limited to a particular set of prescribed values due to the local minima of the potential function, we develop the potential function method here to impose the corresponding restrictions of ice thickness and ice concentration.

We begin by illustrating the potential function method on the transport equation of ice concentration $A$ and note that the case for ice concentration $h$ similarly follows. First, to restrict ice concentration $A$ so that $0\leq A\leq 1$, we define a potential function in a piecewise manner as

where $f_{1}>0$ and $f_{2}>0$ for all $A$, and $\gamma_{1}>0$ and $\gamma_{2}>0$ are parameters chosen so that $f$ has minima on $[0,1]$. For example, if both $f_{1}$ and $f_{2}$ are linear functions, a particular form of $f$ might be

The transport equation (2.2c) is then modified by adding a forcing term given by the gradient of the potential. This has the effect of the ice concentration experiencing a gradient force that tracks down to the physical range $[0,1]$. The resulting equation is given by

Observe that for the piecewise linear case, the forcing term $-f^{\prime}(A)$ is piecewise constant, meaning that the force transition is not continuous. To enable a more desirable smooth transition for this term we will instead choose both $f_{1}$ and $f_{2}$ to be quadratic and define $f$ as

The corresponding forcing term is now given by

which is clearly continuous.

To assess the convergence rate of WENO, we consider a test problem in a domain $\Omega=[0,\,2000]$ km with a known analytical solution.³³3This is, of course, generally not the case as the sea ice model has no known solutions. We construct our test cases by adding appropriate extra forcing terms to the governing equations. Specifically, we introduce to the right hand side of each equation of (2.2) forcing terms which are obtained by plugging into the left hand side terms corresponding to the following proposed solutions:

Observe that the values in (5.1) of $u_{\text{true}}$, $h_{\text{true}}$ and $A_{\text{true}}$ are consistent in magnitude to their corresponding true physical values. Also observe that while we construct extra forcing terms from the left hand side of (2.2a), the wind forcing is canceled out and therefore does not affect the convergence results. Initial conditions for the system are obtained by plugging $t=0$ into (5.1).

The total simulation time is $T=5$ s. The time step $\Delta t=10^{-4}$ s is intentionally chosen to be small enough to ensure that the time discretization error does not affect the convergence rates. It furthermore allows us to conduct convergence tests directly on the VP sea ice model for both the WENO and CD explicit spatial discretization schemes without the usual concern for the stability issue associated with explicit methods.

Table 2 compares the relative $\ell_{2}$ errors for increasing resolutions with each spatial discretization choice. We observe second-order convergence for all three variables for the CD case, which is consistent with the standard convergence analysis results for CD schemes. By employing the WENO scheme in the ideal case, one would expect to obtain sixth-order convergence for the velocity $u$ and fifth-order convergence for both the ice thickness $h$ and ice concentration $A$. However, due to the complexity and non-linearity of the sea ice model, coupled with the fact that the added extra forcing terms are not being updated in the stages of TVRK3 time integration, theoretical accuracy is unlikely to be obtained. We still observe higher-order convergence for all three variables as compared to the CD results. In addition, a direct comparison of the error magnitudes for all three variables indicates that the WENO scheme indeed provides more accurate results than those obtained using CD, noting that WENO appears to be mainly affected by round-off error at $10$ km resolution.

In this example we test the performance of the WENO scheme on a simulation of a sea ice cover with sharp features. To better capture the solution behavior near the discontinuity region while maintaining periodic boundary conditions, the structure of ice is designed such that relatively solid ice covers both ends of the domain and a very thin layer of ice is in the center of the domain. This is realized via a discontinuous setting on the initial conditions of ice thickness and ice concentration given by

For the external forcing in (2.3) we impose uniform constant wind forcing $u_{a}=10$ m/s.

As the WENO and CD schemes yield theoretically different convergence rates, for a direct comparison, we also consider the linear WENO scheme, for which the nonlinear weights are replaced by linear ones of the same order accuracy. Note that this is equivalent to using an upstream centered scheme (upstream in time, centered in space), [13]. Due to the combined stencil, the highest possible order of accuracy is obtained in smooth regions. The results are oscillatory near discontinuities, however.

We now test the model for which no exact solution is known. The main goal of this numerical test is to determine how out-of-range issues, namely $A<0$, $A>1$, or $h<0$, may be effectively mitigated using the potential method described in Section 4.2.

We use the backward Euler time integration with CD spatial discretization and JFNK to solve the momentum equation, along with TVRK3 with CD spatial discretization scheme for transport equations. We note that in choosing to use an implicit time-stepping method for solving (2.2a) we avoid issues concerning stability. Moreover, with regard to the transport problem, we note that the WENO scheme does not yield out-of-range negative solutions for either $A$ or $h$. Hence to determine the efficacy of the potential method we apply the CD scheme (using explicit time-stepping) to the transport equations (2.2b) and (2.2c). As discussed in 4.1, we also use the one-dimensional version of the staggered Arakawa C-grid [1].

We run all the experiments on a domain of length $2000$ km for a total integration time of 6 days. We choose a spatial resolution of $20$ km and $90$ s as the time step. The initial conditions are constant throughout $[0,2000]$ km and given by

We impose Dirichlet boundary conditions so that

Observe that since the ocean is considered to be at rest in this model, the only variable external forcing term is the wind, which is imposed uniformly as a constant given by

For imposing the boundary conditions numerically, we have $u_{0}=0$ (similarly $u_{N}=0$) so that

leading to

for the transport equation (2.2b) (similarly for $\{d(uh)\}_{N-\frac{1}{2}}$). The Dirichlet boundary conditions also yield analogous equations for $A$ at the boundaries in (2.2c).

For comparison purposes we first run the simulation without applying the potential function method. Throughout the 6-day simulation, we find that the largest $A$ value is 1.0540, the smallest $A$ value is -0.1445, and the smallest $h$ value is -0.1606, which are all out of range.

To employ the potential function method, we first simulate the model (without incorporating the potential function method into the transport equations) until the time $T_{1}$ for which $\displaystyle\min_{x}\{A(x,{T_{1}})\}<0$.⁴⁴4We only do this process one time for each out-of-range situation, $A<0$, $A>1$, and $h<0$. The potential function parameters then remain fixed for all time. The potential function variables corresponding to (4.7) are then determined as

While $B_{0}$ is determined directly from the numerical implementation of the scheme, as previously mentioned after (4.6), $a$ is computed as a linear approximation of $\displaystyle\frac{\partial u}{\partial x}$.

Following the discussion in Section 4.2, we then find a uniform bound for $\gamma_{1}$ by computing

The process is similar for determining $\gamma_{2}$ for when $\displaystyle\max_{x}\{A(x,{T_{2}})\}>1$, ${T_{2}}>0$, and in our experiment we obtain

Finally the corresponding range for the parameter $\gamma$ associated with $h<0$ is

Based on these results and the discussion in Section 4.2, we choose

We then run the remainder of the simulation, up until final time $T=6$ days, with the potential function method now incorporated into the transport equations.

The solutions of $u$, $h$, and $A$ are displayed in Figure 3, where the top row depicts the solutions without employing the potential function method, the middle row shows the solutions when the potential function is used, and the bottom row shows the solution plots at the final time. It is evident that using the potential function yields some differences in the solution. The change of color shade in the velocity image around the 2000 km boundary of the region at roughly 30 hours indicates the situation when the non-physical values are detected and potential function starts to take effect. With the potential function method, the velocity remains a relatively large value on a larger portion of the domain and decreases to 0 in a sharper manner towards the end of the boundary. The behavior of the ice thickness near the boundary, especially around 2000 km, dramatically changes due to the application of the potential function method. In particular, the non-physical negative values near the left boundary are replaced by smoother physically meaningful values. Also, the thickness value is much larger near the right boundary, and the ridging effect is more clearly observed. This makes more sense physically for ridging on the ocean-land boundary. For the ice concentration, non-physical negative values near the left boundary and non-physical large values near the right boundary are also appropriately treated by the potential function method.