The finite difference-marker-in-cell (FD-MIC) method is a popularmethod
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 (eta*)
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 (eta*
approximate to 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 instability.