Variational Stokes: free surface for staggered grid finite differences schemes for geodynamic and ice flow modelling

In geodynamic models, the incorporation of a free surface boundary condition is crucial to better understand and resolve, for example, the coupled interaction between the lithosphere, deep mantle, and surface processes. A free surface is similarly important in ice flow models to capture the geometry of the ice sheet and its temporal evolution.

While the implementation of a free surface is relatively straightforward in finite element methods due to their capacity to manage complex and deformable geometries, and boundary conditions, the treatment of a free surface boundary (or any other boundary) that is not aligned with the staggered grid of a finite difference (FD) scheme poses a significant challenge. A common approach adopted by FD geodynamic codes, such as I2/3VIS (Gerya and Yuen, 2007) or LAMEM (Popov and Kaus, 2016), involves the incorporation of an additional rheology layer above the surface, simulating the presence of air. However, the of Stokes solvers is constrained by the viscosity contrast occurring within the domain, typically in the range of 6 to 7 orders of magnitude. Consequently, the viscosity of the air layer is limited to values within the 1e16-1e18 Pa*s range. This approach is not only physically unrealistic, but also leads to an inaccurate topography evolution, and introduces a very strong and sharp viscosity contrast at the rock and air interface, which hinders the convergence of iterative solvers particularly the Accelerated Pseudo-Transient (APT) method employed in this study.

To address these limitations, we propose the implementation of a variational Stokes approach (Larionov et al. 2017), which allows for the incorporation of both real free surface and solid wall boundary conditions. This approach is then combined with either a marker chain or a level set to track evolution of the surface. We demonstrate that this approach greatly improves the convergence rate of the iterative APT solver, as well as demonstrate its accuracy and applicability to geodynamic and ice flow simulations with a set of benchmarks and toy codes.

 

Gerya, Taras V., and David A. Yuen. "Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems." Physics of the Earth and Planetary Interiors 163.1-4 (2007): 83-105.

Popov, Anton, and Boris Kaus. "3D modelling of non-linear visco-elasto-plastic crustal and lithospheric processes using LaMEM." EGU General Assembly Conference Abstracts. 2016.

Larionov, Egor, Christopher Batty, and Robert Bridson. "Variational stokes: A unified pressure-viscosity solver for accurate viscous liquids." ACM Transactions on Graphics (TOG) 36.4 (2017): 1-11.