Stability and error analysis for a diffuse interface approach to an advection–diffusion equation on a moving surface

In this paper we analyze a fully discrete numerical scheme for solving a parabolic PDE on a moving surface. The method is based on a diffuse interface approach that involves a level set description of the moving surface. Under suitable conditions on the spatial grid size, the time step and the interface width we obtain stability and error bounds with respect to natural norms. Furthermore, we present test calculations that confirm our analysis.


Introduction
Let {Γ (t)} t∈[0,T ] be a family of closed hypersurfaces in R n+1 (n = 1, 2) evolving in time. In this paper we consider a finite element approach for solving the parabolic surface PDE equation which models advection and diffusion of a surface quantity u with u(·, t) : Γ (t) → R.
Here, S T = t∈(0,T ) Γ (t) × {t} and v : S T → R n+1 denotes a given velocity field. Furthermore, ∇ Γ is the tangential gradient, Δ Γ = ∇ Γ · ∇ Γ the Laplace Beltrami operator and ∂ • t = ∂ t + v · ∇ denotes the material derivative. Parabolic surface PDEs of the form (1) have applications in fluid dynamics and materials science, such as the transport and diffusion of surfactants on a fluid/fluid interface, [25] or diffusion-induced grain boundary motion, [5]. In these as in several other applications the velocity v is not given but determined through an additional equation so that (1) becomes a subproblem of a more complicated system in which the variable u is coupled to other variables. The analysis and the numerical solution of such systems then naturally requires the development of corresponding methods for (1). We refer to [13] for a comprehensive overview of finite element methods for solving PDEs on stationary and evolving surfaces.
Concerning the numerical methods that have been proposed for (1) one may distinguish between Lagrangian and Eulerian type schemes. The first approach has been pursued by Dziuk and Elliott within their evolving surface finite element method, [8], which uses polyhedral approximations of the evolving hypersurfaces Γ (t). While [8] contains an error analysis in the spatially discrete case, the fully discrete case is investigated in [11,14] and [19]. Optimal L 2 -error bounds are obtained in [12] and a corresponding finite volume approach is proposed and analyzed in [18]. Since the mesh for the discretization of (1) is fitted to the hypersurface Γ (t), a coupling to a bulk equation is not straightforward. This difficulty is not present in Eulerian type schemes, in which Γ (t) is typically described via a level set function defined in an open neighbourhood of Γ (t). In order to discretize the surface PDE in this setting it has been proposed in [1,3] and [27] to extend the surface quantity u to a band around Γ (t) and to solve a suitable (weakly) parabolic PDE in that bulk region using a finite difference method. In [9] and [10], the same idea is used in a finite element context for which the underlying variational formulation is derived with the help of a transport identity. An Eulerian finite element approach that doesn't use an extended PDE is proposed and analyzed in [20] and [21]. The method is based on a weak formulation on the space-time manifold and the finite element space is obtained by taking traces of the corresponding bulk finite elements. The approximation of Γ (t) on which these spaces are defined usually arises from a suitable interpolation of the given level set function describing Γ (t). The resulting discrete hypersurface will in general cut arbitrarily through the background mesh and its location forms one of the main difficulties in implementing the scheme. A different approach of generating the discrete hypersurfaces is pursued in [17], where a discretization of (4) below is combined with the cut finite element technique. Finally, Section 5 in [7] proposes a hybrid method that employs the above-mentioned idea of trace finite elements together with a narrow band technique for the elliptic part of the PDE.
In this paper we are concerned with the diffuse interface approach for solving (1), which was introduced in [22] for a stationary surface and in [16,23] and [26] for evolving surfaces. As in some of the methods described above, the surface quantity u is extended to a bulk quantity satisfying a suitable parabolic PDE in a neighbourhood of Γ (t) and the bulk equation is then localized to a thin layer of thickness with the help of a phase field function (see [15] for a corresponding convergence analysis). Since we are interested in using finite elements, the localized PDE needs to be written in a suitable variational form. Following [16] this is achieved with the help of a transport identity and results in a discretization by linear finite elements in space and a backward Euler scheme in time. The detailed derivation along with an existence result for the discrete solution will be given in Sect. 3. As the main new contribution of our paper we shall derive conditions relating the interface width , the spatial grid size h and the time step τ which allow for a rigorous stability and error analysis. More precisely, we shall prove that the numerical solution is bounded uniformly in L ∞ (L 2 ) and L 2 (H 1 ) over the diffuse interface (see Theorem 1 in Sect. 4) and that it converges with respect to these norms both over the diffuse interface and on the sharp interface with an order see Theorem 2 and Corollary 1 in Sect. 5 respectively. In Sect. 6 we report on results of numerical tests both for n = 1 and n = 2. An advantage of our approach is that in the implementation the evolution of the hypersurfaces is easily incorporated by evaluating the phase field function. We shall employ a function with compact support, namely ρ(x, t) := g( φ(x,t) ), where Γ (t) is the zero level set of φ(·, t) and In view of the evolution of the hypersurfaces the numerical scheme then naturally contains terms in which ρ is evaluated at different times. One of the main challenges in the analysis is to handle the corresponding differences, for which one has to bound integrals that are multiplied by a negative power of (arising from derivatives of ρ) as well as integrals that are not weighted with ρ. We shall deal with these difficulties by introducing an additional stabilization term with extended support that is also used for proving the well-posedness of the scheme. Let us finally remark that a phase field approach involving a phase field function with noncompact support and finite elements has been proposed in [4] for an elliptic surface PDE. Theorem 7 in [4] provides an error estimate in terms of an approximation error and an error due to the phase field representation. The latter decays at a rate O( p ) for some p < 1, while a coupling between and the grid size h is not discussed.

Surface representation and surface derivatives
For each t ∈ [0, T ] let Γ (t) ⊂ R n+1 (n = 1, 2) be a connected, compact and orientable hypersurface without boundary. We suppose that v : S T → R n+1 is a prescribed velocity field of the form Here, ν is a unit normal and V the corresponding normal velocity of Γ (t) and (·, ·) denotes the Euclidian scalar product in R n+1 . Note that the normal part V ν is responsible for the geometric motion of Γ (t), while the tangential part v τ is associated with the transport of material along the surface. We assume that there exists a smooth map Ψ : Let us next introduce the differential operators which are required to formulate our PDE. To begin, for fixed t and a function η : Γ (t) → R we denote by ∇ Γ η = (D 1 η, . . . , D n+1 η) its tangential gradient. Ifη is an extension of η to an open neighbourhood of Γ (t) then Furthermore, Next, for a smooth function η on S T we define the material derivative of η at Ifη is an extension of η to an open space-time neighbourhood, then Our numerical approach will be based on an implicit representation of Γ (t), so that we suppose in what follows that there exists a smooth function φ : Here, Ω ⊂ R n+1 is a bounded domain with Γ (t) ⊂ Ω for all t ∈ [0, T ]. For later use we introduce for t ∈ [0, T ], r > 0 the sets In view of (7) there exist

Extension
Our next aim is to extend functions defined on S T to a space-time neighbourhood. A common approach which is well suited to a description of Γ (t) via the signed distance function consists in extending constantly in the normal direction. In what follows we shall introduce a suitable generalization to the case (7). Consider for P ∈ Γ (0) and t ∈ [0, T ] the parameter-dependent system of ODEs Using a compactness argument it can be shown that there exists 0 < δ < δ 0 so that the solution γ P,t of (9) exists uniquely on (−δ, δ) uniformly in P ∈ Γ (0), t ∈ [0, T ].
Let us next extend the surface differential operators ∇ Γ and ∂ • t . By reversing the orientation of Γ (t) if necessary we may assume that the functions ν : are extensions of the unit normal and the normal velocity respectively. In particular, we define for a function η ∈ C 1 (U δ (t)) its Eulerian tangential gradient by and remark that ( Note that H |Γ (t) is the mean curvature of Γ (t).
Let us also extend the velocity field v to U δ,T . We first extend its tangential part by settingṽ In view of (3) the function v(x, t) := V (x, t)ν(x, t) +ṽ τ (x, t) extends the given velocity field from S T to U δ,T and satisfies In particular, we can use the extended velocity v to define the material derivative for a function η on U δ,T by setting

Phase field approach
Consider for 0 < < 2δ π the function where g ∈ C 1,1 (R) is given by . Furthermore, we obtain from the definition of ∇ φ and (24) The phase field function ρ allows us to approximate the integration over a surface Γ (t) in terms of a volume integral over the diffuse interface. More precisely, for fixed for small > 0, so that we can view 2 π Ω η ρ(·, t) |∇φ(·, t)| dx as an approximation of Γ (t) η dH n . This formula explains the appearance of the weight ρ(·, t) |∇φ(·, t)| in subsequent volume integrals.
In what follows we shall make use of the following continuity properties of ρ.
Using the mean value theorem and (8) we then have i.e. x ∈ U 3 π 4 (t). In order to prove (27) and (28) we first observe that it is enough to verify the estimates for x ∈ U 3 π 4 (t) in view of what we have just shown. There exists ξ between s and t such that by (8). Furthermore, since we see immediately that As a result, Inserting this bound into (29) yields (27). Finally, using again (30) and (8) we obtain Remark 1 Our analysis will work for other profile functions g than the one chosen above as long as they satisfy g ∈ C 1,1 (R) and g(r Profile functions with noncompact support have been used in [4,22] and [26]. However it is not obvious how to extend the analysis presented below to that setting.

Discretization
Suppose that u is a smooth solution of (1). It is shown in Lemma 8 of the "Appendix" that its extension u e satisfies the strictly parabolic PDE where and b lk , c k , d are smooth functions depending on φ and v.
In order to associate with (31) a suitable variational formulation we adapt an idea from [16], which uses an Eulerian transport identity. More precisely, we infer with the help of Lemma 3 in [10], (26) and (31) Here, the last equality follows from integration by parts together with the fact that (∇u e , ∇ρ) = 1 g φ (∇u e , ∇φ) = 0 in view of (15).
Let us first discretize with respect to time and denote by 0 Under a suitable regularity assumption on u we have that |φ R| ≤ C on suppρ so that we neglect the corresponding term when now deriving the spatial discretization from (34).
In what follows we assume that Ω is polyhedral and consider a family (T h ) 0<h≤h 0 of triangulations of Ω with mesh size h = max T ∈T h h T , h T = diam(T ). We assume that the family is regular in the sense that there exists σ > 0 with where r T is the radius of the largest ball contained in T . Let us denote by N h the set of vertices of the triangulation T h . In order to formulate our scheme we require a second phase field function with a slightly larger support, namelỹ For 0 ≤ m ≤ M we then define as well as the finite element space We denote by I m h : Then a) U 3 π Proof a) Let x ∈ D m h , so that there exists y ∈ N h such that |y − x| ≤ h and ρ m (y) > 0. Hence |φ m (y)| < π and the mean value theorem together with (8) yields for s ∈ [max(t m−1 , 0), min(t m+1 , T )] Thenρ(x, t) ≥ cos 2 ( 3π 8 ) and we obtain similarly as above In particular, [I m hρ m ](x) > 0, so that x ∈ D m h . Using the above inequality for t = t m implies b).
Our finite element approximation of (1), (2) now reads : For m = 1, 2 The parameter γ > 0 will be chosen in such a way as to ensure existence and stability for the scheme, see Lemma 5 and Theorem 1 below. b) The term γ τ 2 m Ω I m hρ m (∇u m h , ∇v h ) introduces artificial diffusion into the scheme and will play a crucial role in our analyis. A different form of stabilization is used in [16], Section 2.5. c) Unlike the schemes introduced in [16] our method is not fully practical because we assume that the integrals are evaluated exactly. In Sect. 6 we shall follow [16] in using numerical integration to obtain a fully practical scheme. A nice feature of the resulting method is that the evolution of the hypersurfaces is tracked in a simple way via the evaluation of ρ.
In what follows we shall be concerned with the existence, stability and error bounds for (37). The extension of our analysis to the fully practical method mentioned above is currently out of reach and left for future research. However, the test calculations in Sect. 6 show that the parameter choices suggested by the analysis work well also for the fully practical scheme.

Lemma 4 There exists
Proof To begin, we remark that there exists 0 < h 1 ≤ h 0 and μ > 0 only depending on σ, c 0 , c 1 , c 2 such that for every a ∈ N h ∩ U δ (t) there exists a neighbour is connected it is sufficient to show that for every y ∈ D m h there exists z ∈ Γ (t m ) and a path in D m h connecting y to z. Let us fix y ∈ D m h , say y ∈ T , whereρ m (x) > 0 for some x ∈ T ∩ N h . We assume w.l.o.g. that 0 < φ m (x) < π. In view of (39) there exists a neighbour

Lemma 5 (Existence) Let
There exists τ 0 > 0 such that the scheme (37) has a unique solution u m h ∈ V m h provided that 0 < τ ≤ τ 0 . Proof Since (37) is equivalent to solving a linear system with a quadratic coefficient matrix, it is sufficient to prove that the following problem only has the trivial solution: If we choose τ 0 > 0 so small that 1 2 which implies that u h ≡ 0 on Γ (t m ) and ∇u h ≡ 0 in D m h . According to Lemma 4,D m h is connected, so that we conclude that u h ≡ 0.

Stability bound
The following lemma will be useful in estimating L 2 -integrals that are not weighted by ρ.
It follows from Theorem 4.4 in [8] (extended in a straightforward way to the case of a nontrivial f ) that (1), (2) has a unique solution u which satisfies The following theorem gives a discrete version of this estimate in the phase field setting.
Proof Setting v h = u m h in (37) we find after a straighforward calculation Clearly, while Integrating by parts and abbreviating H m = −∇ · ν m we obtain In order to rewrite I I I we first observe that in view of (22) and (24) ( so that (23), (25) and again (24) imply Inserting (44)-(46) into (43) we infer that We deduce from (28), Lemma 6, (41) and the assumption τ ≤ 2 that if we choose γ ≥ γ 1 := C + 1. Finally, using Taylor expansion and (8) we infer that Inserting the above estimates into (47) we find If τ 1 ≤ τ 0 is sufficiently small we therefore deduce for τ ≤ τ 1 from which we obtain after summation from m = 1, . . . , l and division by that Using Lemma 3 a), (38) and (16) we may estimate Arguing in a similar way for the term involving f e,m we derive The discrete Gronwall inequality yields the bound on max m=1,...,M 1 Ω (u m h ) 2 ρ m |∇φ m |, which combined with (49) implies the second inequality. (14). Then we have for m = 1, . . . , M and t ∈ [t m−1 , t m ]:

Lemma 7 Suppose that (36) holds and let z e be defined by
.
. Standard interpolation theory together with Lemma 3 a) and (16) implies that The second bound follows in the same way using (17).

Theorem 2
Suppose that the solution of (1), (2) satisfies Inserting v h = e m h and following the argument in the proof of Theorem 1 leading to (48) we obtain We now deal individually with the terms S m i , e m h , i = 1, . . . , 8 in (51). Clearly, Here we have used again that τ m ≤ τ ≤ 2 . In a similar way we obtain Inserting the above estimates into (51) we obtain In order to estimate the first term on the right hand side we write e 0 h = (I 0 h u e 0 − u e 0 ) + (u e 0 − u 0 h ) and recall the definition (38) of u 0 h as an L 2 projection: The remainder of the proof follows from (52) and Lemma 7.
Using the result of Theorem 2 we can now also derive an error bound on the surface.

Corollary 1
In addition to the assumptions of Theorem 2 suppose that there exists Proof Let us fix m ∈ {1, . . . , M} and define T m We infer from (8) and (36) and therefore We now argue in a similar way as in [6], page 368. Using an interpolation inequality and an inverse estimate we infer that where the last inequality follows from (54), (8) and the assumption that h T ≥ α , T ∈ T m Γ,h . In a similar way we obtain Therefore, it is not difficult to verify with the help of Taylor expansion that Combining this relation with (58) in the "Appendix" we find that In particular, we infer from (57) that replacingp byp in the extension of v, f and u 0 will not affect the result of Theorem 2. In contrast, it is not straightforward to handle the interpolation terms I m h ρ m and I m−1 h ρ m−1 in (55). Applying a standard interpolation estimate to ρ m − I m h ρ m will result in a term of the form h 2 ρ m H 2 ≈ h 2 2 , which we are currently not able to analyze. The results of our test calculations below however show that the use of the interpolation operator in (55), (56) does not lead to reduced convergence rates. More precisely we investigate the experimental order of convergence (eoc) for the following errors: whereû m (x) = u(p(x, t m ), t m ). We use the finite element toolbox Alberta 2.0, [24], and implement a similar mesh refinement strategy to that in [2] with a fine mesh constructed in D m h and a coarser mesh in Ω\D m h . The linear systems appearing in each time step were solved using GMRES together with diagonal preconditioning. The values of h given below are such that h := max T ∈D m h h T , h T = diam(T ). Remark 4 Although the analysis requires γ > 0, the method works with γ = 0 and produces very similar eocs to the ones displayed in the tables below for γ = 0.01.

2D examples
We set Ω = (−2.4, 2.4) 2 , T = 0.1, and choose γ = 0.01, = 85.33 h as well as a uniform time step τ m = 0.0025ε 2 , m = 1, . . . , M. In all our examples below Γ (t) will be a circle Γ (t) = {x ∈ R 2 | |x − m(t)| = 1} of radius 1 with center m(t) ∈ R 2 . In addition to E 1 , E 2 we shall also investigate the errors appearing in Corollary 1. To do so we choose L > 0 and define the following quadrature points In our computations L = 200 turned out to be sufficient.
Example 1 For our first example we consider the stationary unit circle Γ (t) = Γ = S 1 , t ∈ [0, T ] described as the zero level set of the function φ(x) is a solution of (1), (2) for the velocity field v(x) = π 2 (x 2 , −x 1 ) T , f = 0 and the initial data u 0 (x) = x 1 x 2 . A similar choice of velocity appears in Example 3 in [10]. In Tables 1  and 2 we display the values of E i , i = 1 → 4, together with the eocs.  Tables 3 and 4 where we see eocs that are similar to the ones in Tables 1 and 2. We see that the eoc for E 1 is reducing towards 4, the eocs for E 2 and E 3 are close to 4 and the eoc for E 4 is between 2 and 3 which is better than Theorem 2 predicts. Since E 1 and E 3 approximate L 2 -errors, higher eocs can be expected although a corresponding proof is by no means straightforward and beyond the scope of this paper. The higher eoc for E 2 presumably reflects a superconvergence effect because we consider ∇(I m hû m − u m h ) rather than ∇(û m − u m h ). We expect that E 4 will tend towards 2 if ε, h and τ are reduced further.

3D example
Example 3 Here we consider the first example in Section 7 of [18] in which a family of expanding and collapsing spheres is considered such that Γ (t) = {x ∈ R 2 | |x| = r (t)} where r (t) = 1 + sin 2 (π t), described as the zero level set of φ(x, t) = x 2 1 + x 2 2 + x 2 3 − r (t) 2 . The function u : S T → R, u(x, t) =   For the choice L = 200 the results are displayed in Table 5, where we see eocs close to 4 for E 3 and eocs close to 2 for E 4 . We conclude with Fig. 1 in which we present the approximate u m h at times t m = 0, 0.2, 0.4 plotted on the zero level surface of φ m h .