Phase-field model for Hele-Shaw flows with arbitrary viscosity contrast. I. Theoretical approach

We present a phase-field model for the dynamics of the interface between two inmiscible fluids with arbitrary viscosity contrast in a rectangular Hele-Shaw cell. With asymptotic matching techniques we check the model to yield the right Hele-Shaw equations in the sharp-interface limit and compute the corrections to these equations to first order in the interface thickness. We also compute the effect of such corrections on the linear dispersion relation of the planar interface. We discuss in detail the conditions on the interface thickness to control the accuracy and convergence of the phase-field model to the limiting Hele-Shaw dynamics. In particular, the convergence appears to be slower for high viscosity contrasts.


I. INTRODUCTION
The dynamics of morphologically unstable interfaces is a major problem in nonequilibrium physics from both fundamental and applied points of view. Relevant examples of those are dendritic growth, directional solidification, flows in porous media, flame propagation, electrodeposition or bacterial growth [1]. The so-called Saffman-Taylor problem has played a central role in this context because of its relative simplicity both experimentally and in its theoretical formulation [1,2]. It deals with the motion of the interface between two inmiscible fluids within a Hele-Shaw cell. Due to the highly nonlinear and nonlocal nature of the interfacial dynamics of such systems, analytical understanding is scarce and restricted to high viscosity contrast [3], so in general one relies mostly on numerical work [4][5][6][7][8][9].
From a mathematical point of view, such systems are referred to as moving boundary problems. In practice this implies that one has to keep track of the interface where boundary conditions are applied, and solve a (linear) problem in the bulk which determines in turn the motion of the boundary. This kind of problem has traditionally been addressed in terms of boundary integral methods which reduce the dynamics of the interface to integrodifferential equations. The numerical integration of these equations is quite involved though, particularly for long times, due to stiffness and numerical instability of the equations. In the case of Hele-Shaw flows, boundary integral methods have succesfully been applied [6][7][8][9], although quite sophisticated algorithms have usually been necessary [9].
Recently, the so-called phase-field equations have been proposed in the context of solidification problems as a different approach to the interface dynamics [10][11][12][13][14][15][16][17][18][19]. In the spirit of the well known time-dependent Ginzburg-Landau models [20], the method avoids the tracking of the interface by introducing an auxiliary field (analogous to an order parameter) which locates the interface and whose dynamics is coupled to the other physical fields through an appropriate set of partial differential equations. In this way, there is no boundary condition to explicitely apply at the interface and the whole system is treated as bulk.
This method introduces a mesoscale ǫ which is not present in the original macroscopic equations and gives a finite thickness to the interface. The equations are then chosen in such a way that the original bulk equations and boundary conditions are recovered in the ǫ → 0 limit. Therefore the phase-field equations for a given model are not intended to describe the true mesoscale physics of the system, and are then not unique. In fact there is considerable freedom in choosing a particular form of them, with criteria of either numerical efficiency and convergence [13] or other physical criteria such as thermodynamic consistency [14]. In any case, the nature of the phase-field approach is completely different from the sharp-interface models and therefore the actual numerical advantages and limitations of both are also quite distinct. This makes the two approaches complementary and competitive in different physical situations. A remarkable advantage of the phase-field approach is that it is much simpler to implement satisfactorily from a numerical point of view. On the other hand, the phase-field approach is usually more amenable to generalization, in the sense that it allows to introduce variations and new elements without any major modification of the numerical scheme, for instance in the treatment of fluctuations, liquid crystals [18] and other complex fluids [9]. Finally, the phase-field approach can handle very naturally situations where the sharp interface model is not appropriate, such as for instance topology changes like interface pinching leading to the breakup of bubbles.
In this paper we introduce a phase-field model for Hele-Shaw flows with arbitrary viscosity contrast (or Atwood ratio) c = µ 1 −µ 2 µ 1 +µ 2 . Although in the high contrast limit c = 1 the Hele-Shaw dynamics is quite analogous to the one-sided solidification problem (in the appropiate approximations [8]), the arbitrary viscosity contrast case has been shown to exhibit quite different dynamics than solidification problems, and has in fact opened some interesting questions, particularly concerning the sensitivity of finger competition to viscosity contrast [4][5][6][7]21] and the long time asymptotics of the low viscosity contrast limit [7].
The model presented here is inspired in the vortex-sheet formulation of the problem [4], in which the relevant dynamic variable in the bulk is the stream function. Similar ideas have previously been applied to describe physically diffuse interfaces in the context of steady state selection in thermal plumes [22]. Usually phase-field models are naturally suited to symmetric situations (two-sided models). The present case of Hele-Shaw flow is no exception and becomes most efficient for c = 0. Finite c can also be handled but the model becomes computationally inefficient in the limit c → 1, since this limit must be taken formally after the ǫ → 0 limit. A phase-field model for this one-sided case must differ essentially from the one presented here, such as in the spirit of Ref. [19].
The layout of the rest of the paper is as follows: in Sec. II A we recall the Hele-Shaw macroscopic equations in terms of the stream function, whereas in Sec. II B we present our phase-field equations. We then show in Sec. III that the phase-field equations reduce to the macroscopic ones in the sharp-interface limit. Deviations from that limiting behavior are derived from the phase-field equations themselves to first order in the interface thickness in Sec. IV, and their effect on the linear regime is computed in Sec. V. Finally, a brief summary is given in Sec. VI.

II. THE MODEL
We consider the general case of an interface with surface tension σ between two fluids with distinct viscosities (µ 1 , µ 2 ) and densities (ρ 1 , ρ 2 ) moving in a rectangular Hele-Shaw cell of width W (x-direction) and gap b (z-direction), under an effective gravity g ef f (negative y-direction) and with an injection velocity V ∞ (positive y-direction). Label 1 (2) corresponds to the upper (lower) fluid.

A. Macroscopic equations
Darcy's law is assumed to hold for each fluid, thus defining a certain velocity potential in each bulk, but not on the interface. In contrast, the bulk incompressibility and the continuity of normal velocities on the interface allow us to define its harmonic conjugate, the stream function ψ, even on the interface through u x = ∂ y ψ, u y = −∂ x ψ, where u x , u y are the x, y components of the fluid velocity field u. Then Darcy's law results in a Laplace equation for the stream function (potential flow) and a certain jump for the tangential fluid velocities on the interface, whose value takes into account the Gibbs-Thomson relationship. The fact that the stream function is continuous at the interface makes the use of this variable particularly convenient. The Hele-Shaw equations in stream function formulation [4] can be written in dimensionless form as: which constitute the Hele-Shaw equations. Here r is a coordinate normal to the interface and with origin on it, positive in fluid 1 (0 ± then means on the interface coming from each side), s is arclength along the interface and such that the unit vectors satisfyŝ ×r =x ×ŷ, the subscripts stand for partial derivatives except for v n (s), which is the normal velocity of the interface, positive towards fluid 1, and with κ(s) the interface curvature, positive for a bump into fluid 2. The dynamics are controlled by the two dimensionless parameters We will not be interested in negative values of B (stable configuration) nor c (mirror image interface of −c). So B is a dimensionless surface tension, and can be understood as the ratio between the capillary (stabilizing) force and the driving (destabilizing) force (injection+gravity), and c is the viscosity contrast, which is so far completely arbitrary: 0 ≤ c ≤ 1. This corresponds to having set ourselves in the frame moving with the fluid at infinity (or, equivalently, with the mean interface) and taken W as unit length and 12(µ 1 +µ 2 ) as unit velocity (see [4]). Note that Eqs. (2.1,2.2) can be written together as where δ(r) is the Dirac delta distribution and w ≡ẑ · ( ∇ × u) is the fluid vorticity, which is confined to the interface.

B. Phase-field equations
We put forward the following phase-field model for the above equations with θ being the phase field: andŝ(θ) ≡r(θ) ×ẑ, together with the boundary condition θ(y → ±∞) = ±1, (2.9) so that θ = +1(−1) corresponds to fluid 1 (2). γ(θ), κ(θ) are functionals which generalize the magnitudes defined above for the interface, now to any level-set of the phase-field. If we leave the two last terms aside, Eq. (2.8) is the Cahn-Hilliard equation for a nonconserved order parameter or model A (without noise) in the classification of Ref. [20] of time-dependent Ginzburg-Landau models. The field in this model is known to relax towards a kink solution of a certain width in a short time scale, and then to evolve to minimize the length of the effective interface according to Allen-Cahn law (i.e. with normal velocity proportional to the local curvature). The factor multiplying the laplacian has been choosen to be ǫ 2 for the kink width to be O(ǫ), so that ǫ can be considered the interface thickness, i.e., the small parameter in the asymptotic analysis that will be performed in next section. On the other hand, the ǫ 2 factor in the time derivative ensures that the relaxation towards the kink is much faster than the evolution of the interface. Notice that model A describes the relaxational dynamics of a non-conserved order parameter, whereas our problem is actually non-relaxational and strictly conserved (mass consevation and inmiscibility). The other two term in the phase-field equation will correct this apparent contradiction. In order to cancel out the local Allen-Cahn dynamics of the interface which is buit in model A, we add the term ǫ 2 κ(θ)| ∇θ|. It will be shown that such term cancels out Allen-Cahn law by giving rise, to leading order, to an identical contribution but with opposite sign. With these elements so far, our phase-field relaxes to a kink profile located along an arbitrary interface which (if sufficiently smooth) remains almost completely stationary, regardless of its shape. This is because the dynamical effect of surface tension associated to the Ginzburg-Landau free energy has been removed (up to first order) and the interface has not yet been coupled to the fluid flow, represented by the stream function. This coupling is achieved by adding the last term in Eq. (2.8), which stands for −ǫ 2 u · ∇θ and thus sets the phase-field -and therefore the interface-in the frame moving with the fluid velocity u. This term restores the fully nonlocal dynamics of the Hele-Shaw model. In particular it yields the continuity of normal velocities Eq. (2.3) and reintroduces surface tension, which is contained in the dynamical equation for the stream function through γ(θ).
As for Eq. (2.7), its right hand side is intended to reproduce Eq. (2.6), and therefore also Eqs. (2.1) and (2.2). If the phase-field θ has a kink shape, 1 − θ 2 is a peaked function which, when divided by ǫ, gives rise to the delta distribution for the vorticity. However, this only accounts for the γ in the weight of the delta. The part proportional to the viscosity contrast c must be put apart as the c ∇ · (θ ∇ψ) term because of the non-local character of ψ r (0 + ) + ψ r (0 − ). Finally, the time derivative is multiplied by ǫ to recover the laplacian (and not diffusive) behavior of the Hele-Shaw flow in the sharp-interface limit.
In spite of important differences, the proposed phase field equations Eqs. (2.7,2.8) contain certain similarities to the problem of a thermal plume in a Hele-Shaw cell under gravity [22]. In such a problem there is only one fluid heated from the centre of the channel. The heat diffuses towards the lateral walls, but the temperature profile is not linear, since the fluid density and viscosity decrease with temperature, so that the fluid in the middle of the channel raises because of buoyancy. As a result a so-called plume of hot fluid with a shape similar to the Saffman-Taylor finger, with a stationary upwards velocity and a width close to 1/2 is formed. Outside the plume the fluid is colder, and the transition between the two zones is relatively abrupt, so that one can think in terms of an interface of a certain small thickness. Thus, the equation for the phase field Eq. (2.8) could be thought as a diffusion equation for the temperature in a thermal plume. However, the available equations for that problem hold only for the steady state [22], whereas our phase-field model is intended to describe the whole dynamics. Generalization of the thermal plume equations to include the dynamics is not trivial for non-vanishing viscosity contrast. As a matter of fact, Ref. [22] must restrict itself to low viscosity contrasts -as it is the case in thermal plumes-, whereas we formulate the model for arbitrary viscosity contrast. An interesting difference is the term ǫ 2 κ(θ)| ∇θ| cancelling out Allen-Cahn law. The absence of that term in the thermal plume equations does not prevent the Hele-Shaw steady state equations to be recovered in the sharp-interface limit because of the lower power of ǫ used in the u · ∇θ term, but then Allen-Cahn law arises in the corrections at first order in the interface thickness. In contrast, by means of this ǫ 2 κ(θ)| ∇θ| term we achieve cancellation of the Allen-Cahn law even in such corrections, as we will see in section IV. Finally, another major difference in the case of thermal plumes is the absence of surface tension.

III. SHARP-INTERFACE LIMIT
In order to analyze the small-ǫ behavior of the phase-field equations, Eqs. (2.7,2.8), we expand their fields in powers of ǫ. The expected abrupt variations of these fields through the interface will make it necessary to perform two different expansions. In the interface region (inner region) we rescale the differential operators appearing in these phase-field equations by rewritting them in terms of the streched normal coordinate ρ ≡ r/ǫ (see Appendix). The expansions in the inner region will be matched order by order in powers of ǫ to those in the outer region (in the bulk far from the interface) where the coordinates are not rescaled. The outer and inner expansion are written respectively as a(r, s, t) = a 0 (r, s, t) + ǫa 1 (r, s, t) + ǫ 2 a 2 (r, s, t) + ... (3.1) A(ρ, s, t) = A 0 (ρ, s, t) + ǫA 1 (ρ, s, t) + ǫ 2 A 2 (ρ, s, t) + ... (3.2) where capital letters denote fields written in terms of the rescaled coordinate. This results in the following matching conditions: A 2 (ρ, s, t) = a 2 (0 ± , s, t) + ρa 1,r (0 ± , s, t) + ρ 2 2 a 0,rr (0 ± , s, t) ...
In practice, one does not find explicit solutions for the fields, but some set of equations for them. A sharp-interface model for the small-ǫ dynamics of the phase-field equations, Eqs. (2.7,2.8), is then given by the set of equations obeyed by the outer fields: Those obtained at lowest order in the interface thickness ǫ (O(ǫ 0 )) constitute the ǫ → 0 limit of the phase-field model, which we carry out in this section; whereas those obtained up to O(ǫ) represent what we will call (following Karma and Rappel [13]) a 'thin-interface' model, a model keeping finite interface thickness effects, such as the one derived in Sec. IV.
and iterating we get Due to Eqs. (3.5) and (3.7), θ = ±1 to all orders, and, therefore, the (1 − θ 2 ) term in Eq. (2.7) does not enter this outer limit, whereas the viscosity contrast term in that equation becomes ±c∇ 2 ψ, depending on the phase. Hence, Eq. (2.7) reads in this outer region which implies except for c = 1. Note that we have recovered the sharp-interface Eq. (2.1) in the ǫ → 0 limit. For c = 1, Eq. (2.1) is still recovered in the +1 phase ( viscous fluid), whereas in the -1 phase (inviscid fluid) the stream function turns out to be constant in time to all orders. Although the inviscid fluid does not enter the problem in this limit (see Eq. 2.2), it still has a non-trivial dynamics, since the stream function in it must evolve to keep satisfying Eq. (2.3), and therefore, strictly speaking, we do not really get the right sharp-interface limit for c exactly equal to one. However, the model can be applied to physical high viscosity contrast pairs of fluids. We shall come back to this point in section IV.
Hence we find the Θ 0,s term to vanish, and Eq.

IV. FIRST ORDER CORRECTIONS TO THE SHARP-INTERFACE LIMIT
In the phase field model the interface width and the convergence to the sharp interface limit is controlled by the small but finite value of the parameter ǫ. Then, the question of which value of ǫ is needed to accurately reproduce the actual Hele-Shaw dynamics for given values of the physical parameters B and c arises. This question can be qualitatively answered by noting the distinct roles played by ǫ in the phase-field equations, Eqs. (2.7,2.8): The ǫ factors appearing in ǫ 2 ∇ 2 θ, ǫ 2 κ(θ)| ∇θ| and 1 ǫ 1 2 √ 2 γ(θ)(1 − θ 2 ) all stand for the interface thickness, and this is required to be small compared to the longitudinal length scale |k| −1 of the interface: ǫ|k| << 1.
In contrast, the ǫ in ǫ ∂ψ ∂t has nothing to do with the interface thickness (and we will therefore denote it byǫ from now on), but its aim is to ensure that the stream function is laplacian and not diffusive in theǫ → 0 limit, which commutes with the ǫ → 0 one (the reader can convince himself of this by going through the limit again but now considering ǫ of O(ǫ 0 )):ǫ sets the time scale of the diffusion of the stream function through a given characteristic length of wavenumber k,ǫ (1±c)k 2 (see Eq. 3.8), which must be much smaller than the characteristic growth rate of the interface |ω| −1 , so that the stream function is slaved to the interface:ǫ |ω| k 2 << 1 ± c. We also realize that the viscosity contrast c can be arbitrarily raised, as long asǫ is correspondingly lowered. So our model is valid even for c → 1, as long as this limit is taken formally after the ǫ → 0 one.
The ǫ 2 in ǫ 2 ∂θ ∂t represents the relaxation time of the phase field towards the steady kink solution (see Eq. 2.8), which must be kept well below the interface growth time |ω| −1 for the phase-field to remain close to the kink profile during the interface evolution: ǫ 2 |ω| << 1. This factor must be the same that the one in ǫ 2ẑ · ( ∇ψ × ∇θ) in order to get the macroscopic equation Eq. (2.3). In fact there are at least two distinct powers of ǫ for this relaxation time (ǫ and ǫ 2 ) for which the right sharp-interface limit is achieved, and the corrections which we will compute would also be the same.
To sum up, there are at least two independent small parameters (ǫ andǫ) controlling the limit. When trying to approach macroscopic solutions by means of numerical integration of the phase-field equations, it is very convenient to vary them independently in order to save computing time, since both affect it [23].
A more quantitative answer to the question of the necessary values of ǫ,ǫ to get a given precision can be given by extending the asymptotic analysis of the previous section to first order in the interface thickness ǫ consideringǫ of O(ǫ). Thus, we will obtain a thin-interface model containing the corrections to the limit up to that order in ǫ andǫ.
On the one hand, we found that Ψ 0,s = −v n . We put this into Eq. (3.13) to get the differential equation for Θ 1 : with boundary conditions coming from the matching Eq.
According to the matching Eq. (3.3), the ρ → ±∞ asymptotics of Ψ 1 (ρ) should consist of a finite term, ψ 1 (0 ± ), and a diverging one, ρψ 0,r (0 ± ). For vanishing viscosity contrast the last integral in Eq. (4.2) is ρa 1 (s) and clearly does not contribute to the finite term ψ 1 (0 ± ). Then, since Θ 0 is an odd function of ρ, its integral with respect to ρ will be even, and ψ 1 (0 + )=ψ 1 (0 − ), i.e., the fluid velocity normal to the interface will be continous on it. For non-zero values of c, however, one must compute the integrals in Eq. (4.2), find their ρ → ±∞ asymptotics, and identify ψ 1 (0 ± ) and ψ 0,r (0 ± ). Requiring this latter quantity to be consistent with Eq. (2.2) for ψ 0 , one fixes a 1 (s), and putting this back into the identified ψ 1 (0 ± ) value, one finds where a 2 (s) is another arbitrary function of s. This will give rise to a discontinuity in the fluid velocity: In order to fix ∂ s a 2 (s), we compute the next order (O(ǫ 2 )) of Eq. (2.8) to get (see Eqs. (3.2,A10,A11,A19)): This has the same structure than Eq. (3.13) and an analogue solvability condition applies: Substitution of the expression for Ψ 1 obtained by performing the integrals in Eq. (4.2) into this condition and subsequent computation of the resulting integral fixes ∂ s a 2 (s) so that Finally, to get Ψ 2,ρ (+∞) − Ψ 2,ρ (−∞) we need Eq. (2.7) at (O(ǫ 0 )) (see Eqs. (3.2,A9,A10,A12,A17)): Integrating this from ρ → −∞ to ρ → +∞ we obtain where we have omitted integrals of odd functions of ρ. We use Eq. (3.20) to rewrite the integrand of the remaining integral as − γ 2 Θ 0 + a 1 (s). Θ 0 is an odd function of ρ and does not contribute to the integral, whereas a 1 (s) gives rise to a divergent term of the type [ρ] +∞ −∞ . According to the matching Eq. (3.4), ψ 1,r (0 ± ) corresponds to the finite part of Ψ 2,ρ (±∞), so that we find (4.10) which will leave the jump of the normal derivative of the stream function across the interface unaffected at first order in the kink width. Putting Eqs. (3.8,2.3 and 2.2 for ψ 0 ,4.7,4.10) together, we get an effective sharp-interface model for the dynamics of the θ = 0 level-set up to first order in ǫ andǫ: where Γ ≡ γ + c[ψ r (0 + ) + ψ r (0 − )] is the weight of the vorticity defined in Eq. (2.6) evaluated up to O(ǫ) and g ± (c) = 1 − 1 2 ) ln 1+c 1−c . Note that the desired corrections to the limiting equations Eqs. (2.1-2.3) in Eqs. (4.11) and (4.13) go asǫ and ǫ respectively, and the fact that Eq. (2.2) remains unaffected. Note as well that the correction in ǫ appearing in Eq. (4.13) has nothing to do with an Allen-Cahn law. So the ǫ 2 κ(θ)| ∇θ| term has cancelled this out even in the first order corrections.

V. LINEAR DISPERSION RELATION UP TO FIRST ORDER IN THE INTERFACE THICKNESS
In order to see how such corrections affect some relevant specific situation we compute the linear dispersion relation of a perturbation to the planar interface y(x) = Ae ωt+ikx for Eqs. (4.11-4.13). We make the ansatz ψ(x, y) = a ± Ae ωt+ikx−q ± |y| (5.1) inspired by the actual Hele-Shaw result, where now the coefficient a ± allows for distinct amplitudes in each phase for the stream function to satisfy the discontinuity in the normal velocities of Eq. (4.13), whereas the decay length q ± in the y-direction is set not only by the wavelength of the perturbation 2π/k, but also by the diffusion length in Eq. (4.11), which is also different in each phase. Thus, Eq. (4.11) yields In turn, taking into account that v n = ωAe ωt+ikx and γ = 2iAsign(k)ω 0 e ωt+ikx -where ω 0 = |k|(1 − Bk 2 ) is the actual Hele-Shaw growth rate-, Eq. (4.13) fixes a± to be Finally, Eq. (4.12) requires that the following dispersion relation is satisfied: This consists of the well known Hele-Shaw growth rate multiplied by a factor smaller than 1 carrying the corrections in ǫ,ǫ. We identify the conditions on ǫ,ǫ heuristically derived at the beginning of section IV to control how close this factor is to 1 and in general how close the stream function is to the actual Hele-Shaw one:ǫω/k 2 << 1 ± c (within p ± ) and ǫ|k| << 1 in Eqs. Since these corrections have a stabilizing effect, they could affect the selection of the steady finger width. As a matter of fact, Ben Amar already showed that the ∂ s (ŷ ·ŝ) term of Γ s in Eq. (4.13) on its own was capable of selecting a finger width greater than 1/2 [22]. Then, for small enough values of the physical surface tension (i.e., the physical selection mechanism), for which a width very close to 1/2 should be expected, this term could turn out to control the selection itself, so that an unexpected greater width could be obtained. Of course, this will not be the case if a sufficiently small value of the interface thickness ǫ is used, so that the condition ǫ|k| << 1 is satisfied for the length scale set by the surface tension: ǫ << √ B.

VI. CONCLUSIONS
We have introduced a phase-field model for Hele-Shaw flows with arbitrary viscosity contrast and shown it to yield the proper sharp-interface limit. We have actually found two independent small parameters (ǫ andǫ) and three distinct conditions on them to control the convergence to the sharp-interface limit ǫ,ǫ → 0. In particular,ǫ must be lowered when c is increased. A thin-interface model, i.e. an effective sharp interface model keeping finite-ǫ and -ǫ effects, has been derived for the dynamics of the phase-field model up to first order in both of these parameters. This thin-interface model has then been used to explicitely compute the finite-ǫ and -ǫ corrections to the Hele-Shaw result for a specific situation such as the linear regime, thus suggesting that the single-finger width selection could also be affected by these finite-thickness effects.
In the following paper [23] we perform numerical simulations of the phase field model Eqs. (2.7,2.8), and we explicitly vary the two small parameters ǫ andǫ independently. In this way we both control the simulation accuracy through the conditions mentioned to show how to reproduce the Hele-Shaw dynamics within this method, and explicitly check convergence in the interface thickness. on the interface r and s. To do this, one must precisely define the curvilinear coordinate system and compute its so-called scale factors: Consider the θ = 0 level-set and its intrinsic coordinates s (arclength along it) and r (signed distance to it, positive for a point with θ > 0), so thatŝ ×r =x ×ŷ. Let α be the angle going fromx toŝ. Then κ = α s is the θ = 0 level-set curvature. We introduce X, Y as the values of x, y for a point on the θ = 0 curve with a given value of s. By moving this point infinitessimaly along s we find that these values have changed in dY = ds sin α, dX = ds cos α. Consider also the coordinates of a point x, y with θ = 0 in terms of the values X, Y of its closest neighbour on the θ = 0 level-set and the signed distance between them. Taking into account that α is also the angle going fromŷ tor, one finds x = X − r sin α, y = Y + r cos α. Now one can compute the (positive defined) scale factors h 2 r ≡ x 2 r + y 2 r = 1 =⇒ h r = 1 (A1) h 2 s ≡ x 2 s + y 2 s = (X s − rα s cos α) 2 + (Y s − rα s sin α) 2 = = (cos α − rκ cos α) 2 + (sin α − rκ sin α) 2 = (1 − rκ) 2 =⇒ h s = |1 − rκ| = 1 − rκ (A2) Note that the last equality in Eq. (A2) requires that rκ < 1. In the inner region, where we make use of such formulae, this will hold as long as the interface thickness ǫ is much smaller than the curvature radius at any point of the interface, i.e. not too far from the sharp-interface limit. Otherwise the present analysis would break down, because one could always find a point such that rκ = 1, where h s would vanish, reflecting the fact that the change of coordinates has become ambigous in s.