A second-order accurate in time IMplicit–EXplicit (IMEX) integration scheme for sea ice dynamics
Journal of Computational Physics
11 Current sea ice models use numerical schemes based on a splitting in time 12 between the momentum and continuity equations. Because the ice strength 13 is explicit when solving the momentum equation, this can create unrealis-14 tic ice stress gradients when using a large time step. As a consequence, 15 noise develops in the numerical solution and these models can even become 16 numerically unstable at high resolution. To resolve this issue, we have imple-17 mented an iterated
... ted IMplicit-EXplicit (IMEX) time integration method. This 18 IMEX method was developed in the framework of an already implemented 19 Jacobian-free Newton-Krylov solver. The basic idea of this IMEX approach 20 is to move the explicit calculation of the sea ice thickness and concentration 21 inside the Newton loop such that these tracers evolve during the implicit 22 integration. To obtain second-order accuracy in time, we have also modified 23 the explicit time integration to a second-order Runge-Kutta approach and 24 by introducing a second-order backward di↵erence method for the implicit 25 integration of the momentum equation. These modifications to the code are 26 minor and straightforward. By comparing results with a reference solution 27 obtained with a very small time step, it is shown that the approximate so-28 lution is second-order accurate in time. The new method permits to obtain 29 the same accuracy as the splitting in time but by using a time step that is 30 10 times larger. Results show that the second-order scheme is more than five 31 times more computationally e cient than the splitting in time approach for 32 an equivalent level of error. 33 ⇤ Keywords: sea ice, IMEX method, backward di↵erence, Newton-Krylov 35 method, numerical accuracy 36 To the best of our knowledge, all sea ice model time integration schemes 11 are based on a splitting in time between the momentum and the continuity 12 equations (e.g., [1, 2, 3, 4, 5]). This means that when solving the momentum 13 equation, the thickness distribution (including the amount of open water) is 14 held constant at the previous time level (it, however, varies spatially). Once 15 the velocity field is obtained, the thickness distribution is advanced to the 16 next time level. Furthermore, an operator splitting approach is generally 17 used to separate the change of the thickness distribution associated with 18 advection and the growth/melt related to thermodynamic processes (e.g., 19 [2, 3]). This paper focuses on dynamics and we therefore only discuss the 20 solution of the momentum equation and of the continuity equation without 21 the thermodynamic source terms. 22 23 Current sea ice model numerical schemes su↵er from significant numerical 24 issues. First, as explained by Lipscomb et al. , the splitting in time ap-25 proach leads to noise in the numerical solution and can even make the model 26 numerically unstable. As an illustrative example, consider ice converging 27 toward a coast due to an onshore wind; a stress gradient, associated with 28 an ice strength gradient, develops to oppose the wind stress. When using a 29 large time step with the splitting in time approach, an unrealistically large 30 ice strength gradient can occur. The stress gradient force can then overcom-31 pensate the wind stress and cause an unrealistic reversal of the flow (the ice 32 2 then diverges at the coast). This instability, fundamentally numerical, can 33 be cured by reducing the time step. Unfortunately, this obviously increases 34 the total computational time. Lipscomb et al.  proposed a modification to 35 the ridging scheme in order to mitigate this problem. 36 37 A second numerical issue is related to the solution of the momentum 38 equation. The rheology term, which determines the deformations of the ice 39 cover based on the internal ice stresses, causes the momentum equation to 40 be very nonlinear. Indeed, the VP rheology leads to a large change in the 41 internal stresses when going from a slightly convergent flow to a slightly di-42 vergent one (same idea for shear stresses). The current numerical solvers for 43 the momentum equation, however, have di culties in finding the solution of 44 this very nonlinear problem. There are two main classes of schemes to solve 45 the momentum equation: the implicit solvers, which involve an outer loop 46 iteration (sometimes referred to as Picard iteration, [5, 6, 7]) and the ones 47 based on the explicit solution of the momentum equation using the Elastic-48 VP approach [8, 9]. Both of these approaches, however, lead to a very slow 49 convergence rate [7, 9] if they converge at all [9, 10]. Because of this slow con-50 vergence rate, it is typical to perform a small number of Picard iterations or 51 of subcycling iterations. The approximate solution therefore contains resid-52 ual errors which are carried on in the time integration. 53 54 To resolve this slow convergence rate issue, Lemieux et al.  developed 55 a Jacobian-free Newton-Krylov (JFNK) implicit solver. They showed that 56 the JFNK solver leads to a more accurate solution than the EVP solver  57 and that it is significantly more computationally e cient than a Picard ap-58 proach . Following the work of Lemieux et al. , Losch et al.  have 59 recently developed a parallel JFNK solver for the MIT general circulation 60 model with sea ice . The numerical approaches of Lemieux et al.  and 61 Losch et al. , however, still rely on the splitting in time scheme and are 62 therefore susceptible to exhibit the numerical instability issue. 63 64 It is the purpose of this paper to introduce a fast and accurate time in-65 tegration scheme that resolves the instability associated with the splitting 66 in time approach. One possibility would be to solve fully implicitly the mo-67 mentum and continuity equations. This avenue would imply significant mod-68 ifications to the code and would be quite complex to implement. Instead, 69 the splitting in time issue is cured by using an iterated IMplicit-EXplicit 70 3 (IMEX) approach when solving the momentum and continuity equations. 71 This approach is built around our existing JFNK solver. Basically, the idea 72 is to move the explicit calculation of the thickness distribution inside the 73 implicit Newton loop. We take this approach one step further by modifying 74 the time integration in order to get second-order accuracy in time for the full 75 system. To do so, we introduce a second-order Runge-Kutta scheme for the 76 advection operation and discretize in time the momentum equation using a 77 second-order backward di↵erence (as in ). This paper is inspired by the 78 work of [14, 15] on an iterated IMEX method for radiation hydrodynamics 79 problems. 80 81 The main contribution of this paper is the development and demonstra-82 tion of a first-of-a-kind second-order accurate in time iterated IMEX inte-83 gration scheme for sea ice dynamics. This manuscript also shows the gain 84 in accuracy and computational time of the second-order IMEX method com-85 pared to the common first-order integration scheme based on the splitting in 86 time. 87 88 It is worth mentioning that some authors have recently questioned the 89 validity of the VP rheology. Sea ice models based on a VP rheology do not 90 capture the largest deformations events  and statistics of simulated de-91 formations do no match observations  in both space and time . While 92 some authors propose new and very di↵erent formulations of ice interactions 93 [18, 19], others claim that a VP rheology with modified yield curve and flow 94 rule can adequately represent the sea ice deformations . These new physi-95 cal parameterizations, under evaluation, also lead to very nonlinear problems 96 which would also clearly benefit from the availability of reliable and e cient 97 numerical schemes. 98 99 This paper is structured as follows. Section 2 describes the sea ice mo-100 mentum equation with a VP formulation and the continuity equation. In 101 section 3, the discretization of the momentum and continuity equations and 102 the descriptions of the standard splitting in time and new IMEX integration 103 schemes are presented. In section 4, more information about the model is 104 given. The description of the experiments and the results are outlined in 105 section 5. A discussion and concluding remarks are provided in section 6. 106 107 4 157 by Hibler , this formulation for ⇣ is continuously di↵erentiable. 158 159 694 695 Another improvement would be to replace our di↵usive first-order in space 696 upstream scheme by a more sophisticated advection operator. For example, 697 second-order accuracy in space could also be achieved by using the remap-698 ping scheme of Lipscomb and Hunke  . Note that a stabilization method 699 (di↵erent time-stepping approach) may be required as higher order advection 700 schemes are less di↵usive than a first-order upstream operator. 702 The JFNK solver is remarkably robust in longer simulations (five years). 703 At 40-km resolution, JFNK did not fail for either the SIT or the BDF2-704 IMEX-RK2 integration scheme. At 20-km resolution, convergence was not 705 reached on rare occasions for both integration schemes. With SIT, JFNK 706 had a failure rate as low as 0.027 % while JFNK with the BDF2-IMEX-RK2 707 scheme failed for only 0.025 % of the time levels (this is slightly smaller than 708 for SIT but probably not statistically significant). 709 710 Even though these failure rates are very small and when a failure occurs 711 it usually a↵ects only a few grid cells (not shown), the increase in the failure 712 25 rates with resolution indicates that further work is needed to improve the ro-713 bustness. A more sophisticated approach than the linesearch method might 714 help (e.g. ) but we also suspect that our preconditioning approach might 715 need to be revisited as we refine the grid. 716 717 Indeed, as the spatial resolution increases, the rheology term makes the 718 problem more and more nonlinear. We have observed occasional failures of 719 the preconditioned FGMRES at 10-km resolution for a linear tolerance of 720 0.1. To improve our preconditioning operator, we are currently working on 721 using the MultiLevel (ML) preconditioner from the Trilinos library . It is 722 possible, however, that this might not be su cient and that we might have to 723 reconsider the use of the Picard matrix for the preconditioning step. In other 724 words, our preconditioning matrix might have to be closer to the Jacobian 725 matrix than what the Picard matrix is. 726 727 This study was done using a serial code. Losch et al.  have recently 728 implemented a parallel JFNK solver for sea ice dynamics. They have demon-729 strated that the scaling of JFNK with a similar line relaxation approach for 730 the preconditioner is almost as good as for other solvers (Picard and EVP); 731 in their case for domain decompositions of up to 1000 CPUs. There is no 732 reason to believe that our BDF2-IMEX-RK2 approach would not exhibit 733 similar performances as the additional thickness and concentration calcula-734 tions performed in the Newton loop are explicit and do not require extra 735 communication overheads. Using a di↵erent preconditioner (such as ML) 736 might lead to an improved scalability of JFNK. This is the subject of future 737 work. 738 739 Appendix A: Undamped oscillation with a Crank-Nicolson approach 740 By centering in time (at n 1 2 ) the terms in the momentum equation, a 741 Crank-Nicolson approach also leads to second-order accuracy (not shown).