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

Accurate modeling of sea ice dynamics is critical for predicting environmental variables and is important in applications such as navigating ice breaker ships. Research for both modeling and simulating sea ice dynamics is ongoing, with the most widely accepted model based on the viscous-plastic (VP) formulation introduced by Hibler in 1979. Due to its highly nonlinear features, this model is intrinsically challenging for computational solvers. In particular, sea ice simulations often significantly differ from satellite observations. This study therefore focuses on improving the numerical accuracy of the VP sea ice model. Since the poor convergence observed in existing numerical simulations stems from the nonlinear nature of the VP formulation, this investigation proposes using the celebrated weighted essentially non-oscillatory (WENO) scheme -- as opposed to the frequently employed centered difference (CD) scheme -- for the spatial derivatives in the VP sea ice model. We then proceed to numerically demonstrate that WENO yields higher-order convergence for smooth solutions, and that furthermore it is able to resolve the discontinuities in the sharp features of sea ice covers -- something that is not possible using CD methods. Finally, our proposed framework integrates a potential function method that utilizes the phase field method to naturally incorporates the physical restrictions of ice thickness and ice concentration in transport equations, resulting in a modified transport equations which includes additional forcing terms. Our method does not require post-processing, thereby avoiding the possible introduction of discontinuities and corresponding negative impact on the solution behavior. Numerical experiments are provided to demonstrate the efficacy of our new methodology.