An adaptive Newton multigrid method for a model of marine ice sheets

Guillaume Jouvet, Carsten Gräser
2013 Journal of Computational Physics  
In this paper, we consider a model for the time evolution of three-dimensional marine ice sheets. This model combines the Shallow Ice Approximation (SIA) for the ice deformation, the Shallow Shelf Approximation (SSA) for the basal sliding, and the mass conservation principle. At each time step, we solve a scalar p-Laplace minimization-type problem with obstacle (SIA), a vectorial p-Laplace minimization-type problem (SSA) and a transport equation (mass conservation). The two minimization
more » ... are solved using a truncated nonsmooth Newton multigrid method while the transport equation is solved using a vertex-centred finite volume method. Our approach is combined to an heuristic mesh adaptive refinement procedure to face the large gradients of the solution that are expected between the ice sheet and the ice shelf. As applications, we present some simulations of the Marine Ice Sheet Model Intercomparison Project MISMIP (2D and 3D) and validate our results against an analytic solution (2D) and other participant model results (3D). Further numerical results show that the convergence of our Newton multigrid method is insensitive to local refinements making our overall adaptive strategy fully efficient. friction on the bedrock can be described by or Coulomb laws [36] while perfect sliding occurs where ice is floating. In fact, only the vertical shear (resp. longitudinal) components of the stress tensor are significant in most of shallow ice sheets (resp. shelves) [43] . Simplified shallow ice models -that are mostly used by glaciologists -are derived from such observations. More precisely, simplifications occur by neglecting the terms where the small aspect ratio of ice sheets (in the order of 10 −3 ) appears in the dimensionless Stokes equations. As a result, the Shallow Ice Approximation (SIA) model [21] accounts only for the vertical shear. Similarly, the Shallow Shelf Approximation (SSA) model [36] accounts for the longitudinal stresses and the friction on the bedrock. In the last ten years, more or less simplified mechanical models have been involved in marine ice sheets models. Vertically-integrated models are the most popular ones since they are two-dimensional, which considerably reduces the size and the complexity of the system to solve. Among these models, those based on the SSA are the most common [9, 15, 32, 37, 42, 29] . Some attempts to include the vertical shearing into twodimensional (still vertically-integrated) models have been made by coupling the SIA to the SSA (independently) [43] or in a more sophisticated way using the Schoof-Hindmarsh model [8] . Finally, the most accurate models are three-dimensional like the First Order Approximation (FOA) [5, 6, 29] or the full Stokes model [10, 29] . Let us notice that the SSA, the FOA and the Stokes models have been combined using a domain decomposition in [39] . Among the above contributions, one can distinguish two strategies for the numerical treatment of the GL. Since the grounded and floating domains evolve over time, one can follow the GL exactly or use an approximation. In both cases, the floating domain is determined by a geometrical criterion resulting from Archimedes principle. In the first case, a criterion determines how to move the mesh together with the GL [9, 37] . Each domain (floating or grounded) involves a different model and a boundary condition imposes the continuity of stresses at the GL. Moving a one-dimensional mesh-that corresponds to a two-dimensional ice sheet-is an easy task since the GL is a single point. However, this becomes harder with one more horizontal dimension since the GL is a curve [34] . To overcome such difficulty, the second approach considers a unified model for the ice sheet and for the ice shelf with an implicit description of the GL through the flotation criterion. This criterion designs which points of the mesh are in the floating part and which points are in the grounded part such that different sliding conditions can be applied. Such an approach must be combined with an adaptive mesh refinement procedure to deal with the sharp changes in the ice flow regime that are expected close to GL [34]. The implementation of moving or adaptive meshes is strongly connected to the numerical method, which is used to solve the mechanical equations. Existing shallow ice sheet and ice shelf models are based on the three most popular numerical techniques: finite differences [9,32,42,43], finite volumes [8] and finite elements [15, 10] . It is well known that these techniques do not offer the same flexibility in term of re-meshing. Indeed, finite differences usually require structured meshes. In contrast, finite element and finite volume techniques can deal more easily with locally refined and unstructured meshes. In this paper, we implement an isothermal ice sheets and ice shelves model that is based on a linear superposition of the SIA model and of the SSA model. Such superposition assumes that the modeling of the horizontal and the vertical ice flows can be linearly decoupled. This assumption is physically justified for most ice sheets and ice shelves [43] . A weighted combination of these two models was originally proposed in [7] and later modified into a simple superposition in [43] . There are two ways of superposing SIA and SSA velocity components. The simplest way consists of summing both components such that the mass conservation equation remains a pure advection problem. The second approach, which is adopted in this paper, includes directly the explicit SIA component in the mass conservation, which acts as an additional nonlinear diffusion term due to the gradient involved in the SIA flux [25] . Diffusion and advection are decoupled by a first order splitting algorithm, such that it remains to solve at each time step two p-Laplace problems (one for the SSA and one for the SIA) and one transport equation. On one hand, both p-Laplace problems are approximate using the Truncated Nonsmooth Newton MultiGrid (TNMMG) method [17, 18] , which consists of a nonlinear relaxation scheme combined with Newton-type multigrid corrections. On the other hand, a vertex-centred finite volume method is advocated to solve the convection part of the model. Finally, we implement a heuristic adaptive local mesh refinement procedure to increase the accuracy close to GL. Our model has been shortly presented in [26] for a non-realistic two-dimensional ice sheet, which writes as a onedimensional mathematical problem. The goal of this paper is to extend the model to real three-dimensional marine ice sheets. Adding one horizontal dimension to the model obliges us to revise the meshing procedure. In particular, our solver relies on a geometrical multigrid approach, which requires regular, conform and -most importantly -hierarchical meshes. The construction of such one-dimensional adaptive meshes is straightforward. In contrast, further attention is required to not destroy the properties of a two-dimensional mesh such that efficiency of the geometrical multigrid method remains intact through the adaptive mesh refinement process. The reduction of computational costs by implementing an adaptive refinement procedure is another challenge of three-dimensional models. Finally, this paper dedicates a large part to applications and numerical results. In particular, we apply our method to the test problems proposed by the Marine Ice Sheet Model Intercomparison Project (MISMIP) in two and three dimensions [34, 33] . This allows us to compare our results to the ones of the other participating models and the efficiency of our adaptive re-meshing strategy. In particular, we verify the ability of the adaptive mesh procedure to reproduce a reversible GL after being perturbed in the ice flow and sliding parameters. As an additional outcome, we point out a limitation of the superposed "SSA + SIA" model in describing the GL behaviour, which is a consequence of the uncoupling between the SSA and SIA velocity components. Finally, we provide some numerical results concerning the performance of our nonlinear multigrid solver in term of convergence and computational costs.
doi:10.1016/j.jcp.2013.06.032 fatcat:svbv2jlorfgbvocgjxupekgxxm