On the Quality of Velocity Interpolation Schemes for Marker-in-Cell Method and Staggered Grids

Abstract

The marker-in-cell method is generally considered a flexible and robust method to model the advection of heterogenous non-diffusive properties (i.e., rock type or composition) in geodynamic problems. In this method, Lagrangian points carrying compositional information are advected with the ambient velocity field on an Eulerian grid. However, velocity interpolation from grid points to marker locations is often performed without considering the divergence of the velocity field at the interpolated locations (i.e., non-conservative). Such interpolation schemes can induce non-physical clustering of markers when strong velocity gradients are present (Journal of Computational Physics 166:218–252, 2001) and this may, eventually, result in empty grid cells, a serious numerical violation of the marker-in-cell method. To remedy this at low computational costs, Jenny et al. (Journal of Computational Physics 166:218–252, 2001) and Meyer and Jenny (Proceedings in Applied Mathematics and Mechanics 4:466–467, 2004) proposed a simple, conservative velocity interpolation scheme for 2-D staggered grid, while Wang et al. (Geochemistry, Geophysics, Geosystems 16(6):2015–2023, 2015) extended the formulation to 3-D finite element methods. Here, we adapt this formulation for 3-D staggered grids (correction interpolation) and we report on the quality of various velocity interpolation methods for 2-D and 3-D staggered grids. We test the interpolation schemes in combination with different advection schemes on incompressible Stokes problems with strong velocity gradients, which are discretized using a finite difference method. Our results suggest that a conservative formulation reduces the dispersion and clustering of markers, minimizing the need of unphysical marker control in geodynamic models.

Appendices

Appendix 1: Integration Equations

Fig. 9

Interpolation sketch for staggered grid and $V_y$ velocity component

The general interpolation schemes used in this study are explained in more detail below for one velocity component in 1-D. In Fig. 9, three staggered grid points, $V_y^{i-1,j}, V_y^{i,j}, V_y^{i+1,j}$, are located at ($x_0, x_1, x_2$) points. $(U_{i-1}(x_c^{i-1}), U_i(x_c^{i})$) are grid nodes and $P(x_p)$ is a marker point. We illustrate here interpolation schemes to calculate velocities in the nodes ($U_i, U_{i-1}$) needed for the correction interpolation method. However, the equations below can also be used to calculate marker velocity directly from staggered points (i.e. direct interpolation) or other intermediate velocities steps. The velocity in the nodes, $U_i$, is calculated in the following way:

• (a) Linear interpolation

  $$\begin{aligned} U_i = (1-x_c'^i)(V_y^{i,j} + V_y^{i+1,j}), \end{aligned}$$

  (16)

  where $x_c'^i$ is the local coordinate of the node $U_i$ (i.e. $x_c'^i = (x_c^i-x_1)/\Delta x_i$). The same principle is applied to the other nodes and velocity components. The 3-D tri-linear interpolation scheme is represented by Eq. (5) in the main text.

• (b) Linear interpolation combined with a minmod limiter

  $$\begin{aligned} U_i = V_y^{i,j} - \frac{\displaystyle \Delta x}{\displaystyle 2} minmod\left( \frac{\displaystyle V_y^{i+1,j}-V_y^{i,j}}{\displaystyle \Delta x_i} , \frac{\displaystyle V_y^{i,j}-V_y^{i-1,j}}{\displaystyle \Delta x_{i-1}}\right) , \end{aligned}$$

  (17)
  $$\begin{aligned} U_{i-1} = V_y^{i,j} + \frac{\displaystyle \Delta x}{\displaystyle 2} minmod\left( \frac{\displaystyle V_y^{i+1,j}-V_y^{i,j}}{\displaystyle \Delta x_i} , \frac{\displaystyle V_y^{i,j}-V_y^{i-1,j}}{\displaystyle \Delta x_{i-1}}\right) , \end{aligned}$$

  (18)

  where the minmod slope limiter is defined as:

  $$\begin{aligned} minmod(A,B) = {\left\{ \begin{array}{ll} A, \quad {\text {if}}\; A\cdot B \ge 0& \quad {\text { and }} \quad |A| \le |B| \\ B, \quad {\text {if}}\; A\cdot B \ge 0&\quad \text { and }\quad |A| \ge |B| \\ 0, \quad {\text{if}}\; A\cdot B \le 0&\quad \\ \end{array}\right. } \end{aligned}$$

  (19)

• (c) Quadratic interpolation—calculates one quadratic curve, passing through all three staggered velocity points. As such, the velocity at the node or any point is calculated as:

  $$\begin{aligned} U_i = a x^2 + b x + c, \end{aligned}$$

  (20)

  where $x = x_c^i$ and the coefficients a, b, c are given by:

  $$\begin{aligned} a&= \frac{V_y^{i-1,j}}{(x_0-x_1)(x_0-x_2)}+\frac{V_y^{i,j}}{(x_1-x_0)(x_1-x_2)}+\frac{V_y^{i+1,j}}{(x_2-x_0)(x_2-x_1)}\end{aligned}$$

  (21)
  $$\begin{aligned} b&= \frac{V_y^{i-1,j}-V_y^{i,j}}{x_0-x_1} - a(x_0+x_1)\end{aligned}$$

  (22)
  $$\begin{aligned} c&= V_y^{i-1,j} - a x_0^2 - b x_0. \end{aligned}$$

  (23)

• (d) Spline quadratic interpolation—calculates two quadratic curves, one between $V_y^{i-1,j}$ and $V_y^{i,j}$ and another between $V_y^{i,j}$ and $V_y^{i+1,j}$, and the derivative is continuous at $V_y^{i,j}$. As such, the velocity at the node is calculated as:

  $$\begin{aligned} U_{i-1}&= a_1 x^2 + b_1 x +c_1,\quad \text{with}\;\; x = x_c^{i-1}\end{aligned}$$

  (24)
  $$\begin{aligned} U_i&= a_2 x^2 + b_2 x +c_2, \quad \text{with}\;\; x = x_c^i \end{aligned}$$

  (25)

  and where the coefficients $a_1, b_1, c_1, a_2, b_2, c_2$ are given by:

  $$\begin{aligned} a_1&= 0 \end{aligned}$$

  (26)
  $$\begin{aligned} b_1&= \frac{V_y^{i,j}-V_y^{i-1,j}}{x_1-x_0}\end{aligned}$$

  (27)
  $$\begin{aligned} c_1&= V_y^{i-1,j} - b_1 x_0 \end{aligned}$$

  (28)
  $$\begin{aligned} a_2&= \left( \frac{V_y^{i,j}-V_y^{i+1,j}}{x_1-x_2} - b_1\right) \frac{1}{x_2-x_1}\end{aligned}$$

  (29)
  $$\begin{aligned} b_2&= b_1 - 2 a_2 x_1 \end{aligned}$$

  (30)
  $$\begin{aligned} c_2&= V_y^{i,j} - b_2 x_1 - a_2 x_1^2 \end{aligned}$$

  (31)

Appendix 2: Extended Results

In this section, we present the full results for the experiments performed in Sect 3. Figure 10 shows results for the SolCx benchmark for different resolutions, initial marker distribution and viscosity contrasts. Figure 10a, b shows that marker artifacts occur for all interpolation schemes, with the exception of LinP and CorrMinmod, and are independent of initial marker distribution and resolution. Further tests (Fig. 10c) suggest that artifacts become visible and problematic for $\Delta \eta > 10^3$. Figure 11 shows results for all the interpolation methods for the Rayleigh–Taylor instability test.

Fig. 10

Extended results for the SolCx benchmark

Fig. 11

Extended results for the Rayleigh–Taylor instability test