# Superconvergent second order Cartesian method for solving free boundary problem for invadopodia formation

Olivier Gallinato, Clair Poignard

► To cite this version:

Olivier Gallinato, Clair Poignard. Superconvergent second order Cartesian method for solving free boundary problem for invadopodia formation. Journal of Computational Physics, 2017, 339, pp.412 - 431. 10.1016/j.jcp.2017.03.010. hal-01483484

HAL Id: hal-01483484

https://inria.hal.science/hal-01483484

Submitted on 6 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Superconvergent Second Order Cartesian Method for Solving Free
Boundary Problem for Invadopodia Formation

Olivier Gallinato<sup>a</sup>, Clair Poignard<sup>a,1</sup>

<sup>a</sup>Team MONC, INRIA Bordeaux-Sud-Ouest, Institut de Mathématiques de Bordeaux, CNRS UMR 5251 & Université de Bordeaux, 351 cours de la Libération, 33405 Talence Cedex, France.

**Abstract**

In this paper, we present a superconvergent second order Cartesian method to solve a free boundary problem with two harmonic phases coupled through the moving interface. The model recently proposed by the authors and colleagues describes the formation of cell protrusions. The moving interface is described by a level set function and is advected at the velocity given by the gradient of the inner phase. The finite differences method proposed in this paper consists of a new stabilized ghost fluid method and second order discretizations for the Laplace operator with the boundary conditions (Dirichlet, Neumann or Robin conditions). Interestingly, the method to solve the harmonic subproblems is superconvergent on two levels, in the sense that the first and second order derivatives of the numerical solutions are obtained with the second order of accuracy, similarly to the solution itself. We exhibit numerical criteria on the data accuracy to get such properties and numerical simulations corroborate these criteria. In addition to these properties, we propose an appropriate extension of the velocity of the level-set to avoid any loss of consistency, and to obtain the second order of accuracy of the complete free boundary problem. Interestingly, we highlight the transmission of the superconvergent properties for the static subproblems and their preservation by the dynamical scheme. Our method is also well suited for quasistatic Hele-Shaw-like or Muskat-like problems.

**Keywords:** Finite differences on Cartesian grids, Superconvergence, Interface conditions, Free boundary problem

2000 MSC: 65M06, 65M12, 92C37

**1. Introduction**

Early-stage carcinoma are mostly confined to the epithelium, which is separated from the underlying tissue by a basement membrane composed of dense fibers of extracellular matrix (ECM).

<footnote>Email addresses: olivier.gallinato@inria.fr (Olivier Gallinato), clair-poignard@inria.fr (Clair Poignard)</footnote>

<footnote n="1">¹Corresponding author</footnote>
In order to cross this tight barrier, metastatic cells use a complex internal machinery, named *invadopodia*, which lies on the actin polymerization and that leads to the formation of proteolytic, protrusive and very localized subcellular structures. *Invadopodia* are elongated shapes, which are formed during cell invasion and mesenchymal migration. This phenomenon is the crucial and initiating point in the metastatic process, which is the major cause of death from cancer. The authors and colleagues proposed recently a free boundary problem to model invadopodia and more generally for cell protrusion formation [16]. In this paper, we present a superconvergent second order method on Cartesian grid to solve this quite complex free boundary value problem. Our finite difference method is based on the derivation of appropriate superconvergent schemes: wide-stencils are proposed to reach specific superconvergence properties for the solutions of the Poisson problem. As a result, the global method leads to the second-order accuracy of the moving interface, its normal vector and even its curvature. The accuracy on the curvature will be of major interest in further study for the use of interface regularization techniques in order to model subsequent phenomena involved in cell migration, as myosin-dependant protrusion retraction, for instance.

The purpose of the study is to simulate accurately cell protrusion formation, however the domain of applications can be extended to free-boundary problems arising from physics or biology such as Hele-Shaw, Muskat or quasistatic two-phase Stefan like problems.

### 1.1. Free-boundary problem for invadopodium formation

The invadopodium process relies on a coupled dynamics between the outer and the inner of the cell. Let us present briefly the new model of the phenomenon detailed in [16]. Specific enzymes (MT1-MMPs) produced by the cell membrane degrade the extracellular matrix (ECM), producing ligands that diffuse and bind to membrane receptors. In response, cell generates a signal which diffuses inside the cell and, which leads to actin polymerization: rigid filaments are polymerized, oriented towards the location of the detected ligand. The force exerted by the filaments on the membrane generates a protrusion which grows at the velocity of the filament polymerization. The scheme of the process is given by Figure 1(a), while the geometrical framework is detailed in Figure 1(b). At any time $t$, the cell membrane is parameterized by the map $\gamma(t, \cdot)$ defined on the torus $T = \mathbb{R}/2\pi\mathbb{R}$:

$$ \Gamma_t = \{\gamma(t, \theta), \theta \in T\}. $$

The cell cytoplasm $O_t^i$ is the domain enclosed by $\Gamma_t$ and the ECM is the outer domain

$$ O_t^e = \Omega \setminus O_t^i. $$

Assume the flux of MT1-MMP enzymes $g(t, \cdot)$ be given at any time on the cell membrane. It generates a flux of the degraded matrix (called ligands and denoted by $c^*$) on the cell boundary, and these ligands diffuse in the extracellular medium as described by equations (1a)-(1b). When bound to the cell membrane, the ligands generate a signal $\sigma$, which diffuses inside the cell, as
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Schematic diagram of the molecular interactions involved in our model and geometrical settings. The cell $O_t^i$ is imbedded in the bath $O_t^e$. The whole domain $\Omega$ does not depend on the time variable. It is defined by $\Omega = O_t^e \cup \overline{O_t^i}$</figcaption></figure>

accounted for in equations (1c)–(1d). The cell membrane motion is described by equation (1e).

Degradation of the ECM:

$$ \Delta c^* = 0, \quad x \in O_t^e, \qquad (1a) $$

$$ c^*|_{\partial\Omega} = 0, \quad -\partial_n c^*|_{\Gamma_t} = g|_{\Gamma_t}. \qquad (1b) $$

Generation of the inner signal for actin polymerization:

$$ \Delta\sigma = 0, \quad t \in [0, T], \quad x \in O_t^i, \qquad (1c) $$

$$ \sigma|_{\Gamma_t} = c^*|_{\Gamma_t}. \qquad (1d) $$

Motion of the cell membrane:

$$ \partial_t \gamma(t, \theta) = \nabla \sigma(\gamma(t, \theta)), \quad \theta \in \mathbb{T}, \quad \text{and} \quad \Gamma_t = \{\gamma(t, \theta), \theta \in \mathbb{T}\}. \qquad (1e) $$

The interested reader will refer to [16] for further details about biological phenomena, and modeling hypotheses. The theoretical analysis of the free-boundary problem is also performed in this article. In particular, the well-posedness of the free-boundary problem in Sobolev spaces is precisely proven, provided strictly positive boundary data $g|_{\Gamma_t}$. The proof is based on the explicit characterization of Dirichlet-to-Neumann maps thanks to complex analysis tools. Then, appropriate change of unknowns and the quasilinearization of the problem leads to a parabolic-type equation, which prevents the loss of regularity and thus provides the well-posedness.

## 1.2. Eulerian approach for the numerical formulation

In order to avoid any remeshing of the domain, specific finite difference methods on Cartesian grids regarding the Eulerian formulation of problem (1) is proposed. Such a formulation consists in detecting the cell membrane by a level-set function $\psi(t, x)$ defined in the entire domain $\Omega$, and
such that

$$ \Gamma_t = \{x : \psi(t, x) = 0\}, \quad O_t^i = \{x : \psi(t, x) < 0\}, \quad O_t^e = \{x : \psi(t, x) > 0\}. \qquad (2) $$

Equation (1e) is then changed into a transport equation for the level-set at the velocity $\mathbf{v}$, which derives from the gradient of $\sigma$ as stated in (1e). Here arises a first numerical difficulty, since the velocity of the Lagrangian formalism is only defined on the interface, i.e. $\mathbf{v}|_{\Gamma_t} = \nabla\sigma|_{\Gamma_t}$. An extension of the velocity to the whole domain is thus required in order to solve the transport equation (3a). We thus replace equation (1e) by

$$ \partial_t \psi + \mathbf{v} \cdot \nabla \psi = 0, \quad t \ge 0, \quad x \in \Omega, \qquad (3a) $$

where the velocity $\mathbf{v}$ is extended to $\Omega$ from the data $\nabla\sigma|_{\Gamma_t}$ as follows:

$$ (\nabla\psi \cdot \nabla)\mathbf{v} = 0, \quad \forall x \in \Omega, \quad \text{with} \quad \mathbf{v}|_{\{\psi(t,x)=0\}} = \nabla\sigma|_{\{\psi(t,x)=0\}}. \qquad (3b) $$

**Remark 1** (Comment on the velocity extension). The solution to problem (1a)-(1b)-(1c)-(1d)-(3) depends on the choice of the extension, however only the level-set function would change, but the location of $\Gamma_t$ and the quantities $\sigma$, and $c^*$ would be the same. Therefore the choice of the extension operator is only guided to provide appropriate numerical properties.

As explained in the following, the extension (3b) will provide the good properties to reach second order of accuracy on the interface, its normal and even its curvature. However it is obvious that at any time step, shocks are generated by this extension, and merge with ridges of the level set function, initialized as the signed distance function. If close to the interface, shocks may dramatically decrease the accuracy of the numerical solution. To avoid such drawbacks it is necessary to refine the mesh, for instance thanks to local adaptive meshing such as quadtree or octree. This is far from the scope of the present paper. Therefore, throughout the paper we consider that at any time of simulation the interface $\Gamma_t$ is far from a shock, and the neighborhood of the interface, in which the velocity is continuously extended can therefore be defined as any area without such ridge.

## 1.3. Motivation

Solving accurately the above partial differential equations (PDEs) is crucial in order to provide relevant biological interpretations of the model, however it is very challenging. One of the main difficulty lies in the fact that the velocity of the cell membrane is given as the gradient of the inner signal, which provides a priori a loss of accuracy<fnref n="2" />: usual second order schemes could lead to a first order accuracy of the gradient, and thus the level-set would be not enough accurate.

<footnote n="2">²Such a possible loss of accuracy is similar to the a priori loss of regularity of the continuous solution, as described in [16].</footnote>
The loss of consistency could then be drastically increased by the coupling and the transmission of errors, leading to inconsistent scheme, and thus irrelevant simulations. In what follows, to overcome these difficulties we propose a numerical method based on the use of superconvergence, as defined in Definition 2, and on its transmission from one sub-problem to another.

**Definition 2 (Superconvergent solution).** A numerical approximation $f$ on a grid of a function $f$ of the continuous domain is said superconvergent if, in addition to approximating $f$ with a given order of accuracy, the numerical gradient $\nabla^h f$ of $f$ approximates the gradient $\nabla f$ of $f$ with the same order of accuracy.

If in addition to the gradient, the numerical second order derivatives are also obtained with the same accuracy, then the numerical approximation $f$ is said superconvergent on two levels (to-wards $f$).

In this paper, we numerically highlight that our method makes it possible to compute the solution of the Poisson problem with two levels of superconvergence, under specific conditions. Among them, the accuracy on the boundary data is crucial. Interestingly, our superconvergent method does not require exact data. This is one of the main novelty of the study: specific requirements on the data (the boundary datum (g), possible volumic source terms (f) in the static problems (1c)-(1a) –for invadopodium model $f$ vanishes– and the level set ($\psi$)) are exhibited in Section 3.2 to obtain superconvergent solution to each static subproblem, and we show that these requirements are compatible with the overall dynamical coupling. Indeed, as data do not need to be exact, the transmission of the superconvergence properties from one sub-problem to the other is ensured.

The superconvergent method is based on a stabilized version of the Ghost Fluid Method, introduced by Fedkiw et al. [12]. Depending on the extrapolation level of the ghost values (linear, quadratic or cubic extrapolations can be considered), we obtained different results of superconvergence for the numerical solutions of the Poisson problem and of the overall problem of protrusion formation. In [16], linear extrapolations led to a specific first order method for the Neumann problem (1a)-(1b), and a second order method for the Dirichlet problem (1c)-(1d). The authors numerically showed that the resulting quantities $\sigma$, $c^*$, $\psi$ and even the normal to the interface are first order accurate and suggested the involvement of underlying superconvergence properties related to the methods of resolution of the Poisson equations. Actually, the consistency of the solutions partially results from the first order superconvergence of the solution of the Neumann problem, which is transmitted to the solution of the inner problem thanks to the Dirichlet condition.

Hence, this result encouraged the authors to continue the study of the superconvergence properties and their transmission, to obtain second order numerical approximations of $c^*$, $\sigma$ and $\psi$. The main idea is to use higher extrapolations for the computation of the ghost values. It is worth noting that second order accurate methods to solve separately each problem (1a)-(1b), (1c)-(1d) (for a given location of the level set function $\psi$), and (3) (for a given velocity at the interface) do not necessarily lead to an overall second order accuracy of the coupled problem (1a)-(1b)-(1c)-(1d)-
(3). Actually, the second order accuracy of $c^*$ requires a second order accuracy on the Neumann boundary condition (see [10] for instance), and the second order accuracy of $\sigma$ requires a second order accuracy on the Dirichlet boundary condition. Thus, the interface and its normal should be second order accurate.

Another specificity of our method lies on the velocity extension (3b), which leads to the super-convergence of the level set function on one level. This surprising behavior makes it possible to obtain the second order accuracy of the level set function and the normal vector only from a second order accurate velocity, without any additional superconvergent property. This only requires super-convergence property on one level for $\sigma$. This result can be obtained from quadratic extrapolations of the ghost values for each static sub-problem, which give numerical schemes of Shortley-Weller type [32], for the discretization of the Poisson problems.

Then, the spontaneous following question relies upon the cubically extrapolated ghost values. Obviously, the stencils are wider. However, in addition to $\sigma$, $c^*$, $\psi$ and their first derivatives, their second order derivatives will also be second order accurate in maximum norm, exhibiting a superconvergence property on two levels. This is particularly the case for the interface curvature. This is of major interest for further study about cell migration. Indeed, the protrusion retraction could be accounted for by using a specific retraction velocity $v_{\text{ret}}$ and techniques of interface regularization, which depend on surface tension. For instance, we could consider a well-known interface advection problem with curvature-dependent speed, introduced by Sethian and Osher in [29], and given by

$$ \partial_t \psi + v_{\text{ret}} \cdot \nabla \psi = F(\kappa), \quad t \in [T, T'], \quad x \in \Omega, $$

where $F$ is a given nonlinear function of the curvature $\kappa$. Nonlinearities of $F$ require the second order accuracy on $\kappa$ in order to keep the appropriate properties on $\psi$. This aspect is far beyond the topic of this paper and will be addressed in further studies. However, this justifies the presentation of cubic type schemes in what follows.

The main motivation of the paper is therefore to propose a specific finite difference method on Cartesian grid to compute every unknowns and first and second derivatives of the problem of protrusion formation with a second order accuracy, thanks to a stabilized Ghost Fluid Method and cubic extrapolations of the ghost values for the Poisson problems. It is understood that the reader, if uninterested in the second order accuracy of the curvature for further regularization approach, can easily replace the cubic extrapolations by quadratic extrapolations and keep second order accurate solutions of the problem, except the curvature.

Hence, the use of superconvergence properties can be adjusted to other problems. Actually, superconvergence may be a powerfull property to solve many free boundary problems, for which the interface velocity is implicitly given by one component of the solution. This may be particularly interesting to solve static 2-phase Stefan-type problems [17, 7, 31] or Hele-Shaw-like problems [19].
1.4. Superconvergence: state of the art

To our best knowledge, superconvergence on two levels and superconvergence from perturbed
data have never been mentioned in literature. Several studies have addressed the one order of
superconvergence, especially for finite element [40, 23, 2] and finite volume type methods [8, 3, 38].
Regarding finite difference methods on polygonal domains, Ferreira and Grigorieff have shown
in [13] one level of superconvergence at the second order of accuracy for general elliptic operators
while Li et al. in [26, 24, 25] specifically studied the Shortley-Weller scheme. In [27], Matsunaga
and Yamamoto proved that the Shortley-Weller scheme provides a third order accuracy near the
interface, which can be seen as an essential aspect of the superconvergence property. The Shortley-
Weller approximation, introduced in [32] in 1938, has been widely studied [4, 27, 35, 26, 36] and
is well-known for its superconvergence property, without complete proof until very recently. More
precisely, in [37], Yoon and Min proved the superconvergence in L²-norm at the second order of
accuracy thanks to a discrete version of the divergence theorem, and Weynans, in [34], proposed
a different proof for the maximum norm, involving Green's functions and analysis tools strongly
inspired from [4, 9]. From a numerical point of view, Gibou and colleagues observed in [28] and [6]
the superconvergence of the solution to Poisson equation with variable coefficients on non-graded
adaptive Cartesian grids.

It is worth noting that all these studies consider only exact data. Moreover, they are limited
to one level of superconvergence and to Dirichlet boundary conditions. In addition to these works,
Jomaa et al. achieved in [20] the second order accuracy of the solution to the Robin problem
using the Shortley-Weller method with a quadratic treatment of the boundary condition. For the
Neumann condition, Kreiss and more recently Khatri [22, 21] also pointed out the second order
accuracy of the solution thanks to quadratic or cubic extrapolations for the Neumann condition.
However, none of these articles refer to superconvergence properties of their schemes.

Superconvergent approaches are well-suited for solving models in which one unknown derives
from the solution of a Poisson problem. In these approaches, the specific superconvergent properties
of the numerical solution of the Poisson problem, which can be seen as a discrete version of the
elliptic regularity properties like in [15], prevent the loss of consistency at each derivation. However,
other approaches should be considered. In particular, Bruno and Hoch propose in [5] a Chebyshev
polynomial technique to postprocess the derivative computation and get a limited order-of-accuracy
deterioration compared to the solution. This could be an interesting alternative to explore for
solving the invadopodia model.

Finally, we should mention that in some cases, after a long time, the interface can lose in
regularity or develop pathological configurations. In each situation, the method should be adjusted
in order to avoid instabilities, as this is addressed in [15]. In other cases, the advection governed
by (3a), could lead to instabilities in our numerical approach if the interface develop a kink for
instance. Then, alternative approaches should be considered, as proposed by Friedman in [14] for
Stefan problem.
## 1.5. Outline of the paper

As stated above, the main goal of this paper is to solve Problem (1) with the second order accuracy, using a superconvergent second order accurate Cartesian method for Poisson-type equations.

The schemes used to reach the superconvergent properties are presented in Section 2. There are based on centered schemes and the use of the ghost fluid method [12]. The keypoint lies in the appropriate extrapolations (cubic for the static subproblems, quadratic for the transport of the level set) for the ghost values to obtain the superconvergence. In particular we propose new stencils, partially introduced in [16], which stabilize the standard stencils such as these of Gibou et al. [18] or Cisternino and Weynans [10] for the grid points close to the interface and for the points on the interface.

Section 3 is devoted to a step-by-step study of the overall method: we first highlight the superconvergent properties of the schemes. In particular we show that the use of cubic extrapolations in the static subproblems leads to a solution with two levels of superconvergence. Then, in subsection 3.2.3, we show that superconvergent properties are preserved even with perturbed data. More precisely, we highlight that our method ensures the preservation of the superconvergence properties to the static coupled problem (1a)-(1b)-(1c)-(1d). We finally show that the transport (3), solved thanks to a RK2-LUD scheme [33], do not degrade the two levels of superconvergence provided that the extended velocity is superconvergent on one level.

Again, we emphasize that linear and quadratic extrapolations can similarly be used instead of cubic extrapolations in the discretization of the static part of the model (Poisson problems), resulting in first and second order accuracies of the respective overall methods. As said above, the quadratic approach does not result in the second order accuracy of the curvature. Note that the corresponding solutions of the so-discretized Poisson problem have specific superconvergence properties, which are summarized in the final Table 9 without details. The interested reader will refer to [15] for more informations about those alternative finite-difference superconvergent methods to solve the Poisson problem on Cartesian grid.

## 1.6. Illustration of the main results

The numerical method proposed in this paper makes it possible to solve the invadopodia model given by equations (1a)-(1b)-(1c)-(1d)-(3) with a second order accuracy for the unknowns $\sigma, c^*$, $\psi$ and in addition for the normal vector field $\mathbf{n}$ and the curvature $\kappa$ defined from the level-set $\psi$ by

$$ \mathbf{n} = \frac{\nabla \psi}{|\nabla \psi|}, \quad \kappa = \nabla \cdot (\mathbf{n}). $$

This is summarized in the following numerical simulation, where $g$ is given by

$$ g = 0.1 [2 + \cos(3\pi(x+y)) \cos(\pi(x+0.3))], $$

and $\Gamma|_{t=0}$ is the circle centered in $\mathbf{O}$ of radius 0.25. Figure 2 provides the changes of the cell morphology given by the model, that is the shape of the cell membrane $\Gamma_t$, at different time points.
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Numerical simulation of a change in cell morphology by solving the problem of invadopodia given by (1a)–(1b)–(1c)–(1d)–(3).</figcaption></figure>

The interface velocity, before extension, is computed as the gradient of a harmonic function. As expected, Figure 3 shows a satisfying conservation in the volume.

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Volume conservation.</figcaption></figure>

Table 1 illustrates the accuracy of the numerical method. The quantities $\sigma$, $c^*$, $\psi$, $n$, $\kappa$ are second order accurate in $L^\infty$ norm, by comparison with a numerical reference solution, computed on a fine mesh (500 x 500 points), at $T = 1$ (after 659 time steps).

## 2. Discretizations of the differential operators for the dynamical coupled problem

This section is devoted to the discretization of the differential operators. Throughout the paper, we consider a Cartesian grid, denoted by $G^h$, with the same grid spacing in $x$ and $y$, denoted by $h$. The set of the grid points is denoted by $\mathcal{P}^h$. The continuous domain $\omega$ generically stands either for $O_t^e$, $O_t^i$ or $\Omega$. Let first set the notations in order to present clearly the discretizations.
<table><caption>**Table 1:** Numerical results of accuracy for a simulation of protrusion formation at the final time $T = 1$.</caption><thead><tr><th rowspan="2">Number of points</th><th rowspan="2">Number of time steps</th><th colspan="2">c<sup>*</sup></th><th colspan="2">σ</th><th colspan="2">v</th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td>64</td><td>86</td><td>4.830 × 10<sup>-5</sup></td><td>-</td><td>3.733 × 10<sup>-5</sup></td><td>-</td><td>7.769 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>96</td><td>128</td><td>1.606 × 10<sup>-5</sup></td><td>2.72</td><td>7.129 × 10<sup>-6</sup></td><td>4.08</td><td>2.297 × 10<sup>-4</sup></td><td>3.01</td></tr><tr><td>144</td><td>190</td><td>4.883 × 10<sup>-6</sup></td><td>2.83</td><td>3.237 × 10<sup>-6</sup></td><td>3.02</td><td>1.167 × 10<sup>-4</sup></td><td>2.34</td></tr><tr><td>216</td><td>285</td><td>1.793 × 10<sup>-6</sup></td><td>2.71</td><td>2.384 × 10<sup>-6</sup></td><td>2.26</td><td>6.106 × 10<sup>-5</sup></td><td>2.09</td></tr><tr><td>324</td><td>427</td><td>1.270 × 10<sup>-6</sup></td><td>2.24</td><td>3.032 × 10<sup>-7</sup></td><td>2.97</td><td>2.478 × 10<sup>-5</sup></td><td>2.12</td></tr><tr><th colspan="2"></th><th colspan="2">ψ</th><th colspan="2">n</th><th colspan="2">κ</th></tr><tr><td>64</td><td>86</td><td>3.620 × 10<sup>-4</sup></td><td>-</td><td>8.617 × 10<sup>-3</sup></td><td>-</td><td>5.089 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>96</td><td>128</td><td>1.627 × 10<sup>-4</sup></td><td>1.97</td><td>5.049 × 10<sup>-3</sup></td><td>1.32</td><td>2.365 × 10<sup>-1</sup></td><td>1.89</td></tr><tr><td>144</td><td>190</td><td>8.066 × 10<sup>-5</sup></td><td>1.85</td><td>2.383 × 10<sup>-3</sup></td><td>1.59</td><td>1.381 × 10<sup>-1</sup></td><td>1.61</td></tr><tr><td>216</td><td>285</td><td>3.368 × 10<sup>-5</sup></td><td>1.95</td><td>9.596 × 10<sup>-4</sup></td><td>1.80</td><td>6.337 × 10<sup>-2</sup></td><td>1.71</td></tr><tr><td>324</td><td>427</td><td>1.069 × 10<sup>-5</sup></td><td>2.17</td><td>3.270 × 10<sup>-4</sup></td><td>2.02</td><td>2.050 × 10<sup>-2</sup></td><td>1.98</td></tr></tbody></table>

## 2.1. Notations and definitions of the discrete quantities

We define the numerical level-set function $\psi$ as a perturbation of the restriction of the exact level-set $\psi$ to the grid points $\mathcal{P}^h$. Then, we denote by $\Gamma_t^h$, $\mathcal{O}_t^{i,h}$ and $\mathcal{O}_t^{e,h}$ the numerical interface and the inner and outer numerical domains, respectively:

$$ \Gamma_t^h = \{x \in \mathcal{G}^h : \psi(t, x) = 0\}, \quad \mathcal{O}_t^{i,h} = \{x \in \mathcal{P}^h : \psi(t, x) < 0\}, \quad \mathcal{O}_t^{e,h} = \{x \in \mathcal{P}^h : \psi(t, x) > 0\}. $$

The discrete domain $\omega^h$ generically stands either for $\mathcal{O}_t^{e,h}$, $\mathcal{O}_t^{i,h}$ or $\mathcal{P}^h$. Note that these definitions are the discrete versions of (2).

### 2.1.1. Definitions of the grid, the neighboring and the interface points

Grid points of $\omega^h$ are simply denoted $x_{i,j}$. A neighboring point is a grid point whose one of its closest neighborhood (at least) does not belong to $\omega^h$. Other grid points are called inner points.

For $\theta \in [0, 1]$, the point on the grid located at the distance $\pm h\theta$ in the x-direction (resp. y-direction) of the grid point $x_{i,j}$ is denoted by $x_{i\pm\theta, j}$ (resp. $x_{i, j\pm\theta}$). If the point $x_{i\pm\theta, j}$ (resp. $x_{i, j\pm\theta}$) belongs to the intersection between the grid and the interface $\Gamma^h$, it is named interface point and may be denoted by $x_{i\pm\theta, j}$ (resp. $x_{i, j\pm\theta}$) if necessary.

For any function $u$ defined in the continuous domain $\omega$, we denote by $u_{i\pm\theta_x, j\pm\theta_y}$ its value at the point $x_{i\pm\theta_x, j\pm\theta_y}$. The value of $u$ at the interface points $x_{i\pm\theta_x, j\pm\theta_y}$ is denoted by $u_{i\pm\theta_x, j\pm\theta_y}$.

Since the discretization is performed direction by direction, only the x-direction is presented in the following, and the subscript $j$ is sometimes omitted to simplify notations. Hence $u_i$ stands for $u_{i,j}$, where $j$ is fixed at the same value.
### 2.1.2. Quadratic and cubic extrapolation operators

As said above, the numerical schemes are based on appropriate extrapolations. More precisely, the quadratic and cubic extrapolations, respectively denoted by $\mathbf{Q}$ and $\mathbf{C}$ are defined as

$$ \mathbf{Q}(\mathbf{Z}; \boldsymbol{\theta}) = \sum_{i=1}^{3} a_i^{\mathbf{Q}}(\boldsymbol{\theta}) \mathbf{Z}_i, \quad \forall (\mathbf{Z}; \boldsymbol{\theta}) \in \mathbb{R}^3 \times [-3, 3]^3, \qquad (4a) $$

$$ \mathbf{C}(\mathbf{Z}; \boldsymbol{\theta}) = \sum_{i=1}^{4} a_i^{\mathbf{C}}(\boldsymbol{\theta}) \mathbf{Z}_i, \quad \forall (\mathbf{Z}; \boldsymbol{\theta}) \in \mathbb{R}^4 \times [-4, 4]^4, \qquad (4b) $$

where the vector fields $a^{\mathbf{L}}: [-n_{\mathbf{L}}, n_{\mathbf{L}}]^{n_{\mathbf{L}}} \to \mathbb{R}^{n_{\mathbf{L}}}$ are given by

$$ a_i^{\mathbf{L}}(\boldsymbol{\theta}) = \frac{\prod_{\substack{k=1,\cdots,n_{\mathbf{L}} \\ k \neq i}} \theta_k}{\prod_{\substack{k=1,\cdots,n_{\mathbf{L}} \\ k \neq i}} (\theta_k - \theta_i)}, \quad \text{for } \mathbf{L} = \mathbf{Q}, \mathbf{C}, \qquad (4c) $$

the integer $n_{\mathbf{L}}$ being equal to 3 (resp. 4) for $\mathbf{L} = \mathbf{Q}$ (resp. $\mathbf{L} = \mathbf{C}$).

**Remark 3** (About quadratic and cubic extrapolations). One can wonder why cubic and quadratic extrapolations have to be introduced. First of all, one can notice that in order to get a truncation error on the Laplace operator with a second order of accuracy, it is necessary to compute the ghost values with cubic extrapolation. Such cubic extrapolation adds an error of order $h^3$ to the gradient of the solution, making it possible, if the numerical scheme is superconvergent on two levels, to derive numerically two times the numerical solution, keeping the second order of accuracy of the first and second order derivatives. Quadratic extrapolation is then sufficient for the second numerical derivation to keep this second order of accuracy. Moreover, as noted below, quadratic extrapolation is necessary to prevent numerical instabilities (which appear with cubic extrapolations) for the velocity extension and the transport of the level set.

### 2.1.3. Localization of the level set

The distance $h\theta_x$ of the interface point $x_{i-\theta_x}^\Gamma$ to the neighboring point $x_{i,j}$, if not exactly known, has to be computed. In order to add a second order error at worst on the second order of derivatives of the level-set, we use a fourth order approximation. For such purpose, the level set function $\psi$ is approximated with a third order interpolation polynomial $P$ in the Lagrange form. Denoting by $\alpha$ the normalized distance between $x_i$ and any point $x_{i-\alpha}$ of the x-axis, the polynomial $P$ can be written as a function of $\alpha$:

$$ P(\alpha) = \frac{1}{6}(\psi_{i-2} - 3\psi_{i-1} + 3\psi_i - \psi_{i+1})\alpha^3 + \frac{1}{2}(\psi_{i-1} - 2\psi_i + \psi_{i+1})\alpha^2 + \frac{1}{6}(-\psi_{i-2} + 6\psi_{i-1} - 3\psi_i - 2\psi_{i+1})\alpha + \psi_i. \qquad (5) $$

Solving $P(\alpha) = 0$ thanks to Newton's method provides the desired approximation of $\theta_x$.
## 2.2. Discrete operators to solve the static subproblems (1a)-(1b) and (1c)-(1d)

### 2.2.1. Discrete Laplacian

The second order derivative $\partial_{xx}^h u_i$ of a function $u$ at the grid point $x_i$ is defined at the point thanks to the standard five-point stencil:

$$ \partial_{xx}^h u_i = \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2}. \qquad (6) $$

At any inner grid point, all the neighbors are still in $\omega^h$ and thus the above formula is well defined. However if $x_i$ is a neighboring point, then for instance $x_{i-1} \notin \omega^h$ and thus $u_{i-1}$ is not defined: it is called a ghost value -denoted by $u_{i-1}^G$ for clarity-and it has to be extrapolated from the values at the interface and at the inner points (see Figure 4).

<figure><img src="image_5.png" /><figcaption>Figure 4: Standard and the stabilized methods for the ghost value computation.</figcaption></figure>

The standard cubic Ghost fluid approach [12] consists in using $u_{i-\theta}^G$ and the inner points ($u_{i+k}$) for $k = 0, 1, 2$ (see Figure 4(a))

$$ u_{i-1}^G = C(u_{i-\theta_x}^G, u_i, u_{i+1}, u_{i+2}; \theta), \quad \text{with} \quad \theta = [1 - \theta_x, 1, 2, 3]^T. $$

The drawback is that the coefficients of the extrapolation involve terms<fnref n="3" /> in $1/\theta_x$, which blows up for $\theta_x \to 0$, leading to instabilities as mentioned in [18]. This may affect the second order of accuracy of the solution and also prevent from any superconvergent properties. We thus propose the following extrapolation of $u_{i-1}^G$, in order to stabilize the discrete Laplacian operator:

$$ u_{i-1}^G = C(u_{i-\theta_x}^G, u_{i+1-\theta_x}, u_{i+2-\theta_x}, u_{i+3-\theta_x}; \theta_k)_{k=1, \dots, 4}, \quad \text{with} \quad \theta_k = (k - \theta_x). \qquad (7) $$

The new second order derivative at the neighboring point $x_i$ reads

$$ \partial_{xx}^C u_i = \frac{1}{h^2} \left[ u_{i+1} - 2u_i + \frac{(2-\theta_x)(3-\theta_x)(4-\theta_x)}{6} u_{i-\theta_x}^G - \frac{(1-\theta_x)(3-\theta_x)(4-\theta_x)}{2} u_{i+1-\theta_x} \right. \\ \left. + \frac{(1-\theta_x)(2-\theta_x)(4-\theta_x)}{2} u_{i+2-\theta_x} - \frac{(1-\theta_x)(2-\theta_x)(3-\theta_x)}{6} u_{i+3-\theta_x} \right], \qquad (8) $$

$$ \partial_{xx}^C u_i = \frac{6}{\theta_x(1+\theta_x)(2+\theta_x)h^2} u_{i-\theta_x}^G - \frac{3-\theta_x}{\theta_x h^2} u_i + \frac{2(2-\theta_x)}{(1+\theta_x)h^2} u_{i+1} - \frac{1-\theta_x}{(2+\theta_x)h^2} u_{i+2}. $$

<footnote n="3">³Simple calculations left to the reader lead to the formula for the standard Ghost fluid method with cubic extrapolation:</footnote>
Then, as the values $u_{i+k-\theta_x}$ are not at the grid points, they are extrapolated at the order 4 from the values at the grid points as follows:

$$u_{i+k-\theta_x, j} = C(u_{i,j}, u_{i+1,j}, u_{i+2,j}, u_{i+3,j}; (\theta_x - k + \ell)_{\ell=0,\dots,3}), \quad k \in \{1, \dots, 3\}, \quad (9)$$

and the discretization is finally of the form:

$$\partial_{xx}^C u_i = \frac{1}{h^2} [\alpha(\theta_x) u_i + \beta(\theta_x) u_{i+1} + \gamma(\theta_x) u_{i+2} + \delta(\theta_x) u_{i+3} + \eta(\theta_x) u_{i-\theta_x}^r],$$

Note that even though one more point is necessary compared to the standard method, the new discretization is thus stabilized, since the coefficients do not blow up for $\theta_x \to 0$.

**Remark 4** (Invertibility of the resulting discretization matrix in the 1D Dirichlet case). In contrast to the non-stabilized method, the discretization matrix is not diagonally dominant for a range of values of $\theta_x$. Nevertheless, it is invertible. The main arguments are as follows. Consider the matrices

$$A = (a_{ij})_{1 \le i, j \le n}, \quad a_{ij} = \begin{cases} -2 & \text{if } j = i, \\ 1 & \text{if } j = i \pm 1, \\ 0 & \text{otherwise,} \end{cases}$$

and $B = (b_{ij})_{1 \le i, j \le n}$, with $[b_{11}, b_{12}, b_{13}, b_{14}] = [\alpha(\theta_x) + 2, \beta(\theta_x) - 1, \gamma(\theta_x), \delta(\theta_x)]$, and the other coefficients $b_{ij}$ are 0.

The discretization matrix is therefore of the form $M = A + B = A(I_n + A^{-1}B)$. The matrix A is monotone and it can be easily checked by induction that its inverse matrix is given by

$$A^{-1} = (\bar{a}_{ij})_{1 \le i, j \le n}, \quad \bar{a}_{ij} = \begin{cases} -\frac{i(n-j+1)}{n+1} & \text{if } j \ge i, \\ -\frac{j(n-i+1)}{n+1} & \text{otherwise.} \end{cases}$$

The matrix $I_n + A^{-1}B$ is then a block matrix of the form $\left(\frac{U}{V} \mid \frac{O}{I_n}\right)$. It follows that

$$\det(I_n + A^{-1}B) = \det(U) = \frac{1}{6} [\theta_x^4 + (n-9)\theta_x^3 + (-9n+26)\theta_x^2 + (26n-24)\theta_x - 24n].$$

Deriving with respect to $n$, it is direct that $\det(U)$ is nonzero for any $n \ge 1$ and $\theta_x \in [0, 1]$.

We are confident on the fact that further study may lead to the same conclusion for 2D and Neumann problems.
### 2.2.2. Normal vector to the isolines

The normal vector $\mathbf{n}$ is given thanks to the level-set $\psi$ as

$$ \mathbf{n} = \frac{\nabla \psi}{|\nabla \psi|}, $$

it is thus important to define the first order derivative at the grid points. They are computed thanks to the standard centered scheme:

$$ \partial_x^h u_{i,j} = \frac{u_{i+1,j} - u_{i-1,j}}{2h}, \quad \partial_y^h u_{i,j} = \frac{u_{i,j+1} - u_{i,j-1}}{2h}, \quad \nabla^h = \begin{pmatrix} \partial_x^h \\ \partial_y^h \end{pmatrix}. $$

Since the level-set $\psi$ is defined in the whole domain $\Omega^h$, the above formula makes it possible to compute the normal vector to the isoline at the grid point $x_{i,j}$:

$$ \mathbf{n}_{i,j} = \frac{\nabla^h \psi_{i,j}}{|\nabla^h \psi_{i,j}|}. \qquad (10) $$

*Normal at the interface points.* The normal vector at the interface point is of great interest, especially to solve (1a)–(1b). Since $\psi$ is defined in the whole discrete domain, the normal vector at the interface is extrapolated from the neighboring values. For instance at the interface point $x_{i,j+\theta_y}^\Gamma$, the normal vector $\mathbf{n}_{i,j+\theta_y}$ is extrapolated<fnref n="4" /> by the values at the grid points $(x_{i,j+k})_{k=-1, \dots, 2}$:

$$ \mathbf{n}_{i,j+\theta_y} = C(\mathbf{n}_{i,j+2}, \mathbf{n}_{i,j+1}, \mathbf{n}_{i,j}, \mathbf{n}_{i,j-1}; [\theta_y - 2, \theta_y - 1, \theta_y, \theta_y + 1]), \quad (11) $$

where $\mathbf{n}_{i,j+k}$ are given by (10).

### 2.2.3. Direct and indirect first order derivatives at the interface

It is important to define appropriately the direct and indirect first order interface derivatives –denoted respectively by $\partial_x^\Gamma$ and $\partial_x^{\Gamma,I}$ – especially for Neumann problem. We introduce second order accurate interface derivatives, which are extensions of the first order definitions presented in [16], and which stabilize the stencils of [10, 20]. Note that this stable approach ensures the *stencil continuity*.

In the example of Figure 5(a), the point $x_{i,j+\theta_y}^\Gamma$ is on the y-axis and the indirect x-derivative of a function $u$ at this point computed with the stable ghost fluid method and cubic extrapolation of the ghost value:

$$ \partial_x^{\Gamma,I} u_{i,j+\theta_y} = \frac{u_{i+1,j+\theta_y} - u_{i-1,j+\theta_y}^G}{2h}, \quad (12) $$

where

$$ u_{i-1,j+\theta_y}^G = C((u_{i+k,j+\theta_y})_{k=0, \dots, 3}; [1, 2, 3, 4]). $$

<footnote n="4">⁴The above extrapolation has to be understood component by component, by abuse of notation.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5: Examples of the stencils for the interface normal derivative.</figcaption></figure>

Then the indirect interface derivative reads (see Figure 5):

$$ \partial_x^{\Gamma,I} u_{i,j+\theta_y}^{\Gamma} = \frac{-4 u_{i,j+\theta_y}^{\Gamma} + 7 u_{i+1,j+\theta_y} - 4 u_{i+2,j+\theta_y} + u_{i+3,j+\theta_y}}{2h}, \quad (13) $$

The above values $u_{i+k,j+\theta_y}$, $k = 0, \dots, 3$, have also to be extrapolated. Depending on the orientation of the normal vector at the point $x_{i,j+\theta_y}$, we choose $K \in \{0, 1, 2, 3, 4\}$ such that all the points $(x_{i+k,j+K-\ell})(k, \ell)=0, \dots, 3$ are grid points<fnref n="5" /> of $\omega^h$. Then the values $u_{i+k,j+\theta_y}$ for $k = 0, \dots, 3$ are approached in the y-direction as

$$ u_{i+k,j+\theta_y} = C((u_{i+k,j+K-\ell})_{\ell=0,\dots,3}; (\theta_y - K + \ell))_{\ell=0,\dots,3}). $$

The direct interface derivative is given by

$$ \partial_x^{\Gamma} u_{i-\theta_x,j} = \frac{u_{i+1-\theta_x,j} - u_{i-1-\theta_x,j}^G}{2h}, \quad (14) $$

where

$$ u_{i-1-\theta_x,j}^G = C((u_{i+k-\theta_x,j})_{k=0,\dots,3}; [1,2,3,4]). $$

Using the definition of the cubic extrapolation operator (4b), the direct derivative reads

$$ \partial_x^{\Gamma} u_{i-\theta_x,j}^{\Gamma} = \frac{-4 u_{i-\theta_x,j}^{\Gamma} + 7 u_{i+1-\theta_x,j} - 4 u_{i+2-\theta_x,j} + u_{i+3-\theta_x,j}}{2h}. \quad (15) $$

Here again the values $u_{i+k-\theta_x,j}$, $k = 1, \dots, 4$, are extrapolated from the neighboring grid point values, thanks to (9).

<footnote n="5">⁵In some very specific cases, it is impossible to have 4 aligned grid points of $\omega$: this means that the grid has to be refined and less accurate stencils have to be used. We do not consider such cases.</footnote>
**Remark 5** (Continuity of the stencils). The interesting novelty of stencils (13)-(15) lies in their continuity, in the sense that if the points $x_{i,j+\theta_y}^\Gamma$ and $x_{i-\theta_x,j}^\Gamma$ converge towards each other, the stencils of $\partial_x^{\Gamma,I} u_{i,j+\theta_y}^\Gamma$ and $\partial_x^{\Gamma} u_{i-\theta_x,j}^\Gamma$ are identical. Note that our method seems to have similarities with methods presented in previous works on wave equations [22] or on the evolution of surfactant concentration on a moving interface [21]. In particular, in these previous works the stencils satisfy such a continuity, but the authors do not mention any superconvergence property.

**Remark 6** (Treatment of the Neumann condition at the interface point). If Neumann condition is imposed on the interface the normal derivative is given thanks to the vector at the interface (11) and the above direct and indirect derivatives (13)-(15). For instance, at the interface point $x_{i,j+\theta_y}^\Gamma$, the normal derivative reads:

$$ \partial_{\mathbf{n}}^{\Gamma} u_{i,j+\theta_y} = \begin{pmatrix} \partial_x^{\Gamma,I} u_{i,j+\theta_y} \\ \partial_y^{\Gamma} u_{i,j+\theta_y} \end{pmatrix} \cdot \mathbf{n}_{i,j+\theta_y}. \qquad (16) $$

To simplify notation, let transiently denote by $n_x^\Gamma$ and $n_y^\Gamma$ the x and y components of $\mathbf{n}_{i,j+\theta_y}$. The keypoint to treat the Neumann boundary condition (1b) is that $\partial_{\mathbf{n}}^\Gamma u_{i,j+\theta_y}$ has the following form:

$$ \partial_{\mathbf{n}}^{\Gamma} u_{i,j+\theta_y} = \frac{1}{2h} \left[ \left( \alpha u_{i,j+\theta_y}^{\Gamma} + \sum_{i,j} a_{i,j} u_{i,j} \right) n_x^{\Gamma} + \left( \beta u_{i,j+\theta_y}^{\Gamma} + \sum_{i,j} b_{i,j} u_{i,j} \right) n_y^{\Gamma} \right] := -g_{i,j+\theta_y}^{\Gamma}, $$

where $\alpha, \beta$, $(a_{i,j})$ and $(b_{i,j})$ are the known extrapolation parameters, which depend on $\theta_x$ or $\theta_y$. Then, the interface value of the solution is given by

$$ u_{i,j+\theta_y}^{\Gamma} = - \frac{2h g_{i,j+\theta_y}^{\Gamma} + n_x^{\Gamma} \sum_{i,j} a_{i,j} u_{i,j} + n_y^{\Gamma} \sum_{i,j} b_{i,j} u_{i,j}}{\alpha n_x^{\Gamma} + \beta n_y^{\Gamma}}. \qquad (17) $$

Using the interface derivative discretizations as in [10] is sufficient to achieve the numerical gradient superconvergence. However, our approach stabilizes the method. Interestingly, it avoids the addition of interface unknowns in the discretization matrix as for Cisternino-Weynans' method.

Most of all, the continuous approach is absolutely required to achieve two levels of superconvergence at the order 2 for the solution $u$. Indeed, without stencil continuity, jumps of the values of the numerical second derivative appear, preventing the superconvergence on the second level in maximum norm.

Note that Robin condition can be handled by this method by adding a term to the denominator in relation (17).

In this subsection, we have presented all the operators needed to solve the static subproblems (1a)-(1b) and (1c)-(1d), provided that the level-set function $\psi$ is known at the grid points of the whole domain. In order to solve the coupled invadopodia model, it thus remains to define the extension and the transport operators.
## 2.3. Operators for the motion of the level set

### 2.3.1. Extended velocity, transport and volume conservation

In the Eulerian description, the transport of the level set is given by equation (3a). However, as the velocity is biologically given as the gradient of $\sigma$ in the inner area $O_t^i$, it does not exist in the outer area. In order to avoid the appearance of too high or discontinuous gradients near the interface, a continuous extension $v$ of this velocity is required from $\Gamma^h$ to its neighborhood. Similarly to in [16], the choice of the extension operator is given by (3b), meaning that the velocity is extended constantly along the normal to the isolines of $\psi$. Equation (3b) for velocity extension was originally introduced in [39, 1], and were solved thanks to the fast marching method.

Authors already noticed that the signed distance function was preserved if advected with this kind of extended velocity. This result follows from the derivation of Equation (3a):

$$ \partial_t \nabla \psi + \nabla (v \cdot \nabla \psi) = \partial_t \nabla \psi + (v \cdot \nabla) \nabla \psi + (\nabla \psi \cdot \nabla) v + v \times (\nabla \times \nabla \psi) + \nabla \psi \times (\nabla \times v) = 0 $$

In the right hand side, the third and fourth terms are zero by construction of the extended velocity and by curl operator property, respectively. The fifth term is colinear to the tangential vector. Multiplying by $\nabla \psi$, we therefore obtain

$$ \partial_t (|\nabla \psi|^2) + v \cdot \nabla (|\nabla \psi|^2) = 0 $$

Then, if the level set function is initialized as the signed distance function, we have $|\nabla \psi| = 1$ at the initial time, and the property is preserved at each time. Thus, the expensive procedure of redistanciation can be avoided, since the level set function is automatically reinitialized during advection.

The operators are discretized with a second order accurate truncature error, as described below. Assuming that the velocity is second order accurate, it directly follows that at each time

$$ |\nabla^h \psi| = 1 + \mathcal{O}(h^2). $$

Surprisingly, this conservation of the signed distance function with a second order accuracy seems to result in a superconvergent behavior of the level set function. More precisely, if advected with any second order extended velocity defined by (3b), the level set and its normal vector are both second order accurate. If in addition $\sigma$ has two levels of superconvergence, the extension preserves the superconvergence properties and the velocity therefore has one level of superconvergence. The level set function, its normal vector and its curvature are then second order accurate. These results are reported in Table 8 of Section 3.

Finally, it is worth noting that the extended velocity is not divergence-free. However, the volume is still well-conserved since only the interface velocity, which is extrapolated from the divergence-free velocity $\nabla \sigma$, is implied in volume conservation. A quick proof reads

$$ \partial_t \int_{O_t^i} d\mathbf{x} = \int_{\Gamma_t} \mathbf{v} \cdot \mathbf{n} \, ds = \int_{\Gamma_t} \partial_{\mathbf{n}} \sigma \, ds = \int_{O_t^i} \nabla \cdot \nabla \sigma \, d\mathbf{x} = 0. $$
From the numerical point of view, $\nabla \cdot \nabla \sigma$ is approximated with a second order accuracy in maximum norm, which results in a second order accurate volume conservation and explains the good result of Figure 3.

### 2.3.2. Discretization

The advection equation (3a) is approximated thanks to the standard second order Runge-Kutta method (RK2) for the time derivative and with the linear upwind differencing scheme (LUD), introduced by Warming and Beam [33], for the gradient of $\psi$. Since RK2 scheme in time does not involve any difficulty, we do not present it.

Even if the LUD scheme is well-known, details of a modified LUD scheme is presented below, for the discretization of velocity extension (3b). The modification is introduced for the points close to the interface, in order to take the interface value into account.

First, it is worth noting that the standard LUD discretization can be built from the usual centered scheme for the first derivatives and quadratic extrapolations of the downwind values<fnref n="6" />. In the above subsection 2.2 cubic extrapolations have been used. The choice of cubic extrapolations was justified by Remark 3 in order to have a chance to get two levels of superconvergence of the solutions to the static problems, as it will be investigated in the next section. However, it is known that for propagation equations, such cubic extrapolations of the downwind values generate numerical instabilities, and thus quadratic extrapolations have to be preferred. Let us mention that such quadratic extrapolations will not impact the superconvergence properties, since they are used for the velocity, which is given as the gradient of the solution to (1c)–(1d), and thus it is sufficient to focus on one level of superconvergence on the velocity.

Before giving the numerical schemes of upwind derivatives, let us first define the downwind value.

**Definition 7 (Downwind value).** Let $x_{i,j}$ and $x_{i \pm 1,j}$ be close grid points. The point $x_{i \pm j}$ is said to be downwind in the $x$-direction compared to $x_{i,j}$, if

$$ \pm n_{x_{i,j}} < 0. $$

Downwind points in the $y$-direction are similarly defined.

### 2.3.3. Upwind derivative based on LUD scheme

The upwind discrete derivative operators $\partial_x^{\text{UP}}$ and $\partial_y^{\text{UP}}$, applied to each of the velocity components, has to be upwind, with respect to each of the normal component (with the right sign). Here again we focus on the $x$-direction. At any point grid point $x_i$

$$ \partial_x^{\text{UP}} v_i = \frac{v_{i+1} - v_{i-1}}{2h}. $$

<footnote n="6">⁶Downwind values that appear in transport and velocity extension correspond to the ghost value of the static subproblems.</footnote>
Far from the interface (at the inner grid points), the downwind value is quadratically extrapolated with the 3 closest upwind value. For instance if $x_{i+1}$ is downwind (see Figure 6(a)), then

$$v_{i+1} = Q(v_i, v_{i-1}, v_{i-2}; [1, 2, 3]),$$

and thus $\partial_x^{\text{UP}} v_i$ reads

$$\partial_x^{\text{UP}} v_i = \frac{3 v_i - 4 v_{i-1} + v_{i-2}}{2h}. \qquad (18)$$

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Second order upwind schemes (based on LUD scheme) for velocity extension.</figcaption></figure>

At the point neighboring to the interface (see Figure 6(b)), it is necessary to include the interface value of the velocity, to account for the condition (3b). In order to keep a continuity of the stencil (see Remark 8), as in the previous subsection $v_{i+1}$ is extrapolated as

$$v_{i+1,j} = Q((v_{i+k-\theta_x})_{k=-1,0,1}; (\theta_x - k + 1)_{k=-1,0,1}),$$

and the derivative reads

$$\partial_x^{\text{UP}} v_i = \frac{[(1 + \theta_x)(2 + \theta_x) v_{i+\theta_x} - 2\theta_x(2 + \theta_x) v_{i-\theta_x} + \theta_x(1 + \theta_x) v_{i-1-\theta_x} - 2v_{i-1}]}{4h}. \qquad (19)$$

Then the intermediate points $v_{i\pm 1-\theta_x}$ are quadratically interpolated:

$$v_{i\pm 1-\theta_x} = Q((v_{i+k})_{k=-2,-1,0}; (\pm(1-k) - \theta_x)_{k=-2,-1,0}). \qquad (20)$$

**Remark 8 (Stencil continuity).** Note that this approach implies that in both cases $\theta_x \to 0$ and $\theta_x \to 1$, the scheme tends to be the same as (18). In particular, in the case $\theta_x \sim 1$, it gives the continuity between the stencil close to the interface ($\theta_x < 1$) and the stencil far from the interface ($\theta_x > 1$), which seems to be essential to ensure the preservation of the superconvergence property on the extended velocity. On the other side of the interface, the discretizations are similar, in the opposite direction since the normal component has the opposite sign. The results in subsection 3.3 confirm that the method of extension is second order accurate and preserves the velocity superconvergence.
### 2.4. Algorithm for the numerical computation

We are now ready to provide the algorithm to solve the coupled problem (1a)-(1b)-(1c)-(1d)-(3).

Note that if linear extrapolations of the ghost values give a symmetric matrix for the Dirichlet problem [18], it is not the case for the Neumann problem in general, or for the Dirichlet problem with our stencils. The non-symmetry of the matrix is somehow the cost to get superconvergence properties. In what follows, the non-symmetric matrices are inverted thanks to the BiCGStab method. For each simulation, the computational domain is $D = [-0.5, 0.5] \times [-0.5, 0.5]$.

The computation holds as follows: starting from $\psi^{t_n}$, the solution to the static subproblems $c^{*t_n}$ and $\sigma^{t_n}$ are computed thanks to the operators defined in subsection (2.2). Then the velocity at the interface, derived from $\nabla^T \sigma^{t_n}|_{\Gamma^{h,t_n}}$ is extrapolated constantly along the normals and the membrane is transported by solving (3) on one time step.

The above algorithm makes it possible to solve accurately the dynamical coupled problem (See Figure 2). In particular, the numerical convergence is shown to be at the order 2 for $c^*$, $\sigma$, $\psi$ and even $n$, as provided by table 1, showing the two levels of superconvergence of the unknowns: $\sigma$, $c^*$, and even $\psi$. In the next section superconvergence properties are investigated in more details.

## 3. Numerical study of the superconvergence properties

The aim of this section is to verify that the above schemes provide effectively superconvergent solutions. For such a purpose, it is necessary to compare first and second order derivatives of $\sigma$ and $c^*$ as well as first order derivatives of $v$ in the appropriate domains. The first order derivatives of the velocity have been already defined in the whole domain in subsection 2.3.3, however the first order derivatives of $c^*$ or $\sigma$ have been defined only at the interface. We thus introduce the first order derivatives at the grid points in subsection 3.1. Then we study the superconvergence properties of the static problems in subsection 3.2. We first study the superconvergence of each problem, solved separately, and we then exhibit the preservation of the superconvergence for the static coupled problem (1a)-(1b)-(1c)-(1d). We conclude this subsection by proposing the minimal requirements on the approximation of the data to obtain superconvergence on two levels and transmission of this property. Subsection 3.3 is devoted to the preservation of the superconvergent properties by the dynamics (3).

### 3.1. Numerical derivatives at the grid points of the inner or outer domains

If $x_i$ is an inner points of $\omega^h$, the standard centered stencil is used for the first order derivative $\partial_x^h u_i$ of a function $u$.

If $x_i$ is a neighboring point, then for instance $x_{i-1} \notin \omega$ and thus $u_{i-1}$ is a ghost value. Depending whether $u$ is the solution –namely $c^*$, or $\sigma$– or a first order derivative of the solution, the ghost
value is approached by cubic or quadratic extrapolations respectively. We thus define for $\mathbb{L} = \mathbb{Q}, \mathbb{C}$ the ghost value $u_{i-1}^{\mathbb{L}}$ as<fnref n="7" />

$$u_{i-1}^{\mathbb{L}} = \mathbb{L}((u_{i+k-1-\theta_x})_{k=1,\dots,n_{\mathbb{L}}}; \theta_{\mathbb{L}}), \quad \text{with} \quad \theta_{\mathbb{L}} = (k-\theta_x)_{k=1,\dots,n_{\mathbb{L}}}. \qquad (21)$$

The values $u_{i+k-\theta_x}$ are not at the grid points, and thus they are extrapolated at the order 2 or 3 (as noticed in Remark 3) from the values at the grid points:

$$u_{i+k-\theta_x} = \mathbb{L}((u_{i+\ell})_{\ell=0,\dots,n_{\mathbb{L}}-1}; \theta_{\mathbb{L}}), \quad \text{with} \quad \theta_{\mathbb{L}} = (k-\ell-\theta_x)_{\ell=0,\dots,n_{\mathbb{L}}-1}, \quad \text{for } k=0,\dots,n_{\mathbb{L}}-1,$$

where $\mathbb{L}$ stands for $\mathbb{Q}$, or $\mathbb{C}$. We then define the derivative operator $\partial_x^{\mathbb{L}}$ at the neighboring point $x_i$ as

$$\partial_x^{\mathbb{L}} u_i = \frac{u_{i+1} - u_{i-1}^{\mathbb{L}}}{2h}, \quad \mathbb{L} = \mathbb{Q}, \text{ or } \mathbb{C}.$$

We abusively denote by $\nabla_x^{\mathbb{L}}$ the first order derivative operators, corresponding to $\partial_x^h$ on the inner grid points and $\partial_x^{\mathbb{L}}$ on the grid points neighboring the interface. In addition $\nabla^{\mathbb{L}}$ stands for

$$\nabla^{\mathbb{L}} = \begin{pmatrix} \partial_x^{\mathbb{L}} \\ \partial_y^{\mathbb{L}} \end{pmatrix}, \quad \mathbb{L} = \mathbb{Q}, \text{ or } \mathbb{C}.$$

### 3.2. Superconvergence properties for the static problems

In this subsection, we one of the static subproblems (1a)-(1b) or (1c)-(1d). The set $\omega$ generically denotes either the outer area $O_t^e$ or the inner area $O_t^i$ of problem (1) at a given time $t$. $\Gamma$ denotes the location of $\Gamma_t$ at the time $t$. We consider the generic problem

$$\Delta u = f, \quad \text{in } \omega, \quad u|_{\partial\Omega\cap\partial\omega} = 0, \qquad (22a)$$

$$u|_{\Gamma \cap \partial \omega} = g, \text{ (Dirichlet), or } \partial_n u|_{\Gamma \cap \partial \omega} = g, \text{ (Neumann)}. \qquad (22b)$$

#### 3.2.1. Superconvergence for the Dirichlet problem (numerical test 1)

Consider the Dirichlet problem (22) whose exact solution given by

$$u(x, y) = x^2 \cos(\pi(4x + y)), \quad \text{in } \omega,$$

where $\omega = \{(x, y) : \psi(x, y) < 0\}$, where $\psi$ is given in polar coordinates:

$$\psi(r \cos \theta, r \sin \theta) = r + 0.4 \cos^2(\theta + \frac{\pi}{4}) \sin^2(\theta + \frac{\pi}{4}) - 0.3, \quad \forall (r, \theta) \in \mathbb{R}^+ \times \mathbb{R}/2\pi\mathbb{Z}. \qquad (23)$$

The functions $g$ and $f$ are explicitely inferred from $u$. The tests of numerical accuracy are given in Table 2, in discrete $L^1$-norm and maximum norm.

Note that Table 3 shows that the second-order partial derivatives and the second-order mixed derivatives (in the sense of two successive derivations) are also superconvergent.

<footnote n="7">⁷Remind that $n_Q = 3$ and $n_C = 4$.</footnote>
<table><caption>**Table 2:** Superconvergence of the solution for the Dirichlet problem (numerical test 1).</caption><thead><tr><th rowspan="2">Number of points</th><th colspan="4">∇<sup>C</sup><sub>u</sub></th><th colspan="4">∇<sup>Q</sup> · ∇<sup>C</sup><sub>u</sub></th></tr><tr><th>L<sup>1</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>1</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td>96</td><td>1.980 × 10<sup>-4</sup></td><td>-</td><td>3.034 × 10<sup>-3</sup></td><td>-</td><td>9.181 × 10<sup>-3</sup></td><td>-</td><td>1.096 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>144</td><td>8.667 × 10<sup>-5</sup></td><td>2.04</td><td>1.349 × 10<sup>-3</sup></td><td>2.00</td><td>4.069 × 10<sup>-3</sup></td><td>2.01</td><td>5.618 × 10<sup>-2</sup></td><td>1.65</td></tr><tr><td>216</td><td>3.842 × 10<sup>-5</sup></td><td>2.02</td><td>5.987 × 10<sup>-4</sup></td><td>2.00</td><td>1.805 × 10<sup>-3</sup></td><td>2.01</td><td>2.437 × 10<sup>-2</sup></td><td>1.85</td></tr><tr><td>324</td><td>1.720 × 10<sup>-5</sup></td><td>2.01</td><td>2.670 × 10<sup>-4</sup></td><td>2.00</td><td>8.088 × 10<sup>-4</sup></td><td>2.00</td><td>1.085 × 10<sup>-2</sup></td><td>1.90</td></tr><tr><td>486</td><td>7.656 × 10<sup>-6</sup></td><td>2.01</td><td>1.190 × 10<sup>-4</sup></td><td>2.00</td><td>3.601 × 10<sup>-4</sup></td><td>2.00</td><td>4.957 × 10<sup>-3</sup></td><td>1.91</td></tr></tbody></table>

<table><caption>**Table 3:** Accuracy of the second order partial and mixed derivatives of the solution for the Dirichlet problem (numerical test 1).</caption><thead><tr><th rowspan="2">Number of points</th><th colspan="4">∂<sub>x</sub><sup>Q</sup>∂<sub>x</sub><sup>C</sup><sub>u</sub></th><th colspan="4">∂<sub>y</sub><sup>Q</sup>∂<sub>y</sub><sup>C</sup><sub>u</sub></th></tr><tr><th>L<sup>1</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>1</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td>96</td><td>1.005 × 10<sup>-2</sup></td><td>-</td><td>1.984 × 10<sup>-1</sup></td><td>-</td><td>1.054 × 10<sup>-3</sup></td><td>-</td><td>1.277 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>144</td><td>4.415 × 10<sup>-3</sup></td><td>2.03</td><td>9.466 × 10<sup>-2</sup></td><td>1.82</td><td>4.278 × 10<sup>-4</sup></td><td>2.22</td><td>6.670 × 10<sup>-2</sup></td><td>1.60</td></tr><tr><td>216</td><td>1.950 × 10<sup>-3</sup></td><td>2.02</td><td>3.647 × 10<sup>-2</sup></td><td>2.09</td><td>1.779 × 10<sup>-4</sup></td><td>2.19</td><td>2.725 × 10<sup>-2</sup></td><td>1.90</td></tr><tr><td>324</td><td>8.721 × 10<sup>-4</sup></td><td>2.01</td><td>1.843 × 10<sup>-2</sup></td><td>1.95</td><td>7.812 × 10<sup>-5</sup></td><td>2.14</td><td>1.317 × 10<sup>-2</sup></td><td>1.87</td></tr><tr><td>486</td><td>3.880 × 10<sup>-4</sup></td><td>2.01</td><td>8.546 × 10<sup>-3</sup></td><td>1.94</td><td>3.449 × 10<sup>-5</sup></td><td>2.11</td><td>6.207 × 10<sup>-3</sup></td><td>1.86</td></tr></tbody><thead><tr><th colspan="4">∂<sub>x</sub><sup>Q</sup>∂<sub>y</sub><sup>C</sup><sub>u</sub></th><th colspan="4">∂<sub>y</sub><sup>Q</sup>∂<sub>x</sub><sup>C</sup><sub>u</sub></th></tr><tr><th>L<sup>1</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>1</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td>96</td><td>2.700 × 10<sup>-3</sup></td><td>-</td><td>3.417 × 10<sup>-1</sup></td><td>-</td><td>1.874 × 10<sup>-3</sup></td><td>-</td><td>1.794 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>144</td><td>9.307 × 10<sup>-4</sup></td><td>2.63</td><td>1.655 × 10<sup>-1</sup></td><td>1.79</td><td>6.805 × 10<sup>-4</sup></td><td>2.50</td><td>8.647 × 10<sup>-2</sup></td><td>1.80</td></tr><tr><td>216</td><td>3.419 × 10<sup>-4</sup></td><td>2.55</td><td>7.624 × 10<sup>-2</sup></td><td>1.85</td><td>2.610 × 10<sup>-4</sup></td><td>2.43</td><td>3.601 × 10<sup>-2</sup></td><td>1.98</td></tr><tr><td>324</td><td>1.293 × 10<sup>-4</sup></td><td>2.50</td><td>3.299 × 10<sup>-2</sup></td><td>1.92</td><td>1.089 × 10<sup>-4</sup></td><td>2.34</td><td>1.787 × 10<sup>-2</sup></td><td>1.90</td></tr><tr><td>486</td><td>5.116 × 10<sup>-5</sup></td><td>2.45</td><td>1.447 × 10<sup>-2</sup></td><td>1.95</td><td>4.572 × 10<sup>-5</sup></td><td>2.29</td><td>7.916 × 10<sup>-3</sup></td><td>1.92</td></tr></tbody></table>

### 3.2.2. Superconvergence for the Neumann problem (numerical test 2)

Consider $\Gamma$ implicitly defined as the level 0 of $\psi$ in polar coordinates:

$$ \psi(r \cos \theta, r \sin \theta) = -r - 0.25 \cos^2(\theta + \frac{\pi}{4}) + 0.4, \quad \forall (r, \theta) \in \mathbb{R}^+ \times \mathbb{R}/2\pi\mathbb{Z}. $$

The exact solution of the Neumann problem (22) is given by

$$ u(x, y) = x^2 \sin (2 \pi y), $$

in the outer area $\omega$, where $\psi < 0$. Table 4 shows the expected superconvergent results of the numerical accuracy tests in maximum norm.

<table><caption>**Table 4:** Superconvergence results for the Neumann problem (numerical test 2).</caption><thead><tr><th rowspan="2">Number of points</th><th colspan="2">u</th><th colspan="2">∇<sup>h</sup><sub>u</sub></th><th colspan="2">∇ · ∇<sup>h</sup><sub>u</sub></th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td>96</td><td>4.488 × 10<sup>-5</sup></td><td>-</td><td>1.096 × 10<sup>-3</sup></td><td>-</td><td>1.035 × 10<sup>-2</sup></td><td>-</td></tr><tr><td>144</td><td>1.992 × 10<sup>-5</sup></td><td>2.00</td><td>4.918 × 10<sup>-4</sup></td><td>1.98</td><td>4.633 × 10<sup>-3</sup></td><td>1.98</td></tr><tr><td>216</td><td>8.944 × 10<sup>-6</sup></td><td>1.99</td><td>2.201 × 10<sup>-4</sup></td><td>1.98</td><td>2.069 × 10<sup>-3</sup></td><td>1.99</td></tr><tr><td>324</td><td>3.970 × 10<sup>-6</sup></td><td>1.99</td><td>9.825 × 10<sup>-5</sup></td><td>1.98</td><td>9.222 × 10<sup>-4</sup></td><td>1.99</td></tr><tr><td>486</td><td>1.766 × 10<sup>-6</sup></td><td>1.99</td><td>4.380 × 10<sup>-5</sup></td><td>1.99</td><td>4.107 × 10<sup>-4</sup></td><td>1.99</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Plot of the numerical solution of the Neumann problem (Test-case 2).</figcaption></figure>

### 3.2.3. Minimal requirements on data to ensure superconvergence on two levels and transmission of superconvergence properties

We have numerically highlighted superconvergence properties on two levels for the solution of the Poisson's equation with an exact boundary condition. As the distance to the interface is computed with (5), the interface location is actually known with a fourth order accuracy. This small perturbation does not change significantly the results compared to an exact interface location. However, stronger perturbations would have lead to the loss of superconvergence. For instance, in the case of the Dirichlet problem, any third order perturbation only gives one level of superconvergence, while any second order perturbation does not give any superconvergence property. It is consistent with the fact that, as noted in several studies, the cornerstone of superconvergence lies on the Laplace operator consistancy. Perturbing the data is equivalent to use a pertubed Laplace operator, with exact data. At first glance, a fourth order perturbation on the level set function or on the Dirichlet boundary data implies a perturbed Laplace operator which keeps a second order truncation error, leading to the preservation of the superconvergence properties. A third order perturbation gives a first order truncation error near the interface and logically leads to the same superconvergent behavior on one level, as for the Shortley-Weller discretization. Finally, a second order perturbation implies the inconsistency of the operator near the interface, leading to second order accurate solution without superconvergence, as Gibou et al. For the Neumann problem, the reasoning is similar with one less order of accuracy on the numerical boundary data $g$, as suggested by (17).

In the case of the coupled problem, we cannot expect better than a second order accuracy on the data, especially $g$ and $\psi$, since each subproblem is solved with a second order method. Indeed, the second order solution of one subproblem gives the boundary data of the second problem, that cannot be more accurate than order 2.

In this subsection we exhibit specific second order perturbation of data that lead to superconvergence on two levels. More precisely, we highlight that the data can be only second order accurate, provided that they have themselves some superconvergent properties: this is the transmission of superconvergence properties.
Perturbation of data. Let $\varphi$ be a function of class $C^2(\omega)$, and suppose that its direct approximation $\varphi_{i,j}$ on the Cartesian grid has the following form

$$ \varphi_{ij} = \varphi(\mathbf{x}_{ij}) + \pi(\mathbf{x}_{ij}) h^2 + r_{ij} h^{2+p}, \quad \forall \mathbf{x}_{ij} \in \omega^h, \qquad (24) $$

where $\pi$ is a function of class $C^\infty(\omega)$ and $(r_{ij})_{i,j}$ are random numbers that ensure that derivatives of order higher than $p$ are not second order accurate. Hence $\varphi$ is superconvergent (towards $\pi$) on $p$ levels, for given $p \ge 1$.

**Remark 9** (Superconvergence on $\Gamma^h$). For the seek of simplicity, and by abuse of language, a function $g$ is said superconvergent on $\Gamma^h$ if it is the restriction of a superconvergent scalar field $G$ on $\omega^h$.

In what follows, expression (24) is used to perturb data of the problem and show the transmission of their superconvergent properties to the solution.

Perturbed Poisson's problems (numerical tests 3 and 4). Consider the same problems as subsections 3.2.1–3.2.2. Both Dirichlet and Neumann problems are still solved separately and the data are perturbed by

$$ f_{i,j} = f(x_{i,j}) + \mathcal{R}_f h^2, \quad \text{on } \omega^h, \qquad (25) $$

$$ \psi_{i,j} = \psi(x_{i,j}) + \sin(\pi x y) h^2 + \mathcal{R}_\psi h^p, \quad \text{on } \mathcal{P}^h, \qquad (26) $$

$$ g^\Gamma = g(x^\Gamma) + \cos\left(\frac{\pi}{2} x^\Gamma\right) h^2 + \mathcal{R}_g h^q, \quad \text{on } \Gamma^h, \qquad (27) $$

where $\mathcal{R}_f$, $\mathcal{R}_\psi$ and $\mathcal{R}_g$ are fields of random numbers of $[-5, 5]$. For both problems, we set $p = 4$, which means that the interface location, its normal vector and the curvature are second order accurate. The numerical function $f$ is also second order accurate in order to preserve the second order truncation error of the Laplace operator. The Dirichlet problem is solved by setting $q = 4$ so that the boundary data $g$ is superconvergent on two levels. For the Neumann problem, we just consider one level of superconvergence for the boundary datum by setting $q = 3$. Table 5 shows the results of the numerical accuracy tests. Surprisingly, the solution is superconvergent on two levels for both problems.

The same tests are successively performed with $p = 3$ and then using $q - 1$ instead of $q$ for each problem (see results in Table 6). In each of these numerical tests, one level of superconvergence is lost, which confirms the minimal requirements to achieve two levels of superconvergence.

Effect of the outer-inner coupling (numerical test 5). The consequence of the previous results is the preservation of the superconvergence by coupled problems. Solve now the static coupled problem (1a)-(1b)-(1c)-(1d) with the following continuous data:

$$ \psi(r \cos \theta, r \sin \theta) = r + 0.5 \cos^2 \theta \sin^2 \theta - 0.4, \quad \text{in } \mathcal{P}^h, \qquad (28) $$

$$ g = \sin^2(\pi (x + y)), \quad \text{on } \Gamma^h, \qquad (29) $$

$$ f = 0, \quad \text{in } \omega^h. \qquad (30) $$
<table><caption>**Table 5:** Transmission of superconvergent properties for the Dirichlet problem (numerical test 3, $p = q = 4$) and the Neumann problem (numerical test 4, $p = 4$, $q = 3$).</caption><thead><tr><th></th><th>Number<br>of points</th><th colspan="2">u</th><th colspan="2">∇<sup>Q</sup> · ∇<sup>C</sup><sub>u</sub></th><th colspan="2">∇<sup>Q</sup> · ∇<sup>C</sup><sub>u</sub></th></tr><tr><th></th><th></th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td rowspan="4">Dirichlet<br>problem</td><td>96</td><td>1.621 × 10<sup>-4</sup></td><td>-</td><td>3.161 × 10<sup>-3</sup></td><td>-</td><td>1.174 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>144</td><td>7.287 × 10<sup>-5</sup></td><td>1.97</td><td>1.335 × 10<sup>-3</sup></td><td>2.13</td><td>6.367 × 10<sup>-2</sup></td><td>1.51</td></tr><tr><td>216</td><td>3.248 × 10<sup>-5</sup></td><td>1.98</td><td>5.773 × 10<sup>-4</sup></td><td>2.10</td><td>2.467 × 10<sup>-2</sup></td><td>1.92</td></tr><tr><td>324</td><td>1.444 × 10<sup>-5</sup></td><td>1.99</td><td>2.580 × 10<sup>-4</sup></td><td>2.06</td><td>1.091 × 10<sup>-2</sup></td><td>1.95</td></tr><tr><td></td><td>486</td><td>6.423 × 10<sup>-6</sup></td><td>1.99</td><td>1.151 × 10<sup>-4</sup></td><td>2.04</td><td>5.011 × 10<sup>-3</sup></td><td>1.94</td></tr><tr><td rowspan="4">Neumann<br>problem</td><td>96</td><td>5.734 × 10<sup>-5</sup></td><td>-</td><td>1.105 × 10<sup>-3</sup></td><td>-</td><td>1.092 × 10<sup>-2</sup></td><td>-</td></tr><tr><td>144</td><td>2.594 × 10<sup>-5</sup></td><td>1.96</td><td>4.963 × 10<sup>-4</sup></td><td>1.97</td><td>5.021 × 10<sup>-3</sup></td><td>1.92</td></tr><tr><td>216</td><td>1.165 × 10<sup>-5</sup></td><td>1.97</td><td>2.219 × 10<sup>-4</sup></td><td>1.98</td><td>2.249 × 10<sup>-3</sup></td><td>1.95</td></tr><tr><td>324</td><td>5.235 × 10<sup>-6</sup></td><td>1.97</td><td>9.916 × 10<sup>-5</sup></td><td>1.98</td><td>1.008 × 10<sup>-3</sup></td><td>1.96</td></tr><tr><td></td><td>486</td><td>2.324 × 10<sup>-6</sup></td><td>1.98</td><td>4.419 × 10<sup>-5</sup></td><td>1.98</td><td>4.478 × 10<sup>-4</sup></td><td>1.97</td></tr></tbody></table>

<table><caption>**Table 6:** Loss of one level of superconvergence if the minimal requirements on data are not satisfied.</caption><thead><tr><th rowspan="3">Number<br>of points</th><th colspan="4">∇<sup>Q</sup> · ∇<sup>C</sup><sub>u</sub> in Dirichlet problem</th><th colspan="4">∇<sup>Q</sup> · ∇<sup>C</sup><sub>u</sub> in Neumann problem</th></tr><tr><th colspan="2">p = 3, q = 4</th><th colspan="2">p = 4, q = 3</th><th colspan="2">p = 3, q = 3</th><th colspan="2">p = 4, q = 2</th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td>96</td><td>1.411 × 10<sup>-1</sup></td><td>-</td><td>5.760 × 10<sup>-1</sup></td><td>-</td><td>1.121 × 10<sup>-2</sup></td><td>-</td><td>1.457 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>144</td><td>5.980 × 10<sup>-2</sup></td><td>2.12</td><td>4.744 × 10<sup>-1</sup></td><td>0.48</td><td>5.080 × 10<sup>-3</sup></td><td>1.95</td><td>1.546 × 10<sup>-1</sup></td><td>-0.15</td></tr><tr><td>216</td><td>3.796 × 10<sup>-2</sup></td><td>1.62</td><td>2.873 × 10<sup>-1</sup></td><td>0.86</td><td>3.857 × 10<sup>-3</sup></td><td>1.32</td><td>6.511 × 10<sup>-2</sup></td><td>0.99</td></tr><tr><td>324</td><td>2.341 × 10<sup>-2</sup></td><td>1.48</td><td>2.326 × 10<sup>-1</sup></td><td>0.75</td><td>3.661 × 10<sup>-3</sup></td><td>0.92</td><td>5.471 × 10<sup>-2</sup></td><td>0.81</td></tr><tr><td>486</td><td>1.450 × 10<sup>-2</sup></td><td>1.40</td><td>1.849 × 10<sup>-1</sup></td><td>0.70</td><td>2.561 × 10<sup>-3</sup></td><td>0.91</td><td>4.192 × 10<sup>-2</sup></td><td>0.77</td></tr></tbody></table>

The functions $\psi$ and $g$ are perturbed on the grid with (26) and (27), by setting $p = 4$ and $q = 3$, which are the minimal requirements to get superconvergence on two levels for the solution to the Neumann problem.

The only specificity of the Neumann-Dirichlet coupling lies in the interface transmission between the outer solution $c^*$ and the inner solution $\sigma$, which is performed by reconstructing the interface value of $c^*$ thanks to the formula (17). Figure 8 shows a plot of the numerical solutions computed on a 1000 × 1000 grid and which is used as a reference solution, since we do not know the analytic expression of the exact solutions $c^*$ and $\sigma$. As shown in Table 7, two levels of superconvergence are obtained for both outer and inner solutions, exhibiting thus a preservation of the superconvergence properties by the coupling.

### 3.3. Superconvergence preservation by the dynamics (numerical test 6)

The interface $\Gamma$ at the initial time is the circle of center $(0, 0)$ and radius $0.25$. Let consider an exact velocity $v$, which explicitly depends on the time variable $t$, and whose components are given in the inner area $O_t^i$ by

$$v_x = 0.8 \cos(2\pi(y + 0.35)) \sin(\pi x) \sin(\pi(0.6 - 2t)), \quad (31)$$

$$v_y = -0.5 \cos(\pi x) \sin(\pi(y^2 - 0.35^2)) \sin(\pi(0.5 - 2t)). \quad (32)$$
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8:** Numerical reference solution of the coupled problem (numerical test 5). Left: solution $c^*$ of the Neumann problem. Middle: solution $\sigma$ of the Dirichlet problem. Right: merged solutions.</figcaption></figure>

<table><caption>**Table 7:** Accuracy of the numerical solutions $c^*$ and $\sigma$ and their derivatives for the static part of Problem 1 (numerical test 5).</caption><thead><tr><th rowspan="2">Number of points</th><th colspan="2">u</th><th colspan="2">∇<sup>C</sup><sub>u</sub></th><th colspan="2">∇<sup>Q</sup> · ∇<sup>C</sup><sub>u</sub></th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th></tr></thead><tbody><tr><td rowspan="5">c<sup>*</sup></td><td>96</td><td>7.899 × 10<sup>-5</sup></td><td>-</td><td>2.148 × 10<sup>-3</sup></td><td>-</td><td>1.949 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>144</td><td>3.277 × 10<sup>-5</sup></td><td>2.17</td><td>9.569 × 10<sup>-4</sup></td><td>1.99</td><td>8.776 × 10<sup>-2</sup></td><td>1.97</td></tr><tr><td>216</td><td>1.716 × 10<sup>-5</sup></td><td>1.60</td><td>4.296 × 10<sup>-4</sup></td><td>1.98</td><td>4.874 × 10<sup>-2</sup></td><td>1.45</td></tr><tr><td>324</td><td>7.347 × 10<sup>-6</sup></td><td>2.09</td><td>2.029 × 10<sup>-4</sup></td><td>1.85</td><td>2.046 × 10<sup>-2</sup></td><td>2.14</td></tr><tr><td>486</td><td>2.857 × 10<sup>-6</sup></td><td>2.33</td><td>8.033 × 10<sup>-5</sup></td><td>2.28</td><td>8.871 × 10<sup>-3</sup></td><td>2.06</td></tr><tr><td rowspan="5">σ</td><td>96</td><td>7.685 × 10<sup>-5</sup></td><td>-</td><td>6.497 × 10<sup>-3</sup></td><td>-</td><td>1.291 × 10<sup>0</sup></td><td>-</td></tr><tr><td>144</td><td>3.254 × 10<sup>-5</sup></td><td>2.12</td><td>3.429 × 10<sup>-3</sup></td><td>1.58</td><td>6.496 × 10<sup>-1</sup></td><td>1.69</td></tr><tr><td>216</td><td>1.731 × 10<sup>-5</sup></td><td>1.56</td><td>1.681 × 10<sup>-3</sup></td><td>1.76</td><td>3.311 × 10<sup>-1</sup></td><td>1.66</td></tr><tr><td>324</td><td>7.634 × 10<sup>-6</sup></td><td>2.02</td><td>7.277 × 10<sup>-4</sup></td><td>2.06</td><td>1.538 × 10<sup>-1</sup></td><td>1.89</td></tr><tr><td>486</td><td>2.825 × 10<sup>-6</sup></td><td>2.45</td><td>2.820 × 10<sup>-4</sup></td><td>2.34</td><td>6.656 × 10<sup>-2</sup></td><td>2.06</td></tr></tbody></table>

The initial level set function is perturbed with (26) and $p = 4$, while at each time step the velocity components in $O_t^{i,h}$ are similarly perturbed thanks to

$$v_x = v_x + (1 - x^2) h^2 + \mathcal{R}_{v_x} h^3, \qquad (33)$$

$$v_y = v_y + 2y \sin(\pi x y) h^2 + \mathcal{R}_{v_y} h^3. \qquad (34)$$

The velocity is therefore superconvergent on one level at each time step. The level set function is advected until the final time $T = 0.5$. The errors on $\psi$, $v$ and their derivatives are computed at the final time by comparison with a reference simulation on a 1000 × 1000 grid, and in a tubular area around the level 0 of the level set function so as to avoid the ridges, which are generated by the velocity extension.

Table 8 shows the superconvergent results of the numerical accuracy tests in maximum norm, for the level set function $\psi$, the extended velocity $v$ and their derivatives.

Figure 9 shows plots of the interface location and the $x$-component of the extended velocity at the initial and final times. The point (0,0) in Figure 9(a) and the ridges in Figure 9(b) are singularities, which can disturb the advection or the computation of errors if they are too close to the interface. Note that it does not happen in the simulations presented in this study. Nonetheless,
<table><caption>**Table 8:** Preservation of superconvergence after velocity extension and advetion (numerical test 6).</caption><thead><tr><th rowspan="2">Number of points</th><th colspan="2">L<sup>∞</sup>-error</th><th colspan="2">L<sup>∞</sup>-error</th><th colspan="2">L<sup>∞</sup>-error</th><th colspan="2">L<sup>∞</sup>-error</th></tr><tr><th>order</th><th></th><th>order</th><th></th><th>order</th><th></th><th>order</th><th></th></tr></thead><tbody><tr><td>96</td><td>2.430 × 10<sup>-4</sup></td><td>-</td><td>1.009 × 10<sup>-2</sup></td><td>-</td><td>3.560 × 10<sup>-1</sup></td><td>-</td><td>5.934 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>144</td><td>9.738 × 10<sup>-5</sup></td><td>2.26</td><td>4.272 × 10<sup>-3</sup></td><td>2.12</td><td>1.626 × 10<sup>-1</sup></td><td>1.93</td><td>2.810 × 10<sup>-4</sup></td><td>1.84</td></tr><tr><td>216</td><td>3.976 × 10<sup>-5</sup></td><td>2.23</td><td>1.841 × 10<sup>-3</sup></td><td>2.10</td><td>7.563 × 10<sup>-2</sup></td><td>1.91</td><td>1.240 × 10<sup>-4</sup></td><td>1.93</td></tr><tr><td>324</td><td>1.466 × 10<sup>-5</sup></td><td>2.31</td><td>7.027 × 10<sup>-4</sup></td><td>2.19</td><td>2.822 × 10<sup>-2</sup></td><td>2.08</td><td>4.612 × 10<sup>-5</sup></td><td>2.10</td></tr><tr><td>486</td><td>5.729 × 10<sup>-6</sup></td><td>2.31</td><td>2.693 × 10<sup>-4</sup></td><td>2.23</td><td>1.000 × 10<sup>-2</sup></td><td>2.20</td><td>1.762 × 10<sup>-5</sup></td><td>2.17</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 9:** Interface location (black line) and *x*-component of the extended velocity ($w_x$) at the initial and final times.</figcaption></figure>

this problem can be overcome by occasionally reinitializing the level set function to a distance function [30], especially when the interface moves closer to an area where a ridge is previously appeared. In this case, the *redistanciation* should be performed at the order 4 in order not to degrade the superconvergence properties.

## 4. Conclusion

We have proposed a second order Cartesian method to solve the moving boundary problem of invadopodia. Our method is based upon the understanding and the full use of superconvergent properties on two levels and their transmission. These properties result from a specific method to solve Poisson equation, built thanks to a stable ghost fluid method and cubic extrapolations of the ghost values. This method can handle both Dirichlet and Neumann conditions on the boundary of the domain. Throughout numerical tests on 2D-domains, we have highlighted and classified the superconvergence properties, with respect to the boundary condition. We have also identified the required criteria of accuracy on the boundary data and on the interface location to achieve the superconvergence: the data do not have to be high-order accurate but just only superconvergent themselves, which is particularly suitable for coupling. If data are second order accurate and superconvergent, and despite the nonconsistency of the equivalent perturbed Laplace operator, the
solution, its gradient and the divergence of its gradient (and all the second order derivatives) are
also second order accurate. Moreover, these properties are stable through dynamics thanks to an
appropriate velocity extension and a second order scheme for advection.

The transmission of superconvergence properties is one of the main results of the paper. While superconvergence has been studied and proved in other works for the case of exact data, the transmission of superconvergence properties from second order superconvergent data to the solution is rarely addressed and is still an open issue. We are confident that this study is a step towards the theoretical understanding of this transmission phenomenon.

We emphasize that the method may be useful to many coupled problems similar to invadopodia.
More precisely, this is well suited to solve free boundary problems of Hele-Shaw type, in which the
interface velocity is the gradient of an harmonic function.

Obviously, our method use wide stencils. The parallel implementation of the code –for a three-
dimensional use for instance– is complex to carry out, requiring four communication layers. This is
the cost to get superconvergence properties on two levels, which result in the second order accuracy
of the solutions, including the curvature. This point is of great interest for further modeling of cell
migration processes as protrusion and rear retraction, using interface regularization techniques.

However, if second order accuracy and/or two levels of superconvergence are not required,
our study can be generalized to similar methods, with less accurate extrapolations. Based on a
previous study in [15], results of superconvergence and minimal requirements on data for each case
(linear, quadratic and cubic extrapolations) are summarized, without details, in Table 9. Thus,
each problem that requires superconvergence can be solved with the most suitable method (and the
least computationally expensive as possible), with respect to the desired accuracy and the coupling
degree. In particular, the first order accuracy observed in [16] can be explained by those results.

Finally, we must mention many remaining issues, which are a crucial point for concrete appli-
cations. First, the behavior of the time scheme is not yet well understood, possibly because the
stability condition depends on the spatial discretization. In the cubic case, the discretization is not
only second order accurate but also superconvergent on two levels. This could require a nonlinear
CFL, in the same spirit as in [11]. Using this kind of condition, some instabilities encountered on
very refined meshes disappear. The CFL analysis is therefore essential.

Another issue is related to the appearance of many pathological configurations, with too few
points, which may lead to a local loss of superconvergence. Similarly, the emergence of strong
curvatures and shocks near the interface during advection of the level set can also impact the
results of convergence and the accuracy of solutions. These situations generally arise because of a
too coarse grid, that is not always possible to refine, and because of the fixed nature of Cartesian
grid. For all these sensitive cases, patches are possible but often cumbersome to implement. For
this reason, a real challenge would be to extend this study to adaptive meshes, such as quadtree
and octree meshes.
<table><caption>**Table 9:** Zoology of the Cartesian numerical methods to solve the Poisson equation, and their superconvergence properties, with minimal requirements on data accuracy, based on [15]. L-scheme, Q-scheme and C-scheme are built with linear, quadratic and cubic extrapolations, respectively. The interface is implicitly defined by the level set function $\psi$ and $G$ denotes any lifting of the boundary datum $g$. $L$ stands for $L$ in the case of the L-scheme and $(\tilde{L}, L)$ stands for either $(L, Q)$ or $(Q, C)$ respectively for either the Q-scheme or the C-scheme.</caption><thead><tr><th></th><th colspan="3">Dirichlet</th><th colspan="3">Neumann</th></tr><tr><th></th><th>L-scheme</th><th>Q-scheme</th><th>C-scheme</th><th>L-scheme</th><th>Q-scheme</th><th>C-scheme</th></tr></thead><tbody><tr><th colspan="7">Order of accuracy required on data</th></tr><tr><th>G</th><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><th>∇<sup>L</sup>G</th><td></td><td>2</td><td>2</td><td></td><td></td><td>2</td></tr><tr><th>∇<sup>˜L</sup> · ∇<sup>L</sup>G</th><td></td><td></td><td>2</td><td></td><td></td><td></td></tr><tr><th>ψ</th><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><th>∇<sup>L</sup>ψ</th><td></td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><th>∇<sup>˜L</sup> · ∇<sup>L</sup>ψ</th><td></td><td></td><td>2</td><td></td><td></td><td>2</td></tr><tr><th colspan="7">Order of accuracy of the solution</th></tr><tr><th>u</th><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><th>∇<sup>L</sup>u</th><td>1</td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><th>∇<sup>˜L</sup> · ∇<sup>L</sup>u</th><td></td><td>1</td><td>2</td><td></td><td>1</td><td>2</td></tr></tbody></table>

## Acknowledgements

This study has been supported by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux - CPU (ANR-10-IDEX-03-02) and by the Core program Advanced Research Networks, granted by the JSPS Core.

O.G. and C.P. are partly granted by the Plan Cancer DYNAMO (Inserm 9749) and Plan Cancer NUMEP (Inserm 11099). This study has been carried out within the scope of the Inria Associate Team Num4SEP.
References

[1] D Adalsteinsson and J.A Sethian. The fast construction of extension velocities in level set methods. *J. Comput. Phys.*, 148(1):2–22, January 1999.

[2] Ivo Babuška, Uday Banerjee, and John E. Osborn. Superconvergence in the generalized finite element method. *Numerische Mathematik*, 107(3):353–395, 2007.

[3] Chunjia Bi. Superconvergence of finite volume element method for a nonlinear elliptic problem. *Numerical Methods for Partial Differential Equations*, 23(1):220–233, 2007.

[4] J. H. Bramble and B. E. Hubbard. On the formulation of finite difference analogues of the Dirichlet problem for Poisson’s equation. *Numerische Mathematik*, 4(1):313–327, 1962.

[5] O. Bruno and D. Hoch. Numerical differentiation of approximated functions with limited order-of-accuracy deterioration. *SIAM Journal on Numerical Analysis*, 50(3):1581–1603, 2012.

[6] Han Chen, Chohong Min, and Frédéric Gibou. A Supra-Convergent Finite Difference Scheme for the Poisson and Heat Equations on Irregular Domains and Non-Graded Adaptive Cartesian Grids. *Journal of Scientific Computing*, 31(1–2):19–60, 2007.

[7] Han Chen, Chohong Min, and Frédéric Gibou. A numerical scheme for the Stefan problem on adaptive Cartesian grids with supralinear convergence rate. *Journal of Computational Physics*, 228(16):5803–5818, 2009.

[8] So-Hsiang Chou and Xiu Ye. Superconvergence of finite volume methods for the second order elliptic problem. *Computer Methods in Applied Mechanics and Engineering*, 196(37–40):3706–3712, 2007. Special Issue Honoring the 80th Birthday of Professor Ivo Babuška.

[9] Philippe G. Ciarlet. Discrete maximum principle for finite-difference operators. *aequationes mathematicae*, 4(3):338–352, 1970.

[10] Marco Cisternino and Lisl Weynans. A parallel second order Cartesian method for elliptic interface problems. *Communications in Computational Physics*, 12:1562–1587, June 2012.

[11] Erwan Deriaz. Stability conditions for the numerical solution of convection-dominated problems with skew-symmetric discretizations. *SIAM J. Numer. Anal.*, 50(3):1058–1085, 2012.

[12] Ronald P. Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method). *Journal of Computational Physics*, 152(2):457–492, 1999.

[13] José A. Ferreira and Rolf D. Grigorieff. On the supraconvergence of elliptic finite difference schemes. *Applied Numerical Mathematics*, 28(24):275–292, 1998.
[14] A. Friedman. *Variational principles and free-boundary problems*. Pure and applied mathematics. Wiley, 1982.

[15] Olivier Gallinato. *Modeling of cancer phenomena and superconvergent methods for the resolution of interface problems on Cartesian grid*. Theses, Université de Bordeaux, November 2016.

[16] Olivier Gallinato, Masahito Ohta, Clair Poignard, and Takashi Suzuki. Free boundary problem for cell protrusion formations: theoretical and numerical aspects. *Journal of Mathematical Biology*, pages 1–45, 2016.

[17] Frédéric Gibou and Ronald Fedkiw. A Fourth Order Accurate Discretization for the Laplace and Heat Equations on Arbitrary Domains, with Applications to the Stefan Problem. *J. Comput. Phys.*, 202(2):577–601, January 2005.

[18] Frédéric Gibou, Ronald P. Fedkiw, Li-Tien Cheng, and Myungjoo Kang. A Second-Order-Accurate Symmetric Discretization of the Poisson Equation on Irregular Domains. *Journal of Computational Physics*, 176:205–227, 2002.

[19] Thomas Y. Hou, Zhilin Li, Stanley Osher, and Hongkai Zhao. A hybrid method for moving interface problems with application to the heleshaw flow. *Journal of Computational Physics*, 134(2):236 – 252, 1997.

[20] Z. Jomaa and C. Macaskill. The ShortleyWeller embedded finite-difference method for the 3D Poisson equation with mixed boundary conditions. *Journal of Computational Physics*, 229(10):3675–3690, 2010.

[21] Shilpa Khatri and Anna-Karin Tornberg. An embedded boundary method for soluble surfactants with interface tracking for two-phase flows. *Journal of Computational Physics*, 256:768 – 790, 2014.

[22] Heinz-Otto Kreiss, N. Anders Petersson, and Jacob Yström. Difference approximations of the neumann problem for the second order wave equation. *SIAM Journal on Numerical Analysis*, 42(3):1292–1323, 2004.

[23] Nick Levine. Superconvergent Recovery of the Gradient from Piecewise Linear Finite-element Approximations. *IMA Journal of Numerical Analysis*, 5(4):407–427, 1985.

[24] Zi-Cai Li, Hsin-Yun Hu, Qing Fang, and Tetsuro Yamamoto. Superconvergence of Solution Derivatives for the ShortleyWeller Difference Approximation of Poisson’s Equation. II. Singularity Problems. *Numerical Functional Analysis and Optimization*, 24(3-4):195–221, 2003.
[25] Zi-Cai Li, Hsin-Yun Hu, Song Wang, and Qing Fang. Superconvergence of solution derivatives of the Shortley-Weller difference approximation to Poisson's equation with singularities on polygonal domains. *Applied Numerical Mathematics*, 58(5):689-704, 2008.

[26] Zi-Cai Li, Tetsuro Yamamoto, and Qing Fang. Superconvergence of solution derivatives for the Shortley-Weller difference approximation of Poisson's equation. part I: smoothness problems. *Journal of Computational and Applied Mathematics*, 151(2):307-333, 2003.

[27] Nami Matsunaga and Tetsuro Yamamoto. Superconvergence of the Shortley-Weller approximation for Dirichlet problems. *Journal of Computational and Applied Mathematics*, 116(2):263-273, 2000.

[28] Chohong Min, Frédéric Gibou, and Hector D. Ceniceros. A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids. *Journal of Computational Physics*, 218(1):123-140, 2006.

[29] S. Osher and J.A. Sethian. Fronts Propagating with Curvature Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations. *Journal of Computational Physics*, 79(1):12-49, 1988.

[30] Stanley Osher and Ronald Fedkiw. *Level set methods and dynamic implicit surfaces*, volume 153 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 2003.

[31] Joseph Papac, Frédéric Gibou, and Christian Ratsch. Efficient symmetric discretization for the Poisson, heat and Stefan-type problems with Robin boundary conditions. *Journal of Computational Physics*, 229(3):875-889, 2010.

[32] George H. Shortley and Royal Weller. The Numerical Solution of Laplace's Equation. *Journal of Applied Physics*, 9(5):334, 1938.

[33] R. F. Warming and Richard M. Beam. Upwind second-order difference schemes and applications in aerodynamic flows. *AIAA Journal*, 14(9):1241-1249, 2016/07/06 1976.

[34] Lisl Weynans. A proof in the finite-difference spirit of the superconvergence of the gradient for the Shortley-Weller method. Inria Research Report RR-8757. https://hal.inria.fr/hal-01176994.

[35] Tetsuro Yamamoto. *Symbolic Algebraic Methods and Verification Methods*, chapter A New Insight of the Shortley-Weller Approximation for Dirichlet Problems, pages 245-253. Springer Vienna, Vienna, 2001.

[36] Gangjoon Yoon and Chohong Min. A review of the supra-convergences of Shortley-Weller method for Poisson equation. *J. KSIAM*, 18(1):51-60, 2014.
[37] Gangjoon Yoon and Chohong Min. Convergence Analysis of the Standard Central Finite Difference Method for Poisson Equation. *Journal of Scientific Computing*, pages 1–16, 2015.

[38] Tie Zhang and Shuhua Zhang. The gradient superconvergence of the finite volume method for a nonlinear elliptic problem of nonmonotone type. *Applications of Mathematics*, 60(5):573–596, 2015.

[39] Hong-Kai Zhao, T. Chan, B. Merriman, and S. Osher. A variational level set approach to multiphase motion. *Journal of Computational Physics*, 127(1):179 – 195, 1996.

[40] Milős Zlámal. Superconvergence and reduced integration in the finite element method. *Mathematics of Computation*, 32(143):663–685, 1978.