Improving Numerical Accuracy for the Viscous-plastic Formulation of Sea Ice

By Tongtong Li, Anne E. Gelb, and Yoonsang Lee

Sea ice dynamics play a critical role in shaping the ice cover in polar regions. The accurate modeling and simulation of such dynamics are crucial for the prediction of environmental variables, which in turn find a wide range of practical applications — such as the navigation of icebreaker ships. As such, many researchers have extensively studied modeling of sea ice dynamics. For example, W.D. Hibler introduced the viscous-plastic (VP) sea ice model in 1979 [2]. This model, which consists of a nonlinear momentum equation and two transport equations, forms the basis of the most widely accepted model in the present day. Yet due to the momentum equation’s highly nonlinear features that stem from the VP material law, the VP model is intrinsically challenging to solve numerically. And now that finer mesh resolutions are available, it is increasingly apparent that solutions obtained by traditional low-order numerical methods are not well resolved; indeed, they exhibit a significant discrepancy with obtainable satellite observations. This discrepancy motivates us to improve the numerical accuracy of sea ice representation based on the one-dimensional (1D) VP model. With the 1D results in hand, we anticipate that our approach will prove useful in the more complicated two-dimensional case.

To tackle the difficulties that arise from the nonlinearity in the momentum equation for the Hibler model, researchers have designed several approaches that generally fall into two categories. First, the employment of implicit time-stepping schemes ensures numerical stability, but implementing iterative methods for nonlinear problems is notoriously difficult. Multiple investigations have proposed methods to improve the efficiency of implicit numerical solvers and their corresponding preconditioners. For example, the celebrated Jacobian-free Newton-Krylov (JFNK) method is computationally faster than the standard Picard solver and suitable for parallelization [7]. But the JFNK solver is not robust, as its failure rate increases with more refined spatial resolution. 

On the flip side, the second category of approaches entirely avoids implicit methods by adding an artificial elastic term to the VP constitutive equation, therefore resulting in the so-called elastic-viscous-plastic (EVP) model [3, 4]. Although the relaxed stability condition enables the use of explicit time stepping, notable differences between the EVP approximate solution and the reference solution become relatively more distinct with finer resolution [6]. Finally, neither the JFNK nor the EVP approach explicitly address the physical constraints for ice thickness and ice concentration that are determined by the transport equations. Specifically, the numerical solver should preserve the non-negativity of the ice thickness and maintain a range of \([0,1]\) for the ice concentration.

We propose numerical techniques to help mitigate the limitations of current numerical solvers for the Hibler model [8]. In particular, we demonstrate the advantages of higher-order methods and use the weighted essentially non-oscillatory (WENO) scheme as a prototype [9]. The WENO method achieves higher-order accuracy in smooth regions while sharply resolving discontinuities in an essentially non-oscillatory fashion. These properties are highly desirable for the sea ice model because (i) low-order methods yield poor convergence in the case of smooth solutions and (ii) sea ice covers often contain sharp features. 

Figure 1. Simulation of sea ice with sharp features on viscous-plastic (VP) model. 1a. Velocity plot using weighted essentially non-oscillatory (WENO) scheme at one hour. 1b. Velocity plot using linear WENO at 2,000 seconds. 1c. Velocity plot using centered-difference (CD) scheme at 2,332 seconds. Figure courtesy of [8].

We verify WENO’s superior performance for both cases as compared to the more common centered-difference (CD) scheme and linear WENO scheme.¹ We observe that WENO yields more accurate and higher-order convergence for smooth solutions and also corroborate WENO’s capacity to resolve discontinuities (see Figure 1). The latter involves designing the ice structure in the model so that relatively solid ice covers both ends of the domain, with a very thin layer of ice in the center. Figure 1 shows that only WENO maintains a sharp, non-oscillatory solution for the velocity in each sub-region; the linear WENO velocity profile includes oscillations, while a sharp overshoot occurs in the CD case. Although we plot the WENO solution at the final time of one hour, we respectively present the solutions for the linear WENO and CD at 2,000 and 2,332 seconds, as the oscillations eventually cause these solutions to blow up.

Next, we introduce the potential function method to further ensure that our computational technique preserves the physical constraints on ice thickness and concentration. The double-well potential function in the phase field method—wherein the resulting equation is limited to a particular set of prescribed values due to the potential function’s local minima—motivates this idea [1, 5]. For both the ice thickness and concentration variables, we define corresponding potential functions in a piecewise manner so that their minima fall on physically meaningful ranges. We then modify the transport equation accordingly by adding a forcing term that is given by the gradient of the potential. As a result, the variable experiences a gradient force that tracks down to the physical range.

Figure 2 numerically illustrates the benefits of the potential function method by comparing the solutions of ice velocity, ice thickness, and ice concentration at the final simulation time (144 hours) both with and without the potential method. Noticeable distinctions in the results include the observation that the value of velocity remains larger for a greater portion of the domain, then sharply decreases to zero near the boundary when using the potential method. In addition, the ice thickness’ behavior with the potential method changes dramatically when approaching the right boundary. This outcome represents a ridging effect, which is consistent with ocean-land boundary behavior. The nonphysical negative solution at the left boundary is also replaced by something that is more physically meaningful. Finally, the potential function method effectively mitigates nonphysical negative values near the left boundary and nonphysical large values near the right boundary for the ice concentration.

Figure 2. An example of the potential function method with solution plots at 144 hours. 2a. Solution plot for velocity. 2b. Solution plot for ice thickness. 2c. Solution plot for ice concentration. Figure courtesy of [8].

In summary, our research explores the limitations of the current methodology for solving the VP sea ice model — namely, poor convergence and out-of-range issues for both ice concentration and thickness. To improve the numerical solvers’ performance, we propose the use of higher-order methods, with WENO as a prototypical example. We provide a case study for the 1D Hibler model and verify that WENO exhibits better numerical convergence properties than commonly employed numerical algorithms. Moreover, WENO can resolve discontinuities and sharp features that may occur in sea ice covers. With regard to the out-of-range issue, our investigation develops and implements a potential function method that naturally incorporates the physical restrictions of ice thickness and ice concentration in the transport equations.

The current work is restricted to a 1D case study, since its numerical convergence properties are relatively straightforward. Future investigations will employ both higher-order methods and the potential function approach in more realistic environments, as well as within the data assimilation framework.

¹ The linear WENO method is WENO with “linear weights”, which could provide the highest possible order of accuracy when the solution is smooth. We use it here to represent more traditional higher-order methods that require smoothness in the underlying solution.

References
[1] Fix, G. (1983). Phase field models for free boundary problems. In A. Fasano & M. Primicerio (Eds.), Free boundary problems: Theory and applications (Vol. 2, pp. 580-589). Boston, MA: Pitman Advanced Publishing Program.
[2] Hibler, W.D. (1979). A dynamic thermodynamic sea ice model. J. Phys. Oceanogr., 9(4), 815-846.
[3] Hunke, E.C. (2001). Viscous-plastic sea ice dynamics with the EVP model: Linearization issues. J. Comput. Phys., 170(1), 18-38.
[4] Hunke, E.C., & Dukowicz, J.K. (1997). An elastic-viscous-plastic model for sea ice dynamics. J. Phys. Oceanogr., 27(9), 1849-1867.
[5] Langer, J.S. (1986). Models of pattern formation in first-order phase transitions. In G. Grinstein & G. Mazenko (Eds.), Directions in condensed matter physics (pp. 165-186). Singapore: World Scientific. 
[6] Lemieux, J.-F., Knoll, D., Tremblay, B., Holland, D., & Losch, M. (2012). A comparison of the Jacobian-free Newton-Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: A serial algorithm study. J. Comput. Phys., 231(17), 5926-5944.
[7] Lemieux, J.-F., Tremblay, B., Sedláček, J., Tupper, P., Thomas, S., Huard, D., & Auclair, J.P. (2010). Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton-Krylov method. J. Comput. Phys., 229(8), 2840-2852.
[8] Li, T., Gelb, A., & Lee, Y. (2022). Improving numerical accuracy for the viscous-plastic formulation of sea ice. Preprint, Arxiv:2206.10061.
[9] Liu, X.-D., Osher, S., & Chan, T. (1994). Weighted essentially non-oscillatory schemes. J. Comput. Phys., 115(1), 200-212.

	Tongtong Li is a postdoctoral researcher in the Department of Mathematics at Dartmouth College working on a Department of Defense/Office of Naval Research-funded Multidisciplinary University Research Initiatives (MURI) project #N00014-20-1-2595 entitled "The Integrated Foundation of Sensing, Modeling, and Data Assimilation for Sea Ice Prediction." She is interested in numerical analysis and solutions of partial differential equations, as well as data assimilation with applications in fluid dynamics and sea ice modeling.
	Anne E. Gelb is John G. Kemeny Parents Professor of Mathematics at Dartmouth College. She is a numerical analyst focusing on high order methods for signal and image restoration, classification, and change detection for real and complex signals from temporal sequences of collected data.
	Yoonsang Lee is an assistant professor in the Department of Mathematics at Dartmouth College. His research focuses on applied mathematics and computational issues in prediction and uncertainty quantification of complex dynamical systems.