Discretization errors and free surface stabilization in the finite difference and marker-in-cell method for applied geodynamics: A numerical study

The finite difference–marker‐in‐cell (FD‐MIC) method is a popular method in thermomechanical modeling in geodynamics. Although no systematic study has investigated the numerical properties of the method, numerous applications have shown its robustness and flexibility for the study of large viscous deformations. The model setups used in geodynamics often involve large smooth variations of viscosity (e.g., temperature‐dependent viscosity) as well large discontinuous variations in material properties (e.g., material interfaces). Establishing the numerical properties of the FD‐MIC and showing that the scheme is convergent adds relevance to the applications studies that employ this method. In this study, we numerically investigate the discretization errors and order of accuracy of the velocity and pressure solution obtained from the FD‐MIC scheme using two‐dimensional analytic solutions. We show that, depending on which type of boundary condition is used, the FD‐MIC scheme is a second‐order accurate in space as long as the viscosity field is constant or smooth (i.e., continuous). With the introduction of a discontinuous viscosity field characterized by a viscosity jump (η*) within the control volume, the scheme becomes first‐order accurate. We observed that the transition from second‐order to first‐order accuracy will occur with only a small increase in the viscosity contrast (η* ≈ 5). We have employed two methods for projecting the material properties from the Lagrangian markers onto the Eulerian nodes. The methods are based on the size of the interpolation volume (4‐cell, 1‐cell). The use of a more local interpolation scheme (1‐cell) decreases the absolute velocity and pressure discretization errors. We also introduce a stabilization algorithm that damps the potential oscillations that may arise from quasi free surface calculations in numerical codes that employ the strong form of the Stokes equations. This correction term is of particular interest for topographic modeling, since the surface of the Earth is generally represented by a free surface. Including the stabilization enables physically meaningful solutions to be obtained from our simulations, even in cases where the time step value exceeds the isostatic relaxation time. We show that including the stabilization algorithm in our FD stencil does not affect the convergence properties of our scheme. In order to verify our approach, we performed time‐dependent simulations of free surface Rayleigh‐Taylor ins