{'fraction_non_alphanumeric': 0.052881686365297914, 'fraction_numerical': 0.020519069087364446, 'mean_word_length': 4.014786753024564, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A decoupled, stable, and linear FEM for a phase-field model of variable density two-phase incompressible surface flow

The paper considers a thermodynamically consistent phase-field model of a two-phase flow of incompressible viscous fluids. The model allows for a non-linear dependence of fluid density on the phase-field order parameter. Driven by applications in biomembrane studies, the model is written for tangential flows of fluids constrained to a surface and consists of (surface) Navier-Stokes-Cahn-Hilliard type equations. We apply an unfitted finite element method to discretize the system and introduce a fully discrete time-stepping scheme with the following properties: (i) the scheme decouples the fluid and phase-field equation solvers at each time step, (ii) the resulting two algebraic systems are linear, and (iii) the numerical solution satisfies the same stability bound as the solution of the original system under some restrictions on the discretization parameters. Numerical examples are provided to demonstrate the stability, accuracy, and overall efficiency of the approach. Our computational study of several two-phase surface flows reveals some interesting dependencies of flow statistics on the geometry.


Introduction
Due to many applications in the physical and biological sciences, binary (two-phase) flows have attracted increasing interest in the past decades. The phase-field method has emerged as a powerful and effective technique for modeling multi-phase problems in general. This method describes the physical system using a set of field variables that are continuous across the interfacial regions separating the phases. The main reasons for the success of the phase-field methodology are the following: it replaces sharp interfaces by thin transition regions called diffuse interfaces, making front-tracking unnecessary, and it is based on rigorous mathematics. The classical diffuse-interface description of binary incompressible fluid with density-matched components is the so-called Model H [22], which is a Navier-Stokes-Cahn-Hilliard (NSCH) system. Several generalization to the case of variable density components have been presented [2,4,6,7,8,14,15,29,42]. Some of these models relax the incompressibility constraint to quasi-compressibility and not all models satisfy Galilean invariance, local mass conservation or thermodynamic consistency. In this paper, we focus on the generalization of Model H first presented in [2]. This model is thermodynamically consistent when the density of the mixture depends linearly on the phase-field variable (concentration or volume/surface fraction). We present an extension of the model in [2] that is thermodynamically consistent for a general monotone relation of density and phase-field variable.
Driven by applications in the analysis of lateral organization of plasma membranes (see [45,46] for more context), we adopt this new model to simulate two-phase flow dynamics on arbitrary-shaped closed smooth surfaces. While there exist many computational studies of multi-component fluid flows in planar and volumetric domains (see, e.g., [19,28,31,44,47] for some recent works), the number of papers where NSCH systems are treated on manifolds is limited. This can be explained by the fact that solving equations numerically on general surfaces poses additional difficulties that are related to the discretization of tangential differential operators and the approximate recovery of (possibly) complex shapes. In [30], the authors resort to a stream function formulation, decouple the surface Navier-Stokes problem from the surface Cahn-Hilliad problem, and use a parametric finite element approach. Yang and Kim [43] experimented with a finite difference method on staggered marker-and-cell meshes for a NSCH system embedded in the mesh.
In the present paper, we study a geometrically unfitted finite element method for the simulation of two-phase incompressible flow on surfaces. Our approach builds on earlier work on a unfitted FEM for elliptic PDEs posed on surfaces [36] called TraceFEM. Unlike some other geometrically unfitted methods for surface PDEs, TraceFEM employs sharp surface representation. The surface can be defined implicitly and no knowledge of the surface parametrization is required. This method allows to solve for a scalar quantity or a vector field on a surface, for which a parametrization or triangulation is not required. TaceFEM has been extended to the surface Stokes problem in [32,33] and the surface Cahn-Hilliard problem in [45]. One of the main advantages of TraceFEM is that it works well for PDEs posed on evolving surfaces, including cases with topological changes [27]. Although we will not benefit from this advantage in this paper (as the surfaces are fixed), it will become handy in a future modeling of membrane deformation and fusion for multicomponent lipid bilayers.
Although technical, the numerical analysis of diffuse-interface models for two-phase fluids with matching densities can be largely built on well-established analyses for the incompressible Navier-Stokes equations and Cahn-Hilliad equations alone. See, for example, [12] for the convergence study of a coupled implicit FE method, [25] for the analysis of a decoupled FE scheme, and [21] for the analysis of a linear and second-order in time FE method. The design of energy stable, efficient, and consistent discretizations of diffuse-interface models for two-phase fluids with non-matching densities turns out to be more challenging.
Following the ideas from [18] for the incompressible Navier-Stokes equations with variable density, in [40] Shen and Yang suggested a stable time-stepping method for a NSCH system with the momentum equation modified based on an inconsistent mass conservation equation. In [41], the same authors introduced a time-stepping scheme for the thermodynamically consistent model from [2]. This scheme (which is first order and linear, and it decouples fluid and phase equations) introduces some extra terms in the momentum equation and modifies the advection flux for the order parameter to achieve energy stability. In the same paper, the authors present the stability analysis for a semi-discrete case. In [48], the approach from [41] was extended with the application of a scalar auxiliary variable method [38] and the stability was studied for a finite difference scheme on uniform staggered grids. A fully discrete FE method for the model from [2] was analyzed in [17], subject to the coupled treatment of the phase and fluid equations and a modification of the momentum equation. In [13], this analysis is extended to a posteriori estimates and adaptivity. The question of building provably stable and computationally efficient methods for thermodynamically consistent discretization models of two-phase fluids with variable densities remains largely open.
The present study contributes to the field with a first order in time, linear, and decoupled FE scheme for our generalization of the model from [2]. One advantage of our scheme is that it neither modifies the momentum equation nor alters the advection velocity. We prove that our scheme is stable under a mild time-step restriction. Moreover, the analysis tracks the dependence of the estimates on , the critical model parameter that defines the width of the transition region between phases. The fact that the model is posed on arbitrary-shaped closed smooth surfaces and the use of an unfitted finite element method for spatial discretization add some further technical difficulties to our analysis. However, apart from these extra technical details, the general line of arguments can be followed for the simplified to Euclidian setting or extended to other spatial discretization techniques.
Numerical analysis is even less developed for quasi-incompressible diffuse-interface models due to their additional complexity, e.g. a more subtle nonlinearity. Here, we would like to mention the analysis of a semi-discrete method (coupled, but resulting in linear systems at each time step) in [42].
The rest of the paper is organized as follows. In Sec. 2, we state our generalization of the model from [2] and its variational formulations. A decoupled linear FEM based on the application of TraceFEM to our model is presented in Sec. 3, while the analysis is reported in Sec. 4. In Sec. 5, we report several numerical results, including a convergence test and the simulation of the Kevin-Helmholtz and Rayleigh-Taylor instabilities. Sec. 6 collects a few concluding remarks.

Problem definition
On surface Γ we consider a heterogeneous mixture of two species with surface fractions c i = S i /S, i = 1, 2, where S i are the surface area occupied by the components and S is the surface area of Γ. Since S = S 1 +S 2 , we have c 1 + c 2 = 1. Let c 1 be the representative surface fraction, i.e c = c 1 . Moreover, let m i be the mass of component i and m is the total mass. Notice that density of the mixture can be expressed as ρ = m S = m1 S1 S1 S + m2 S2 S2 S . Thus, Similarly, for the dynamic viscosity of the mixture we can write where η 1 and η 2 are the dynamic viscosities of the two species.
In order to state the model posed on surface Γ, we need some preliminary notation. Let P = P(x) := I − n(x)n(x) T for x ∈ Γ be the orthogonal projection onto the tangent plane. For a scalar function p : Γ → R or a vector function u : Γ → R 3 we define p e : O(Γ) → R, u e : O(Γ) → R 3 , extensions of p and u from Γ to its neighborhood O(Γ). The surface gradient and covariant derivatives on Γ are then defined as ∇ Γ p = P∇p e and ∇ Γ u := P∇u e P. These definitions are independent of a particular smooth extension of p and u off Γ. On Γ we consider the surface rate-of-strain tensor [20] given by The surface divergence operators for a vector g : Γ → R 3 and a tensor A : Γ → R 3×3 are defined as: with e i the ith standard basis vector in R 3 . The classical phase-field model to describe the flow of two immiscible, incompressible, and Newtonian fluids is the so-called Model H [22]. One of the fundamental assumptions for Model H is that the densities of both components are matching. Several extensions have been proposed to account for the case of non-matching densities; see Sec. 1. Here, we restrict our attention to a thermodynamically consistent generalization of Model H first presented in [2]. For surface based quantities the model reads: Here, u is the surface averaged tangential velocity u = cu 1 + (1 − c)u 2 , density is given by (2.1) and viscosity by (2.2), σ γ is line tension and µ is the chemical potential defined in (2.7). Force vector f , with f · n = 0, is given. In (2.6) M is the mobility coefficient, while in (2.7) f 0 (c) is the specific free energy of a homogeneous phase and is a parameter related to the thickness of the interface between the two phases. We recall that in order to have phase separation, f 0 must be a non-convex function of c. For a comprehensive discussion of the related Navier-Stokes and Cahn-Hilliard equations on surfaces we refer, for example, to [23,10,9]. Problem (2.4)-(2.7) has only one additional term with respect to Model H: the middle term at the righthand side in eq. (2.4). Notice that this term vanishes in the case of matching densities since dρ dc = (ρ 1 − ρ 2 ). However, the term is crucial for thermodynamic consistency when the densities do not match and it can be interpreted as additional momentum flux due to diffusion of the components driven by the gradient of the chemical potential.
In practice, we are interested in more general relations than (2.1) between ρ and c, since depending on the choice of f 0 (c) (and because of numerical errors while computing) the order parameter may not be constrained in [0, 1] and so ρ and η based on (2.1) may take physically meaningless (even negative) values. Since we do not see how to show the thermodynamic consistency of (2.4)-(2.7) for non-linear ρ(c), we propose a further modification. Without the loss of generality let ρ 1 ≥ ρ 2 . For a general dependence of ρ on c, it is reasonable to assume that ρ is a smooth monotonic function of c, i.e. dρ dc ≥ 0 (for ρ 1 ≥ ρ 2 ), and so we can set dρ dc = θ 2 . (2.8) Then, we replace (2.4) with The updated model (2.9),(2.5)-(2.7) obviously coincides with (2.4)-(2.7) for ρ(c) from (2.1), but exhibits thermodynamic consistency for a general monotone ρ-c relation as we show below. The consistency is preserved if M is a non-negative function of c rather than a constant coefficient. Thermodynamic consistent extensions of (2.4)-(2.7) for a generic smooth ρ(c) (no monotonicity assumption) were also considered in [1,3] for the purpose of well-posedness analysis. Those extensions introduce yet more term(s) in the momentum equation, so for computational needs and numerical analysis purpose we opt for (2.9). From now on, we will focus on problem (2.9), (2.5)-(2.7). This is a Navier-Stokes-Cahn-Hilliard (NSCH) system that needs to be supplemented with the definitions of mobility M and free energy per unit area f 0 . Following many of the existing analytic studies, as well as numerical studies, we assume mobility to be constant (i.e., independent of c). As for f 0 , a classical choice is the Ginzburg-Landau double-well potential where ξ defines the barrier height, i.e. the local maximum at c = 1/2 [11]. We set ξ = 1 for the rest of the paper.
For the numerical method, we need a weak (integral) formulation. We define the spaces We define the Hilbert space V 0 T as an orthogonal complement of E in V T (hence V 0 T ∼ V T /E), and recall the surface Korn's inequality [23]: In (2.12) and later we use short notation · for the L 2 -norm on Γ. Finally, we define L 2 0 (Γ) := { p ∈ L 2 (Γ) | Γ p ds = 0 }. To devise the weak formulation, one multiplies eq. (2.9) by v ∈ V T , eq. (2.5) by q ∈ L 2 (Γ), eq. (2.6) by v ∈ H 1 (Γ), and eq. (2.7) by g ∈ H 1 (Γ) and integrates all the equations over Γ. For eq. (??) and (2.6), one employs the integration by parts identity, which for a closed smooth surface Γ reads: here κ is the sum of principle curvatures. Identity (2.13) is applied to the second term in (2.6) (i.e. g = cu), which leads to no contribution from the curvature term since u is tangential, and to the third term in (2.6) (i.e. g = M ∇ Γ µ), which makes the curvature term vanish as well. For a similar reason (component-wise), the curvature term vanishes also when identity (2.13) is applied to the diffusion term in (2.9). The weak (integral) formulation of the surface NSCH problem (2.9), (2.5)-(2.7) reads: Find (u, p, c, µ) ∈ Remark 2.1 By testing (2.14) with v = u, (2.15) with q = p, (2.16) with v = µ, and using (2.7), one readily obtains the following energy equality: In other words, model (2.9), (2.5)-(2.7) is thermodynamically consistent, i.e. the system is dissipative for Indeed, the first two terms in eq. (2.14) tested with v = u can be handled as follows: Using generic ρ = ρ(c), (2.8), and (2.5) we get With the help of (2.6) , this yields which will cancel with the second term at the right-hand side in eq. (2.14) tested with v = u. Thus, from eq. (2.14) tested with v = u we obtain the following equality: where we have also used eq. (2.15) tested with q = p. From eq. (2.16) tested with v = µ and multiplied by σ γ , we get: The first term on the right-hand side can be handled using (2.7) as follows: Plugging this into (2.20) and then using what we get in (2.19), we obtain the energy balance (2.18).

Numerical method
For the discretization of the variational problem (2.14)-(2.17) we apply an unfitted finite element method called trace finite element approach (TraceFEM). To discretize surface equations, TraceFEM relies on a tessellation of a 3D bulk computational domain Ω (Γ ⊂ Ω holds) into shape regular tetrahedra untangled to the position of Γ. We start with required definitions to set up finite element spaces and variational form. A few auxiliary results will be proved. Then we proceed to the fully-discrete method by introducing a splitting time discretization, which decouples (2.14)-(2.17) into one linear fluid problem and one phase-field problem per every time step.
Surface Γ is defined implicitly as the zero level set of a sufficiently smooth (at least Lipschitz continuous) function φ, i.e. Γ = {x ∈ Ω : φ(x) = 0}, such that |∇φ| ≥ c 0 > 0 in a 3D neighborhood U (Γ) of the surface. The vector field n = ∇φ/|∇φ| is normal on Γ and defines quasi-normal directions in U (Γ). Let T h be the collection of all tetrahedra such that Ω = ∪ T ∈T h T . The subset of tetrahedra that have a nonzero intersection with Γ is denoted by T Γ h . The grid can be refined towards Γ, however the tetrahedra from T Γ h form a quasi-uniform tessellation with the characteristic tetrahedra size h. The domain formed by all tetrahedra in T Γ h is denoted by Ω Γ h . On T Γ h we use a standard finite element space of continuous functions that are piecewise-polynomials of degree k. This bulk (volumetric) finite element space is denoted by V k h : In the trace finite element method formulated below, the traces of functions from V k h on Γ are used to approximate the surface fraction and the chemical potential. Our bulk velocity and pressure finite element spaces are either Taylor-Hood elements on Ω Γ h , or equal order velocity-pressure elements These spaces are employed to discretize the surface Navier-Stokes system. Surface velocity and pressure will be represented by the traces of functions from V h and Q h on Γ. In general, approximation orders k (for the phase-field problem) and m (for the fluid problem) can be chosen to be different.
Assumption 3.1 We assume that integrals over Γ can be computed exactly, i.e. we do not consider geometry errors.
In practice, Γ has to be approximated by a (sufficiently accurate) "discrete" surface Γ h ≈ Γ in such a way that integrals over Γ h can be computed accurately and efficiently. For first order finite elements (m = 1, k = 1), a straightforward polygonal approximation of Γ ensures that the geometric approximation error is consistent with the finite element interpolation error; see, e.g., [36]. For higher order elements, numerical approximation of Γ based on, e.g., isoparametric trace FE can be used to recover the optimal accuracy [16] in a practical situation when parametrization of Γ is not available explicitly.
There are two well-known issues related to the fact that we are dealing with surface and unfitted finite elements: 1. The numerical treatment of condition u · n. Enforcing the condition u h · n = 0 on Γ for polynomial functions u h ∈ V h is inconvenient and may lead to locking (i.e., only u h = 0 satisfies it). Instead, we add a penalty term to the weak formulation to enforce the tangential constraint weakly.
2. Possible small cuts of tetrahedra from T Γ h by the surface. For the standard choice of finite element basis functions, this may lead to poorly conditioned algebraic systems. We recover algebraic stability by adding certain volumetric terms to the finite element formulation.
To make the presentation of the finite element formulation more compact, we introduce the following finite element bilinear forms related to the Cahn-Hilliard part of the problem: The second term in (3.3) and (3.4) take care of the issue of the small element cuts (the above issue 2) [16]. These terms are consistent up to geometric errors related to the approximation of Γ by Γ h and n by n h in the following sense: any smooth solution µ and c can be extended off the surface along (quasi)-normal directions so that n · ∇µ = 0 and n · ∇c = 0 in Ω Γ h . We set the stabilization parameters as follows: For the stability of a numerical method, it is crucial that the computed density and viscosity stay positive. The polynomial double-well potential does not enforce c to stay within [0, 1] interval and hence ρ and η may eventually take negative values, if one adopts the linear relation between c and ρ in (2.1) or between c and η in (2.2). Numerical errors may be another reason for the order parameter to depart significantly from [0, 1]. Therefore, assuming without loss of generality that ρ 1 ≥ ρ 2 and η 1 ≥ η 2 , we first replace (2.1) and (2.2) with the following cut-off functions that ensure density and viscosity stay positive: Note that unlike some previous studies, we clip the linear dependence (2.1) and (2.2) only from below. The resulting convexity of ρ(c) plays a role in the stability analysis later. Nevertheless, (3.5) is not completely satisfactory since θ 2 = δρ δc is discontinuous, while we need θ to be from C 1 . To this end, we approximate ρ(c) from (2.1) by a smooth monotone convex and uniformly positive function. In our implementation we let θ 2 = ρ1−ρ2 2 (tanh(c/α) + 1), with α = 0.1, and ρ(c) = c 0 θ 2 (t)dt + ρ 2 . Later we make use of the decomposition of a vector field on Γ into its tangential and normal components: u = u + (u · n)n. For the Navier-Stokes part, we introduce the following forms: In (3.6), τ > 0 is a penalty parameter to enforce the tangential constraint (i.e., to address the above issue 1), while β u ≥ 0 in (3.6) and β p ≥ 0 in (3.9) are stabilization parameters to deal with possible small cuts. They are set according to [24]: If one uses equal order pressure-velocity trace elements instead of Taylor-Hood elements, then form (3.9) should be replaced by The second term in (3.7) is consistent since the divergence of the true tangential velocity is zero. To avoid differentiation of the projector operator, one may use the identity ∇ Γ u = ∇ Γ u − (u · n)H to implement the a-form and c-form. For the analysis, we need the identity for the form (3.7) given in the following elementary lemma.
Proof: Using the definition of the covariant gradient in terms of tangential operators, ∇ Γ u = P(∇u)P, and the integration by parts (2.13), we compute for the first integral term in (3.7) From this equality and using ∇ Γ ρ = dρ dc ∇ Γ c we obtain: We recall that ρ = ρ − dρ d c c and substitute (3.11) into the definition of the form (3.7). After straightforward computations, we arrive at the result: For the numerical experiments in Sec. 5, we also add to bilinear form (3.6) the grad-div stabilization term [35], γ Γ div Γ u div Γ v ds. This term is not essential, in particular for the analysis in Sec. 4, but we find it beneficial for the performance of the iterative algebraic solver and for the overall accuracy of the solution. We set the grad-div stabilization parameter γ = 1.
At time instance t n = n∆t, with time step ∆t = T N , ζ n denotes the approximation of generic variable ζ(t n , x). Further, we introduce the following notation for a first order approximation of the time derivative: The decoupled linear finite element method analyzed and tested in this paper reads as follows. At time step t n+1 , perform: for At each sub-step of the scheme, a linear problem (Chan-Hilliard type system at step 1 and linearized Navier-Stokes system at step 2) has to be solved. This allows us to achieve low computational costs, while the scheme is provably stable under relatively mild restrictions. Moreover, the results of numerical experiments in Sec. 5 do not show that any restrictions on the discretization parameters are required in practice (Remark 4.3 discusses what arguments in our analysis require these restrictions).
Before proceeding with analysis, we note that the above scheme does not modify the advection velocity in (3.13) for the purpose of analysis, unlike some other linear decoupled schemes for the NSCH equations found in the literature. We also avoid another common helpful trick to prove energy stability of the variable density NSCH, namely the modification of the momentum equation based on a mass conservation condition of the form ∂ t ρ+u·∇ Γ c− div Γ (M ∇ Γ µ) = 0, which follow from (2.6) and (2.1). This modification, however, is not consistent if a non-linear relation between ρ and c is adopted; also it introduces several extra terms in the finite element formulation.

Analysis of the decoupled Finite Element method
For the analysis in this section, we assume no forcing term, i.e. f n+1 = 0 for all n. To avoid technical complications related to handling Killing vector fields (see, e.g. discussion in [5]), we shall also assume that the surface does not support any tangential rigid motions, and so V 0 T = V T . We further use the a b notation if inequality a ≤ c b holds between quantities a and b with a constant c independent of h, ∆t, , and the position of Γ in the background mesh. We give similar meaning to a b, and a b means that both a b and a b hold.
The following lemma will be useful in the proof of the main result, which is reported in Theorem 4.2.
For v h satisfying Γ v h ds = 0, the factor v h H 1 (Γ) on the r.h.s. can be replaced by ∇ Γ v h L 2 (Γ) .
Proof: Denote by p(x) ∈ Γ, x ∈ Ω, the closest point projection on Γ. Since Γ is smooth, there is a tubular neighborhood of Γ such that x = p(x) + sn, x ∈ U , defines the local coordinate system (s, y), with y = p(x) and |s| = dist(x, Γ). We can always assume Denote by Ω Γ h a "reachable from Γ" part of Ω Γ h in the following sense: for any x ∈ Ω Γ h the interval (x, p(x)) is completely inside Ω Γ h . Let g ± (y) = ± sup{s ∈ R : (y ± sn, y) ⊂ Ω Γ h }, y ∈ Γ, where g ± (y) are piecewise smooth and g ± L ∞ (Γ) h. Thanks to the co-area formula and the smoothness of Γ, it holds Substituting this in (4.4) and using Ω Γ h ⊂ Ω Γ h we get Letting u = u h ∈ V k h and applying FE inverse inequality to treat the last term, we arrive at (4.5) Next, we need the following technical result from Lemma 7.9 in [16]: There is δ h, depending only on the shape regularity of the mesh T Γ h such that for any T ∈ T Γ h there exists a ball B δ (T ) ⊂ T ∩ Ω Γ h of radius δ. Since on every tetrahedron u h is polynomial function of a fixed degree by the standard norm equivalence argument, we have with C depending only on δ and the shape regularity of T . Raising both parts of this inequality to power p, summing over all T ∈ T Γ h , raising both parts to power 1/p and using B δ (T ) ⊂ Ω Γ h we get We now use this in (4.5) and apply another FE inverse inequality to get (4.6) For u ∈ H 1 (Γ), recall that the Sobolev embedding theorem implies u ∈ L p (Γ), p ∈ [1, ∞), and This result follows from the corresponding embedding theorem in R 2 by standard arguments based on the surface local parametrization and partition of unity. Using this in (4.6), we obtain the estimate (4.7) Letting p = | ln h| proves the result in (4.1). Applying the Poincaré inequality for the function u h satisfying Γ u h ds = 0 proves the second claim of the lemma. Following [39], we modify the double-well potential in (2.10) for the purpose of analysis so that it is C 2 smooth but has quadratic growth for large |c|. Straightforward calculations give us the following expression for f 0 (c) with sufficiently large but fixed α > 1: Function f 0 (x) satisfies the following Lipschitz condition with L = 3α 2 −1 and growth condition |f 0 (x)| ≤ L|x|.
∆t (a(η n ; u n h , u n h ) + a µ (µ n , µ n ) + s h (p n h , p n h )) ≤ K, (4.10) Proof: We use induction in N to prove (4.10). For N = 0, the estimate is trivial and provided by the initial condition. For the induction step, assume that it holds with N = n.
in (3.14) and adding the two equations together, we get With the truncated Taylor expansion f 0 (c n+1 n+1 t | 2 f 0 (ξ n ), and (4.8) we get: Since the second term at the right-hand side has a negative sign, we let γ c be large enough in order to cancel it with the 1 2 of the last term on the left-hand side of (4.11). We obtain After re-arranging terms, multiplying by ∆t, and dropping some non-negative terms on the left-hand side, we obtain the following inequality The first two terms on the right-hand side of (4.13) are bounded due to the induction assumption. To handle the third term on the right-hand side, we let c 0 = |Γ| −1 Γ c n+1 h ds and use Lemma 4.1 that yields 2 by the definition of the a c form in (3.4). This, the Cauchy inequality, and the equality u n h · ∇ Γ µ n+1 h = u n h · ∇ Γ µ n+1 h help us with the following estimate: (4.14) In last inequality in (4.14), we used the induction assumption to estimate u n h and the fact that ρ n ≥ ρ 2 by definition of the ρ(c) function. Letting v h = 1 in (3.13) we have We conclude that |c 0 | in (4.14) can be bounded by a constant depending only on the initial data. Then, using Young's inequality in (4.14) we get the following bound with some constants C 1 , C 2 independent of h, ∆t, , and position of Γ in the background mesh. Using this back in (4.13) with ∆t satisfying assumptions of the theorem, and applying (4.10) for N = n for the remainder terms on the right-hand side of (4.13), we get: with some C > 0 independent of h, ∆t, and the position of Γ in the mesh. We need this bound later in the proof.
Let v h = u n+1 h in (3.15) and q h = −p h in (3.16), and add the two equations together. We also apply Lemma 3.1 to handle the c-form. This brings us to the following equality: By adding ± 1 2∆t (ρ n+1 ) 1 2 u n+1 h 2 and re-arranging terms, we get Denote by P h : H 1 (Γ) → V k h an H 1 -type projection into V k h given by the a µ bilinear form: By standard analysis of the TraceFEM for the Laplace-Beltrami problem, e.g. [37], we have Since ρ is a convex function of c and dρ dc ≥ 0, we have Using (4.18) and eq. (3.13) with v h = P h |θ n+1 u n+1 h | 2 , for the terms on the second line of (4.16) we obtain with Next, we use (4.19) to simplify (4.16) as follows: After adding (4.12) multiplied by σ γ to (4.20) and dropping some non-negative terms on the left-hand side, we arrive at It remains to estimate the terms on the right-hand side of (4.21). We handle the first term by invoking the result of Lemma 4.1 as follows: With the help of ∆t ≤ C| ln h| −1 for sufficiently small C and 1 we obtain Now, we proceed with the terms in I h . For the first terms in I h , the Cauchy-Schwarz inequality and estimate (4.17) for the L 2 (Γ) norm of e h gives: To estimate the right-most factor, we need the inequalities which follow by applying Lemma 4.1 componentwise and then using the Korn inequality (2.12) and the fact that η n+1 is uniformly bounded from below (see the definition in (3.5)). Thanks to 0 ≤ θ 2 1 and 0 ≤ dθ 2 dc = d 2 ρ dc 2 1, we have Using this together with the a priori bound for [c h ] n+1 t from (4.15) in (4.25) and using the assumption h ≤ c| ln h| −1 ∆t for sufficiently small c, we have Each term on the right hand side can be treated individually by invoking Lemma 4.1, a priori bound from (4.15), and induction assumption: Using these estimates back into (4.27) and with the assumption h ≤ c| ln h| − 3 2 |∆t| − 1 2 for sufficiently small c, we arrive at the estimate Finally, we substitute (4.24) and (4.28) in (4.21) to obtain This completes the induction step.

Remark 4.3
From the proof we see that the assumption h ≤ c| ln h| −1 min{∆t, | ln h| − 1 2 |∆t| 1 2 } results from the fact that at the discrete level we cannot test the transport equation for the order parameter with v h = |u n+1 h | 2 , and so we have to project |u n+1 h | 2 (or |u n+1 h | 2 for surfaces) in the finite element space V k h and handle the resulting inconsistency. If the finite element velocity space is such that |v h | 2 ∈ V k h for v h ∈ V h , then the upper bound on h is not needed in the analysis (in the surface case we still would need to handle |u n+1 h | 2 − |u n+1 h | 2 , but this is possible due to the control over n · u n+1 h 2 that we have thanks to penalty term in the TraceFEM formulation). An example, when |v h | 2 ∈ V k h holds, is P 1 -P 1 stabilized velocity-pressure element combined with P 2 element for the order parameter and chemical potential.

Numerical results
After checking the convergence orders of the method described in Sec. 3, we present a series of numerical tests to study well-known two-phase fluid flows (the Kelvin-Helmholtz and Rayleigh-Taylor instabilities) on surfaces. Thanks to our method, we can investigate the effect of line tension on such instabilities. In addition, for the Rayleigh-Taylor instability we assess the effect of viscosity and surface shape.
For all the tests, we choose P 2 -P 1 finite elements for fluid velocity and pressure and P 1 -P 1 finite elements for surface fraction and chemical potential.

Convergence test
We proceed with checking the spatial accuracy of the finite element method described in Sec. 3 with a benchmark test. The aim is to validate our implementation of the method. For this purpose, we consider the two-phase fluid system on the unit sphere centered at the origin. The surface is characterized as the zero level set of function φ(x) = x 2 − 1 and is embedded in an outer cubic domain Ω = [−5/3, 5/3] 3 . We choose Van der Waals "tanh" exact solution for the surface fraction and solenoidal Killing vector field for velocity: Nonzero CH equation forcing term is computed from (2.6). We set = 0.05. Fig. 1 (leftmost panel) displays c(0, x). For this test, we consider fluids with matching densities and viscosities: ρ 1 = ρ 2 = 1 and η 1 = η 2 = 1. In addition, we set M = 0.05 and σ γ = 0. The time interval of interest is t ∈ [0, 1], during which the initial configuration c 0 is rotated by 180 • . See Fig. 1. Notice that by setting σ γ = 0 the NSCH system one-way coupled: phase-separation is affected by the fluid flow, but not vice versa. The initial triangulation T h of Ω consists of eight sub-cubes, where each of the sub-cubes is further subdivided into six tetrahedra. Further, the mesh is refined towards the surface, and ∈ N denotes the level of refinement, with the associated mesh size h = 10/3 2 +1 . For the purpose of numerical integration, we apply several "virtual" levels of refinements for the tetrahedra cut by the mesh and integrate our bilinear forms over a piecewise planar approximation of Γ on this virtual grid. This allowed us to apply a standard quadrature rule and reduce the geometric error in our convergence test. The time step was refined together with the mesh size according to ∆t = 1/(25 · 2 −2 ). For this test, we used BDF2 for the time discretization of the Cahn-Hilliard problem at Step 1, instead of BDF1 as reported in (3.13). This is why the time step is refined linearly. Table 1 reports the H 1 and L 2 errors for the velocity and the L 2 error for the order parameter at the end of the time interval, i.e. t = 1, for levels = 3, 4, 5. For each mesh, Table 1 gives the number of sublevels (virtual levels) used for more accurate numerical integration. We observe slightly better than expected convergence rates, which might be the effect of the interplay between interpolation and geometric error reduction.  Table 1: H 1 and L 2 errors for the velocity and the L 2 error for surface fraction at t = 1 for levels = 3, 4, 5 and rate of convergence.

The Kelvin-Helmholtz instability
While the Kelvin-Helmholtz (KH) instability is a classical example of two-phase fluid flow in planar or volumetric domains, the number of numerical studies on curved surfaces is limited. The KH instability arises when there is a difference in velocity at the interface between the two fluids and a perturbation is added to the interface. This perturbation eventually makes the interface curl up and generates a vortex strip. Here, we will simulate the KH instability on a sphere and investigate the effect of varying line tension.
To design this experiment, we follow what done in [26,24,34]. The initial velocity field is given by the counter-rotating upper and lower hemispheres with speed approximately equal 1 closer to equator and vanishing at the poles. The velocity field has a sharp transition layer along the equator, where the perturbation is added. See, e.g., [34] for details on the perturbation. The initial surface fraction is given by where = 0.01. Also for this test, we consider fluids with matching densities and viscosities: ρ 1 = ρ 2 = 1 and η 1 = η 2 = 10 −5 . In addition, we set M = 0.01. We consider time interval [0, 10]. We select mesh level = 6 (see mesh description for the previous test). We choose ∆t = 1/640. Fig. 2 and 3 show the evolution of surface fraction and vorticity for three different values of line tension: σ γ = 0, 0.01, 0.1. The evolution of both quantities does not vary significantly when going from σ γ = 0 to σ γ = 0.01, although some differences can be noticed from t = 4.5 on. Changing to σ γ = 0.1 produces more evident differences, starting already from t = 1.375. When σ γ = 0, the NS part of the CHNS system is two-way coupled to the CH part. So, the differences are significant both for surface fraction and vorticity.

The Rayleigh-Taylor instability
The Rayleigh-Taylor (RT) instability occurs when a gravity force is acting on a heavier fluid that lies above a lighter fluid. As the RT instability develops, "plumes" of the lighter fluid flow upwards (with respect to the gravitational field) and "spikes" of the heavier fluid fall downwards. We will simulate the RT instability on a sphere and on a torus with the aim of investigating the effect of the geometry. In addition, we will vary line tension and fluid viscosity. We take two fluids with densities ρ 2 = 3, ρ 1 = 1 and matching viscosities η 1 = η 2 = η, which will be specified for each test. The initial surface fraction is given by where = 0.025 and z rand is a uniformly generated random number from the range (−0.1 , 0.1 ). The role of the perturbation generated by z rand is to onset the RT instability. We set M = 0.0025. Let us start with the sphere. We select mesh level = 5 (see mesh description for the convergence test) and set ∆t = 0.1. Next, we consider an assymetric torus with constant distant from the center of the tube to the origin R = 1 and variable radius of the tube: r min = 0.3 ≤ r(x, y) ≤ r max = 0.6, with r(x, y) = r min + 0.5(r max − r min )(1 − x √ x 2 +y 2 ). We characterize the torus surface as the zero level set of function φ = (x 2 + y 2 + z 2 + R 2 − r(x, y) 2 ) 2 − 4R 2 (x 2 + y 2 ). The torus is embedded in an outer domain Ω = [−5/3, 5/3] 3 , just like the sphere. We also selected same mesh level (l = 5) and same time step (∆t = 0.1) as for the sphere. We set the line tension to σ = 0.025 and vary the viscosity: η = 10 −2 , 10 −1 , 1. Fig. 5 displays the evolution of the surface fractions for these three values of viscosity. First, we observe that in all cases the instability develops more slowly on the "skinny" side of the torus. See second column in Fig. 5. The fact that geometry has a considerable effect on the surface RT instability is also clear when one compares the results on the sphere and the torus for the same values of σ γ and η, i.e. the top row in Fig. 5 with the bottom row in Fig. 4. In particular, notice that while the heavier fluid reaches the bottom of the sphere around t = 30 (Fig. 4, bottom second-last panel), the two fluids are still very much mixed on the torus at t = 160 (Fig. 5, top left panel). We need to increase the viscosity value to 1 to be able to see most of the heavier fluid at the bottom of the torus at t = 160 (Fig. 5, bottom left panel), although that is still far from being settled.

Conclusions
In this paper, we presented an extension of a well-known phase field model for two-phase incompressible flow, and we applied and analyzed an unfitted finite element method for its numerical approximation. The advantage of our model is that it is thermodynamically consistent for a general monotone relation of density and phase-field variable. Because of our interest in biomembranes, this Navier-Stokes-Cahn-Hilliard type system is posed on an arbitrary-shaped closed smooth surface. Although in this paper we assumed the surfaces to be rigid, our long term goal is to study two-phase incompressible flow on evolving shapes. In fact, biological membranes exhibit shape changes that need to be accounted for in a realistic model. This need dictated our choice for the numerical approach, which is a versatile geometrically unfitted finite element method called TraceFEM. In order to reduce the computational cost, the discrete scheme we proposed decouples the fluid problem (a linearized Navier-Stokes type problem) from the phase-field problem (a Cahn-Hilliard type problem with constant mobility) at each time step. An attractive feature of our scheme is that the numerical solution satisfies the same stability bound as the solution of the original system under some restrictions on the discretization parameters.
We validated our implementation of the proposed numerical scheme with a benchmark problem and applied it to simulate well-known two-phase fluid flows: the Kelvin-Helmholtz and Rayleigh-Taylor insta-bilities. We investigated the effect of line tension on such instabilities. For the Rayleigh-Taylor instability, we also assessed the effect of viscosity and surface shape, which plays an important role in the evolution of the instability.