Stéphane Popinet
email: stephane.popinet@upmc.fr

Numerical Models of Surface Tension

Keywords: well-balanced, Eulerian, levelset, volume-of-fluid, curvature, height-function, implicit, stability

Numerical models of surface tension play an increasingly important role in our capacity to understand and predict a wide range of multiphase flow problems. The accuracy and robustness of these models have improved markedly in the past twenty years, so that they are now applicable to complex, three-dimensional configurations of great theoretical and practical interest. In this review I attempt to summarise the most significant recent developments in Eulerian surface tension models, with an emphasis on well-balancing, curvature estimation, stability and implicit timestepping as well as test cases and applications. The advantages and limitations of various models are discussed, with a focus on common features rather than differences. Several avenues for further progress are suggested.

INTRODUCTION

Natural philosophers have noticed the effects of surface tension since the time of Aristotle, who recorded that flat pieces of iron or lead could float on the surface of water (Aristotle 350BC). A solid understanding of the phenomenon had to wait until 1805 when Young and Laplace independently published a comprehensive theory of capillarity inspired from the speculations of their predecessors (including Newton, von Segner and Monge: for a history of surface tension see the article by [START_REF] Maxwell | Scientific Papers of James Clerk Maxwell[END_REF] and its commentary by [START_REF] Pomeau | 5th Warsaw school of statistical physics, chap. Surface tension: from fundamental principles to applications in liquids and in solids[END_REF]). Besides its connection to modern understanding of the molecular nature of matter, surface tension is closely associated with major developments in differential geometry (by Clairaut, Monge and especially [START_REF] Gauss | Principia generalia theoriae figurae fluidorum in statu aequilibrii[END_REF] who formalised the analysis of minimal surfaces).

The coupling between volumetric fields and the geometry of surfaces is indeed at the heart of theoretical and numerical models of surface tension. From an analytical perspective, taking into account the (nonlinear) boundary conditions imposed on a surface which is itself part of the unknown solution presents formidable difficulties, so that analytical solutions are usually restricted to small surface deformations and/or static configurations (many of which were already obtained by Young, Laplace, Gauss and later Plateau).

A natural way to obtain numerical approximations for surface tension is to use boundaryconforming discretisations, where a Lagrangian description of volumetric fields is constructed so that the boundaries of volume elements coincide with the moving interface or free surface. Imposing the boundary or jump conditions given by surface tension is then relatively straightforward (see e.g. [START_REF] Fyfe | Surface tension and viscosity with Lagrangian hydrodynamics on a triangular mesh[END_REF] for an early example). The price to pay for this conceptual simplicity is the geometric complexity of Lagrangian methods, where the mesh needs to adapt to the underlying deformation of space.

Due to the necessity to deal with very large deformations, Eulerian field descriptions are a natural choice for fluid mechanics (in contrast with solid mechanics). This can be coupled with either a Lagrangian or an Eulerian representation of the interface. Lagrangian interface representation methods include the pioneering Marker-And-Cell method of [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], the immersed-boundary method of [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], the front-tracking method of Tryggvason and collaborators [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF][START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF] or the marker technique of [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]. Eulerian interface representations include essentially the Volume-Of-Fluid (Scardovelli & Zaleski 1999), levelset [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF] and phase-field methods [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF].

Surface (interface) and volume (fields) representations are then coupled through: a) interface kinematics i.e. the transport by the Eulerian velocity field of the Lagrangian or Eulerian interface description. b) interfacial dynamics i.e. changes in material properties (density and viscosity) and boundary or jump conditions associated with the interface. For Lagrangian interface representations transport is simple and accurate, however difficulties arise for large deformations and especially breakup and coalescence. Eulerian interface descriptions on the other hand can deal transparently with changes of topology but lead to more complex transport schemes. The solution of this kinematic problem has seen major progress in the past twenty-five years, with the development of higher-order geometric VOF methods which guarantee non-diffusive, sharp interface motion [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF], conservative levelset [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF], Feng et al. 2011) or coupled VOF-levelset methods [START_REF] Sussman | A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF].

The situation for interfacial dynamics, and surface tension in particular, is more complex and a wide range of methods or their combinations are available, often with limited information on their relative merits. The aim of this review is to highlight the most sig-nificant developments in Eulerian surface tension models in the past twenty years, with a specific emphasis on well-balancing, curvature estimation techniques, stability and implicit timestepping as well as test cases and applications.

SURFACE TENSION FORMULATIONS

Let us consider the Navier-Stokes equations for incompressible flow with surface tension

∂tρ + u • ∇ρ = 0 ρ (∂tu + u • ∇u) = ∇ • µ ∇u + ∇ T u -∇p + f σ ∇ • u = 0
with ρ and µ the variable density and viscosity, u the velocity, p the pressure and f σ the surface tension force per unit volume. The equation of continuity requires the solution of the kinematic problem mentioned in the introduction i.e. (non-diffusive) interfacial transport.

The assumption of incompressibility means that pressure loses its thermodynamic definition and reduces to the constraint necessary to impose a divergence-free velocity. This can be used to rearrange body forces and will be used later in this review.

Integral Formulation

Perhaps the most natural way to derive an expression for f σ is to consider the forces acting on a two-dimensional curve under tension, as was done by Young in 1805. In the case of interfaces between fluids (this is different for thin membranes), two-dimensional tension is simply a force per unit length tangential to the curve, which can be expressed as σt, with t the unit tangent vector and σ the surface tension coefficient. If we now consider an elementary volume Ω, intersected by the curve in two points A and B, the total tension force acting on Ω is

Ω f σ = B A σdt = σBtB -σAtA (1)
where the second integral is along the interface, and where σA and σB are the (possibly different) surface tension coefficients at A and B. The resultant of surface tension forces on the control volume thus reduces to the sum of tensions at the entry and exit points of the interface.

From a numerical perspective this formulation has several advantages: a) It only involves low-order derivatives of the geometry which should lead to accurate numerical estimates ; b) In a manner similar to flux-based integration of the divergence operator in finite-volume schemes, the contributions of surface tension forces on neighboring control volumes cancel out exactly (the directions of the unit tangent vectors in (1) are simply reversed). This can be interpreted as the divergence of a surface stress and is related to the Continuum-Surface-Stress formulation of [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF]. This ensures exact local and global momentum conservation for surface tension. In particular the net force exerted on the fluid by a closed interface is exactly zero. This is not the case for the volumetric force formulations discussed later in this review.

This formulation has been used in combination with a Lagrangian representation of the interface by [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] who obtained very accurate results (compared to methods available at the time), in particular for Laplace balance and capillary waves.

This required a consistent finite-volume discretisation of the pressure gradient, taking into account the pressure jump across the interface at the intersection points.

In the context of front-tracking techniques, the generalisation of (1) to three dimensions [START_REF] Weatherburn | Differential geometry of three dimensions[END_REF]) has been used by [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF] to define a globallyconservative interfacial force. This has not been combined however with the corresponding finite-volume integration on the Eulerian grid (rather the surface force is spread over the Eulerian grid using a smoothing kernel), so that the final scheme shares some of the properties of the volumetric force formulations of the next section.

Volumetric Formulation

Using the first Frenet formula for parametric curves dt = κnds, with κ the curvature, n the unit normal and s the curvilinear coordinate, gives (for a constant surface tension coefficient)

Ω f σ = B A σdt = B A σκnds = Ω σκnδs
with δs a surface Dirac δ-function which is non-zero only on the interface (see [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF] appendix B for a detailed derivation). Note that in the case of variable surface tension, tangential (Marangoni) stresses need to be added to this formulation (while they are naturally taken into account in the integral formulation).

Volumetric formulations are based on a numerical approximation of the surface Dirac function which allows a direct evaluation of the volumetric force σκnδs. This approach can be traced back to the original "immersed boundary" method of [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], although this link is not always acknowledged in subsequent articles. Most methods belong to this category independently of the way the interface is represented. The idea is to use the relation

σκδsn = σκ∇H(x -xs)
where H is the Heaviside function and xs is the position of the interface. The next step is to choose a suitable numerical approximation H of the Heaviside function which will typically depend on a (small) parameter such that lim →0

H = H

The small parameter is a length scale related to the characteristic thickness of the interface. Depending on this choice, which is linked to the type of interface representation, different methods can be constructed.

For the Continuum Surface Force (CSF) method of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], the interface is represented through the volume fraction field c and one simply sets H = c with = ∆, the mesh size. In the original CSF paper c is replaced with a filtered or smoothed version c which increases the characteristic interface thickness . This is not necessary at this stage however, as it is mostly related to a different issue: the estimation of interfacial curvature, which will be addressed later.

The original method of [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] is designed for elastic membranes, rather than interfaces, and relies on an explicit Lagrangian description of the membrane. [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] extended this approach to interfaces, using an indicator function reconstructed from the Lagrangian description, as an approximation of the Heaviside function. [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] proposed to use a representation of the interface as the zero contour of a levelset function φ. The Heaviside function is then approximated as a smooth function such as

H (x -xs) = H (φ(x)) =      0 if φ(x) < - 1 if φ(x) > 1+φ/ +sin(πφ/ )/π 2 otherwise (2)
In the common case where φ is chosen to be the signed distance to the interface, is the characteristic interface thickness. See [START_REF] Engquist | Discretization of Dirac delta functions in level set methods[END_REF] for a detailed discussion. Another type of surface tension discretisation usually associated with levelset interface representations is the Ghost Fluid Method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF]. Although this may not be obvious at first, this method can also be recast as a volumetric force formulation. Rather than including the surface tension as a surface force, Fedkiw et al. proposed to directly include the corresponding jump condition in the discretisation of the pressure gradient operator. This is essentially the "immersed interface" approach of [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF], itself a formalisation of Peskin's "immersed boundary" method.

Let us consider a simple second-order discretisation of a one-dimensional pressure gradient operator, which can be written

(∇p) i-1/2 = pi -pi-1 ∆
with ∆ the grid spacing. The Ghost Fluid Method (GFM) proposes to modify this operator locally in order to take into account the pressure jump induced by surface tension. This leads to the following scheme

(∇ p) i-1/2 = 1 ∆ p - i -pi-1 if the interface is in[xi-1 : xi] pi -pi-1 otherwise with the ghost fluid value p - i = pi ± (σκ) i-1/2
, where the sign of the jump depends on the orientation of the interface. This can be rewritten as

(∇ p) i-1/2 = (∇p) i-1/2 -(σκδ) i-1/2 with δ i-1/2 = ±1/∆ if the interface is in [xi-1 : xi] 0 otherwise (3)
Note that δ i-1/2 is indeed a consistent approximation of a Dirac delta function. The corresponding approximation of the Heaviside function is

Hi = 1 if xi is "inside the interface" 0 otherwise (4) and δ i-1/2 in (3) is δ i-1/2 = (∇H) i-1/2 = Hi -Hi-1 ∆
The GFM is indeed naturally suited to a levelset representation of the interface since the "inside the interface" condition in (4) is then φi > 0. Note however that it is also applicable to a VOF representation of the interface. The standard CSF approximation of the Dirac delta

δ i-1/2 = (∇c) i-1/2
can simply be replaced with δ i-1/2 = (∇H) i-1/2 where H is the equivalent of the GFM approximation (4)

Hi = 1 if ci > 0.5 0 otherwise
To sum up, the CSF method (typically combined with a VOF interface representation), the smoothed Heaviside method (typically combined with levelset or front-tracking) and the GFM method (typically combined with levelset) can all be summarised as σκδsn = σκ∇H with the different approximations for H summarised in Table 1.

Method Heaviside function H CSF (VOF) c

Smooth (Levelset/Front-Tracking) The corresponding graphs for H and δ (normalised by ) are represented in Figure 1. The fact that the approximation of the Heaviside function for the Ghost-Fluid method looks exact has led to the claim that the GFM leads to a "sharp" interface representation, in contrast with the other methods. It is clear however that the approximation of the Dirac function is not significantly sharper for the GFM than for the other methods. In fact all methods lead to a characteristic interface thickness of order = ∆. This is obvious for the Dirac approximations but is also true for the "sharp" Heaviside approximation of the GFM method, because this approximation is insensitive to a shift of ± /2 of the actual position of the interface.

     0 if φ < - 1 if φ > [1 + φ/ + sin(πφ/ )/π]/2 otherwise GFM (Levelset/VOF) 1 if φ > 0(levelset) or c > 0.5(VOF) 0 otherwise

WELL-BALANCED METHODS

A classic and non-trivial problem when designing numerical schemes for partial differential equations is ensuring that specific equilibrium solutions of the (continuous) equations are recovered by the (discrete) numerical scheme. This problem arises for a wide range of equations and applications: hyperbolic systems of conservation laws with source terms (LeVeque 1998), hydrostatic balance in σ-coordinate ocean or atmosphere models [START_REF] Mesinger | On the convergence and error problems of the calculation of the pressure gradient force in sigma coordinate models[END_REF], topographic source term in shallow-water equations [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]) and more. Numerical schemes which recover these equilibrium solutions are often called well-balanced methods.

In the case of surface tension, setting u = 0 in the Navier-Stokes equations, we get the equilibrium condition -∇p + σκnδs = 0 which is verified for

[p] = σκ, κ = constant,
where [p] is the jump of the pressure across the (spherical) interface. This is the well-known Laplace's relation between pressure and surface tension in a drop in equilibrium. If we assume that each of the (non-zero) terms in the Navier-Stokes equations are approximated with second-order accuracy, the leading-order discrete equilibrium condition can be written

∇p + Cp∆ 2 = σκnδs + Cσ∆ 2
where Cp and Cσ are scheme-and solution-dependent factors which control the pressure gradient and surface tension errors respectively. If special care is not taken, these errors have no reason to cancel and the best that one can hope is that the equilibrium condition is verified to within the accuracy of the scheme i.e. asymptotically to within O(∆ 2 ). While this may be acceptable for some applications, for example when surface tension is not the dominant force and/or when the interface is far from this equilibrium solution, in practice two-phase interfacial flows often include drops or bubbles close to Laplace's equilibrium. For many years numerical methods were not able to recover this simple equilibrium solution and the resulting quasi-stationary velocity pattern became known as spurious or parasitic currents, of varying intensity depending on the method. Interestingly the volumetric force formulation leads to a very simple condition for discrete equilibrium, as first pointed out by [START_REF] Renardy | PROST: a parabolic reconstruction of surface tension for the volumeof-fluid method[END_REF]. Indeed we have the (exact) relation -∇p + σκnδs = -∇p + σκ∇H = 0 which, in the case of constant σ and κ, can be approximated as -∇ (p -σκH) = 0 where ∇ is a numerical approximation of the gradient (assumed to be a discrete linear operator). The exact discrete numerical solution, which guarantees exact balance between surface tension and pressure in the case of constant κ, is then simply

p = σκH + constant, ( 5
)
where H is one of the approximations in Table 1.

For a staggered, MAC discretisation where velocity components are stored at mid-points, and pressure at the centers of a regular Cartesian grid, a simple well-balanced scheme is

- pi+1 -pi ∆ + (σκ) i+1/2 Hi+1 -Hi ∆ (6)
Unfortunately this was not the scheme adopted by Brackbill et al. in the original CSF article of 1992 (for which H = c). They instead chose to compute the mid-point surface tension force as the average of cell-centered values. This means that different operators are used to compute the gradient of the pressure (at mid-points) and volume fraction (at centers) and thus breaks the well-balanced property. The motivation for this choice was probably related to the desire to re-use the normal components also used to compute κ = ∇ • n at cell centres.

We thus see that an explanation for the relative confusion around the issue of spurious currents is linked to the contradictory requirements of well-balancing and accurate curvature estimation. Well-balancing requires that the gradient of a (usually) smooth function (the pressure) be estimated using the same discrete operator as that used to estimate the gradient of a discontinuous function (the Heaviside approximation/volume fraction). It is also wellknown that such a gradient estimate will lead to a poor approximation of nδs and an even poorer estimate of the curvature κ = ∇ • n. On the other hand, when looking for a better estimate of κ, one may be tempted to also use the better estimate of nδs associated with it, unfortunately this will break the well-balanced property. This is what happened for the original CSF formulation and many other early implementations.

It is thus necessary to decouple the estimate of nδs and the estimate of κ. Note that this is what happens naturally for levelset methods. In these methods the normal is estimated as n = ∇φ and the curvature as κ = ∇ • n, however nδs is estimated as ∇H where H is computed using either a smooth function or the GFM (see Table 1).

The reasoning on well-balancing applied to surface tension, can (and should) be applied to any system for which a source term can be balanced by the pressure gradient. If this source term can be expressed as the gradient of a potential, a consistent discretisation of the pressure and potential gradients will lead to a well-balanced scheme.

If we assume that a body force can be expressed as the gradient of a potential φ times a material property f , we have

-∇p + f ∇φ = -∇p + ∇(f φ) -φ∇f = -∇p -φ∇f with p = p -f φ.
If we now assume that the material property jumps from a constant value f0 to a constant value f1 across the interface, we have

f (H) = (f0 -f1)H + f1 = [f ]H + f1
where H is the Heaviside function associated with the interface. This gives

-∇p + f ∇φ = -∇p -[f ]φ∇H = -∇p -[f ]φnδs (7)
The body force can thus be replaced by an interfacial force equal to the potential times the jump of material properties across the interface. Surface tension itself can be reformulated using this general relation as

-∇p -σH∇κ = -∇p + σκnδs (8)
with p = p -σκH (using φ = κ and f = -σH in ( 7)). This has the advantage of giving a naturally-balanced scheme for a constant curvature, without any constraint on the choice of the discrete gradient operators (aside from ∇ constant = 0). See [START_REF] Ghidaglia | Capillary forces: A volume formulation[END_REF] for a detailed derivation.

Similarly the acceleration of gravity can be reformulated as 7)). The hydrostatic pressure is thus substracted analytically which is advantageous for ensuring exact hydrostatic balance across the interface and also simplifies boundary conditions for open domains (see e.g. [START_REF] Wroniszewski | Benchmarking of Navier-Stokes codes for free surface simulations by means of a solitary wave[END_REF] for an application to breaking waves).

-∇p + ρg = -∇p -[ρ]g • xnδs with p = p -ρg • x (using φ = g • x and f = ρ in (
More complex models can also be simplified using relation (7). For example [START_REF] Mahady | A numerical approach for the direct computation of flows including fluid-solid interaction: Modeling contact angle, film rupture, and dewetting[END_REF] propose to discretise the disjoining pressure induced by a potential of the form φ = (h /y) m where y is the distance to the substrate, as

-∇p + K∇φ = -∇p -φ[K]nδs
which avoids dealing with the divergence of the potential close to the substrate (this divergence is balanced analytically by the corresponding divergence of the pressure p).

CURVATURE ESTIMATION

The accurate estimation of interface curvature is of course central to the performance of volumetric surface tension force models (for integral formulations only the interface normal is required). Although the details depend on the particular technique used for interface representation, one can essentially distinguish two classes of methods: a) direct derivation from the implicit representation of the interface, b) discrete differential geometry operators applied to an explicit description of the location of the interface (provided directly by the interface tracking representation or reconstructed from an implicit description of the interface).

Levelset Methods

Method a) is the primary argument in favour of levelset interface representations. In the ideal case where the levelset function φ is the signed distance to the interface, the curvature is simply obtained using the relations

n = ∇φ and κ = ∇ • n (9)
The levelset function φ is smooth and discrete schemes can be easily constructed to estimate both n and κ to any order of accuracy. This is the approach adopted in early articles where levelset is used to represent interfaces and surface tension [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. Although useful schemes can be obtained with this approach, it has an important limitation which is evident when considering the well-balanced property of the previous section. Let us consider the simple well-balanced discretisation

- pi+1 -pi ∆ + (σκ) i+1/2 H(φi+1) -H(φi) ∆ , ( 10
)
where H is a suitable approximation of the Heaviside function (e.g. ( 2) or the GFM approximation in the last row of Table 1). A simple estimate of the curvature is

κ i+1/2 = [∇ • (∇φ)] i+1/2 (11)
where the divergence and gradient operators are suitable discrete approximations. To recover Laplace's equilibrium relation ( 5), the estimate of the curvature needs to be constant for a spherical droplet. Note however that the relations (9) as well as their discretisation (11) define a curvature everywhere in the domain, not just on the interface (φ = 0). The curvature given by ( 9) is indeed the curvature of the local contour, not that on the interface.

In the case of a spherical interface, (11) will thus return an estimate of the curvature which will not be constant and Laplace's equilibrium will not be guaranteed. The deviation from the equilibrium solution will be controlled by the distance between the interface and the points where the curvature is estimated (i.e. at i + 1/2) times the value at these locations of the gradient of the Heaviside function. This distance obviously scales like the thickness of the interface (i.e. essentially the grid size) and will thus tend toward zero with spatial resolution, which makes the method consistent but not well-balanced. This limitation of "naive" levelset methods has been recognised for example by [START_REF] Sussman | A stable and efficient method for treating surface tension in incompressible two-phase flow[END_REF] who advocate the use of a height-function curvature calculation (described later in this section), rather than relations (9), even in the case of a levelset description of the interface.

A closely related issue is the explanation of the improved performance of the GFM relative to the classical "smoothed Heaviside" levelset formulation. In the original article of [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] the decrease in spurious currents observed with the GFM is attributed to the "sharper" surface tension representation provided by the GFM, a conclusion which is contested in the latter article by [START_REF] Francois | A balancedforce algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF], in connection with their "Sharp-Surface-Force" (SSF) model (this model is actually identical to the Dirac formulation (3) of the GFM). The improvement is not due to the difference in Dirac approximations summarised in Table 1, but to the interpolation formula proposed by Kang et al. to estimate the pressure jump (i.e. the curvature) at the location of the interface (based on the earlier work by [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF]):

κ i+1/2 = κi|φi+1| + κi+1|φi| |φi| + |φi+1|
i.e. a distance-weighted average of the cell-centered values of the curvature, which are computed using (9). In the case of a spherical interface, it is clear that this formulation will lead to values of the curvature in (10) which are much closer to a constant than e.g. the simple average κ i+1/2 = (κi + κi+1)/2 often used for the classical levelset formulation. This has recently been confirmed by [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks[END_REF] who showed that the GFM/SSF levelset formulation leads to exact well-balancing, in contrast with the classical levelset formulation which gives significant spurious currents.

The benefits of ensuring a close-to-constant curvature estimate in the case of a spherical interface are also obvious for the alternative surface tension formulation (8), which has a built-in well-balanced property (irrespective this time of the specific choice of gradient operators).

The main drawback of levelset formulations is the lack of discrete volume/mass conservation. This can be minimised by employing frequent redistancing/re-initialization of the levelset function. Unfortunately this redistancing step tends to perturb the curvature and prevents the system from reaching the constant curvature required to guarantee exact balance (see [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks[END_REF].

Smoothed Volume Fraction

The normal and curvature relations (9) can in principle be applied to any field which provides an approximation of the interface position as a local isocontour. This was the observation which lead [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] to propose the original CSF model where the curvature is approximated using

n = ∇c |∇c| and κ = ∇ • n
where c is a "smoothed" (i.e. diffused) version of the sharp volume fraction field. The smoothing was adopted to try to circumvent the difficulty of differentiating a discontinuous function. While results do improve when stronger diffusion is applied, several studies have since showed that this method is not consistent for curvature i.e. errors on curvature tend to increase with spatial resolution (see e.g. Wiliams et al. (1998), [START_REF] Cummins | Estimating curvature from volume fractions[END_REF] or [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]). This inconsistency of the curvature estimate, often combined with the lack of well-balancing of the gradient terms described previously, explains most of the difficulties encountered with early CSF implementations. For many years this has been a strong argument in favour of levelset methods which, even for their naive not wellbalanced versions, are at least consistent and give much better results for surface tension. Although this argument is often repeated in recent publications, the smoothed volume fraction method should be considered obsolete since much better alternatives now exist, as described next.

Height-Functions

The height-function method relies on the simple observation that one can always define a local coordinate system in which a surface is described as the graph of a function. The simplest case is that of a nearly-horizontal interface which can simply be described by y = hy(x) and for which the curvature is given by

κ = h y 1 + h 2 y
Provided the values of hy (the "height-function") are known exactly at discrete locations, one can easily derive discrete schemes to estimate κ at any order of accuracy. This reasoning is easily generalised to any number of spatial dimensions. Of course an issue arises when the slope of the interface tends to infinity and the graph becomes multivalued. One can then simply switch to a different local representation of the interface which is well-behaved such as x = hx(y).

This simple idea has a long history, possibly starting with [START_REF] Poo | A computational method for determining curvatures[END_REF]. It was first described comprehensively in the context of two-phase flows by [START_REF] Sussman | A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles[END_REF] and later analysed in detail by [START_REF] Cummins | Estimating curvature from volume fractions[END_REF]. The key to the success of this method is having access to sufficiently accurate discrete values of the height function. Indeed, as double differentiation is required to estimate curvature, the height function needs to be known with higher than second-order accuracy in order to get a consistent (i.e. converging) estimate of the curvature. In the case of volume-of-fluid methods, if we assume that the discrete volume fractions ci,j are known exactly (for example using volume fractions obtained by analytical integration based on the exact interface [START_REF] Bnà | Numerical integration of implicit functions for the initialization of the VOF function[END_REF])), we can write without approximation

x i+1/2 x i-1/2 hi(x)dx = ∆hi = j=∞ j=-∞ ci,j + constant ( 12
)
where hi is the exact average height function value for a given column. The ±∞ limits of the sum on the right-hand-side simply indicate that vertical summation of the volume fractions is performed on the entire column. A symmetric relation is of course derived for integration along rows, which allows to perform the coordinate system rotation described earlier. Once these exact average values are computed one can use simple differencing such as

κi = ( hi-1 -2 hi + hi+1)/∆ 2 1 + [( hi+1 -hi-1)/(2∆)] 2
to obtain second-order accurate estimates of the curvature. Both second-and fourth-order convergence (using five-points stencils) have been demonstrated in practice for circular interfaces [START_REF] Cummins | Estimating curvature from volume fractions[END_REF][START_REF] Sussman | High-order techniques for calculating surface tension forces[END_REF][START_REF] Bornia | On the properties and limitations of the height function method in two-dimensional Cartesian geometry[END_REF]. Generalisation of the method to three-dimensions is straightforward.

Besides its conceptual simplicity and good convergence properties, the method has several advantages which are particularly relevant in the context of the well-balanced property discussed earlier: it naturally defines curvatures "on the interface", it is insensitive to the exact distribution of volume fractions within a column: this is important for robustness since volume-of-fluid schemes often lead to the formation of small interfacial fragments. Using these small fragments as markers of interface position can lead to erroneous curvature estimates (for example when using the interpolation techniques described in section 4.4). The spatial averaging performed by the H-F method naturally avoids these artefacts, more generally, the connection between interface representation (and transport) and surface tension is very direct when combining volume-of-fluid, balanced surface force discretisation such as (6) and height-function curvature estimation: this tight coupling minimises the possibility of uncontrolled numerical modes and gives a robust method.

Several issues must be addressed when implementing the method in practice however. The summation limits on j in expression (12) need to be specified. Early implementations only used a fixed number of cells (typically seven) in each column. If one considers a nearly horizontal interface, described using a non-diffusive VOF method, only three cells in the vertical direction are necessary to obtain a consistent value for the average height. On the other hand some interface configurations may require up to nine cells in the vertical to obtain consistent height functions (see e.g. [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] for a discussion and examples). More generally, while consistent height-function approximations can always be obtained when the spatial resolution is high enough (i.e. when the product κ∆ is small enough), things become more complicated as soon as this product is larger than around 1/5.

A first useful step is to use a variable stencil height which adjusts automatically to the local topology of the interface, in order in particular to ensure that vertical summation is only applied across a single interface. This also permits to use optimal stencils (e.g. 3 instead of 7 cells) in the simpler cases. See [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] and Lopez et al. (2009) for examples of this approach.

As spatial resolution decreases further, consistent height functions are increasingly difficult to obtain: when one reaches κ∆ 1, at best a single value of the height function can be constructed in each direction and differentiation to obtain curvature is no longer possible. Although one cannot reasonably expect accurate solutions at such low resolutions, it is important to ensure that the overall method remains robust across these transitions. This can be achieved by switching progressively to the interpolation methods described in the next section. A complete example of this approach is described in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. An alternative approach is proposed by [START_REF] Owkes | A mesh-decoupled height function method for computing interface curvature[END_REF] who use a rotated height-function stencil.

Differential Geometry of Discrete Surfaces

The second class of methods for curvature estimation relies on an explicit, often local, discretisation of the surface. It is the natural method to use for front-tracking where the interface is defined explicitly, for example using Lagrangian vertices connected by triangular facets.

The approximation of surfaces using discrete elements has many applications besides interfacial flows, for example 3D laser scanning, mesh compression and denoising, minimal surfaces and geodesics. Approximations of differential quantities such as Gaussian and mean curvatures, principal directions, etc. has recently seen a surge of interest in the general context of computational geometry. An interesting review of the connections between the continuous mathematical concepts and their discrete equivalents is given by [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF]. An important result from this article is the derivation of a discrete equivalent of the mean curvature normal operator which is directly related to surface area minimization. This discrete operator mimics several of the important properties of the continuous operator; in particular it minimises a discrete version of the surface energy. Note that the the techniques described in this article have not, as far as I am aware, been applied yet to front-tracking codes.

Besides this discrete formulation, two main types of techniques exist depending on whether a globally continuous approximation is required or whether a local approximation is sufficient. In global methods, local surface patches are constructed and continuity conditions are imposed between patches, which leads to the inversion of a system for the entire surface. Splines and B-splines belong to this family. In the context of interfacial flows, spline curves have been used for example by [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] to connect Lagrangian markers in two dimensions, while B-splines have been used by [START_REF] Torres | The point-set method: front-tracking without connectivity[END_REF] to obtain a global interface representation, also in two dimensions. These methods are accurate and benefit from the robustness provided by a global reconstruction (which can be seen as minimising a global surface energy functional i.e. the bending energy for 2D splines [START_REF] Birkhoff | Piecewise polynomial interpolation and approximation[END_REF]), however they are also expensive and complex, particularly in three dimensions.

In practice local surface approximations are much more common. A standard technique is to approximate the surface by least-square minimisation of the quadratic form

i wi(xi • Axi + n • xi + b) 2 , ( 13
)
with n a unit vector approximating the normal and A a symmetric matrix. Further constraints can be used to reduce the number of free parameters. For example, [START_REF] Renardy | PROST: a parabolic reconstruction of surface tension for the volumeof-fluid method[END_REF] argue that one can get a consistent, second-order approximation by imposing that the axis of the paraboloid defined by (13) be aligned with n i.e. An = 0. Popinet ( 2009) uses a similar argument to reformulate the approximation problem in a local coordinate system aligned with n. A sufficient number of points xi are then chosen in a local neighborhood to make the system (13) invertible. The coefficient wi are optional ad-hoc weights accounting for differences in accuracy of the positions. Besides the normal n, the mean curvature is then given by κ = 2 tr(A). With the possible exception of triangulated surfaces of known connectivity [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF], it is clear that such local approximations are not nearly as simple or systematic as levelsets or even height-functions. While the least-square minimisation is not particularly complex or computationally expensive, the logic of point (and weight) selection can be complicated and somewhat fragile. Besides its application to front-tracking, variants of this method have been used in particular in:

• The PROST method of [START_REF] Renardy | PROST: a parabolic reconstruction of surface tension for the volumeof-fluid method[END_REF]: in this article it is not the discrete interface locations which are fitted through ( 13), but the volume fractions given by the intersection of the quadratic surface with each cell of a local stencil. The resulting minimisation problem is non-linear and very expensive to solve numerically, but provides accurate normal and curvature estimates directly from the volume fraction field. Ad-hoc weighting is used to avoid computational modes and to increase the robustness of the method. • The generalised Height-Function method of [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]: as described earlier, when κ∆ < 1/5 the number of consistent discrete heights may not be enough to allow differentiation. It is then necessary to use other types of surface approximations. Ideally switching from one approximation (e.g. height-function) to the next (e.g. least-square minimisation) should be done while minimising the potential "jumps" in the estimated curvature. A first step is to obtain interface positions by combining horizontal and vertical "heights" and applying (13). If the number of interface positions obtained in this way is too small, one can switch as a last resort to interface positions given directly by the VOF geometric reconstruction. [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] and many subsequent applications of the method to complex problems have shown that this method is both accurate and robust. The computational cost is limited, since the relatively expensive least-square problems need to be solved only in marginal cases. The main cost is rather the added code complexity, which is quite manageable however (see for example [START_REF] Popinet | Basilisk[END_REF] for a complete open-source implementation).

Mixed Methods

The relative advantages and drawbacks of various methods has led to techniques seeking to combine several interface representations. The primary motivation was the improvement of interface kinematics, accuracy and mass conservation properties in particular, but accuracy of surface tension representation was also an important factor. One of the better known mixed method is the CLSVOF method of [START_REF] Sussman | A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF] , which couples a geometric VOF representation (for conservative transport) with a levelset (for simple curavture estimation). Many variants of this approach have since been developed and successfully applied to complex interfacial flow problems. The main weakness of this technique is common to all mixed methods: it is difficult to switch between interface representations without loss of accuracy on the interface position. These errors can become dominant when computing curvature. For example, [START_REF] Cummins | Estimating curvature from volume fractions[END_REF] have analysed the curvature errors induced by their Reconstructed Distance Function (RDF) procedure, which is close to that employed in CLSVOF methods. Using a second-order-accurate RDF reconstruction, they demonstrate that, as expected, normals estimated from the RDF converge at a first-order rate, while curvature errors do not converge with resolution. More generally, one can expect that this reconstruction error will affect most types of redistanciation algorithms used to ensure the consistency of levelset representations. This should be particularly clear when considering the well-balanced cases discussed earlier: in these cases redistanciation (or CLSVOF coupling) will lead to a continuous injection of curvature errors which will prevent reaching an equilibrium curvature. This assumption has been confirmed for the spurious current test cases discussed later in this review.

Another common class of mixed interface representations combines fronttracking/Lagrangian particles with levelset/isosurfaces/VOF. From a kinematic point-ofview the goal is to combine the high-accuracy of Lagrangian advection schemes with the topological flexibility of levelset/VOF interface representations. The improved kinematic properties of these schemes have been demonstrated in [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF] or [START_REF] Hieber | A Lagrangian particle level set method[END_REF] for particles + levelset and [START_REF] Aulisa | A mixed markers and volume-of-fluid method for the reconstruction and advection of interfaces in two-phase and free-boundary flows[END_REF] for particles + VOF. The Level Contour Reconstruction Method (LCRM) of [START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF], uses (temporary) interface markers reconstructed from the isosurface of the indicator (i.e. levelset) function. This eliminates the need to keep track of surface connectivity (i.e. topology) of standard front-tracking techniques. The curvature is computed using the local triangulated isosurface reconstruction and surface tension is then implemented using balanced CSF applied to the indicator function. Although the convergence of curvature errors is not discussed in [START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF] (one can expect the same difficulties as for the RDF method), good results are obtained for the spurious currents test case.

STABILITY AND IMPLICIT TIMESTEPPING

As first discussed briefly by [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], a time-explicit discretisation of the surface tension term should lead to a stability constraint of the form

∆t < (ρ1 + ρ2)∆ 3 4πσ ≡ ∆tσ, ( 14
)
with ∆t the timestep and ρ1 and ρ2 the densities on either side of the interface. The physical justification is that the timestep must be small enough to resolve the fastest capillary waves in the system which are obtained for the wavenumber k = π/∆ and verify the dispersion relation cσ = σk/(ρ1 + ρ2) . A more formal stability analysis can be found for example in [START_REF] Sussman | A stable and efficient method for treating surface tension in incompressible two-phase flow[END_REF]. Note that the assumptions which lead to the 4π constant in ( 14) are debatable, see [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF] or [START_REF] Denner | Numerical time-step restrictions as a result of capillary waves[END_REF] for a detailed discussion.

Explicit schemes for the transport of interfaces are subject to the standard CFL constraint ∆t < ∆/|u| ≡ ∆t adv . The ratio of these two stability constraints is thus

∆tσ ∆t adv = (ρ1 + ρ2)|u| 2 ∆ 4πσ = √ We∆
where We∆ is the cell Weber number which estimates the ratio of inertial to surface tension forces. In the absence of viscosity, the minimum characteristic scale (i.e. radius of curvature) of interfaces is due to the balance of inertial and surface tension forces. This implies that a well-resolved numerical simulation will necessarily verify We∆ 1, so that the minimum characteristic scale of interfaces (e.g. minimum droplet size) is (much) larger than the mesh size. This means that the capillary timestep restriction ∆tσ must always be more restrictive than the CFL constraint for interface advection. For example if we consider a one millimeter air bubble rising in water with a terminal velocity of 0.1 m/s, we get

∆tσ ∆t adv = √ We∆ = √ WeDN -1/2 ≈ 0.1N -1/2
where WeD is the bubble Weber number (based on diameter D) and N is the number of grid points per bubble diameter. For a moderate resolution of N = 10 the capillary timestep is thus 32 times smaller than the advection timestep.

One might wonder whether viscosity can relax the stability condition. This could be particularly relevant for small-scale flows, such as those occurring in microfluidics devices, and has been studied in this context by [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF]. They conclude and demonstrate numerically that, for capillary-driven Stokes flows, the relevant stability criterion is ∆t < max µ∆ σ , ρ∆ 3 σ where the first term corresponds to a CFL condition built with the Stokes velocity σ/µ. The combined criterion thus becomes advantageous whenever the Stokes velocity is smaller than the speed of the shortest capillary waves i.e. σ/µ < σ/ρ∆ which gives

∆ < lµ = µ 2 ρσ ,
where lµ is the visco-capillary length, which only depends on the fluid properties. As a reminder lµ ≈ 14 × 10 -9 m and σ/µ ≈ 72 m/s for an air/water interface and lµ ≈ 2.5 cm and σ/µ ≈ 4.5 cm/s for an air/glycerol interface. Note that microfluidics devices typically use mineral oils which have viscosities and surface tensions comparable to that of water and channel widths of the order of ten microns. A well-resolved simulation will have a spatial resolution of order ∆ ≈ 0.1 microns which is still much larger than the visco-capillary length. It is thus not obvious that the relaxed Stokes criterion is beneficial for these applications, unless one considers significantly more viscous fluids and/or smaller channels. In practice, explicit surface tension stability is thus a stringent constraint for microfluidics applications, see [START_REF] Ling | Droplet migration in a Hele-Shaw cell: Effect of the lubrication film on the droplet dynamics[END_REF] for a detailed discussion.

Semi-Implicit Schemes

Time-implicit discretisations are a standard way to relax timestep restrictions. Fast modes are damped rather than resolved and do not constrain the overall stability. One of the first article discussing time-implicit schemes for surface tension is by [START_REF] Bänsch | Finite element discretization of the Navier-Stokes equations with a free capillary surface[END_REF] who proposes a variational formulation coupled with a Lagrangian finite-element discretisation of the interface (see also the interesting review by [START_REF] Buscaglia | Variational formulations for surface tension, capillarity and wetting[END_REF]). This is inspired by earlier work on the computation of minimal surfaces through discretisation of the mean curvature flow equation [START_REF] Dziuk | An algorithm for evolutionary surfaces[END_REF]). The extension of these ideas to an Eulerian interface description was presented by [START_REF] Hysing | A new implicit surface tension implementation for interfacial flows[END_REF]. The scheme is derived starting from the standard differential geometry relation [START_REF] Weatherburn | Differential geometry of three dimensions[END_REF])

∆xs = κn
with xs = xδs the Lagrangian interface position and ∆ the Laplace-Beltrami operator (or surface Laplacian). The surface tension force can then be written

f σ = σκnδs = σ∆xsδs
The idea is then to use the interface position predicted at the next timestep by

xs,n+1 = xs,n + ∆tun+1,
to obtain a semi-implicit discretisation of the surface tension force as

f n+1 = σ (∆xs) n+1 δs,n = σ (∆xs) n δs,n+∆tσ (∆u) n+1 δs,n = σ(κnδs)n+∆tσ (∆u) n+1 δs,n
The scheme is not fully implicit in particular because the timelevel is not taken into account for δs (which is always defined at time n). The scheme is equivalent to the addition of a surface viscosity proportional to ∆tσ which will dampen fast capillary waves and lead to stabilisation. Hysing demonstrates gains of one order of magnitude in timestep compared to standard explicit schemes. The method is also applied in a Finite-Volume/VOF context by [START_REF] Raessi | A semi-implicit finite volume implementation of the CSF method for treating surface tension in interfacial flows[END_REF] who find similar stability properties.

More recently [START_REF] Sussman | A stable and efficient method for treating surface tension in incompressible two-phase flow[END_REF] proposed to estimate the curvature at n+1 using the mean curvature flow equation

∂txs = σκn = σ∆xs (15)
rather than the full coupled Navier-Stokes system. The stationary solutions of this equation minimise the surface energy, which is clearly a desirable property when considering the stability of integrations with large timesteps, for which they obtain the stability condition ∆t ∆(ρ1 + ρ2) 2π ( 16)

They then numerically demonstrate improved stability compared to the standard discretisation, with increases in timestep of the same order as for the method of Hysing et al. These two equations are problematic however, since neither are dimensionally consistent. The analysis can be fixed by noting that (15) is not an equation describing the evolution of an interface under the effect of surface tension, despite its connection with minimal surfaces. Indeed dimensional consistency implies that the coefficient σ in (15) has dimensions L 2 T -1 i.e. that of a diffusion coefficient. The mean curvature flow equation ( 15) is simply a surface diffusion equation which will filter high-frequency surface modes and therefore stabilise the solution. This is well known in the computer graphics community where ( 15) is used to "denoise" surface meshes [START_REF] Desbrun | Implicit fairing of irregular meshes using diffusion and curvature flow[END_REF].

The scheme can be reformulated consistently as the levelset evolution equation [START_REF] Chopp | Computing minimal surfaces via level set curvature flow[END_REF] ∂τ φ = κ

where the pseudo-time τ has the dimension of a length squared. A filtered curvature can then be defined as

κ = 1 Λ 2 Λ 2 0 κdτ = φ(Λ 2 ) -φ(0) Λ 2
where φ(Λ 2 ) is obtained by advancing ( 17) from zero to Λ 2 , with Λ a characteristic smoothing length. The amount of smoothing required is found by considering the linear stability of the resulting scheme, which gives the stability condition

Λ > ∆t σ ∆(ρ1 + ρ2) = ∆tc∆,
with c∆ the capillary wave speed for the (smallest) wavelength ∆. The correct interpretation of the scheme of Sussman & Ohta is thus stabilisation by diffusive surface smoothing over the characteristic travel distance of the fastest capillary waves: Λ = ∆tc∆. The schemes of Sussman & Ohta and Hysing et al. thus work in a similar manner: added surface damping filters high-frequency modes and thus stabilises the solution. An important difference between the two schemes however, is that the method of Sussman & Ohta directly filters the interface position/curvature while the scheme of Hysing filters the (surface) velocity field. In particular, the scheme of Hysing will not affect equilibrium shapes for which u = 0, but the scheme of Sussman & Ohta will. More generally, these schemes are not very scale-selective filters (because they are both based on low-order differential operators) i.e. they will also significantly dampen lower-frequency modes which do not restrict stability.

Pushing this approach further, one can devise, at least formally, near optimal filtering schemes. If we consider the simpler case of a one-dimensional interface defined through its graph η(x, t), the Fourier transform of the corresponding curvature is given by κ

(k, t) = -k 2 η(k, t)
with η(k, t) the Fourier transform of the interface position and where we have assumed a vanishing interface slope. An optimal filtered curvature can then be defined in Fourier space as κ(k, t) = min 1, ρ1 + ρ2 σk 3 ∆t 2 κ(k, t) Computing the inverse Fourier transform and using the resulting filtered curvature will then ensure stability of the explicit scheme. The filtering is optimal since only the necessary (mode-dependent) amount of damping is added. This scheme works well in practice, for example using FFTs for periodic graphs in one or two dimensions, however generalising it to more complex topologies seems difficult. For front-tracking interface representations spectral mesh processing could be a solution. See the course by Levy et al. (2010) for an interesting introduction.

TEST CASES

Test cases are important for the development and assessment of new numerical schemes. Publicly-accessible automated test suites (see Basilisk (2013) andGerris (2003) for examples), cross-referenced with journal articles are also an excellent way of ensuring reproductibility and independent peer-review of the numerical results. In this section I will try to point out a minimal set of test cases for surface tension models as well as their shortcomings and common pitfalls.

Laplace's Equilibrium and Spurious Currents

Laplace balance between surface tension and pressure gradient provides a trivial equilibrium solution which is nevertheless difficult to reproduce numerically, leading to the production of numerical artefacts, the so-called spurious or parasitic currents. This was first observed and discussed in detail in the context of lattice Boltzmann methods for two-phase flows [START_REF] Gustensen | Lattice-Boltzmann studies of multiphase flow through porous media[END_REF]) and many variants of this test case have since appeared in the litterature.

An important issue is that of the timescale required to reach equilibrium. There are two natural timescales in this simple system: the period of oscillation scaling like Tσ = ρD 3 /σ with D the droplet diameter and the viscous dissipation timescale Tµ = ρD 2 /µ. The ratio of these two timescales is Tσ/Tµ = µ/ √ ρσD, the Ohnesorge number. To reach the asymptotic regime corresponding to the equilibrium solution, one then needs to make sure that the simulations are run on a timescale (much) larger than either of these two timescales. Detailed parameters for such a setup are provided for example in [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF] or [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF].

Note that many studies have been published which do not verify this condition of asymptotic convergence. For example, [START_REF] Francois | A balancedforce algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] used a version of this test where a few or even a single timestep are performed before measuring the amplitude of "spurious currents". In the case of well-balanced schemes, the cause of spurious currents after a few timesteps is only the deviation from constant of the initial curvature computed by the scheme. This deviation is better characterised by, for example, a convergence test on the curvature estimate for a spherical interface (see e.g. [START_REF] Cummins | Estimating curvature from volume fractions[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]), without running the risk of confusing several properties of the scheme (i.e. well-balancing versus curvature estimation).

On the other hand, if the test respects the asymptotic conditions t > Tµ and t > Tσ, one expects a consistent, well-balanced method to converge toward an interfacial shape which will ensure exact equilibrium (i.e. u = 0 to within machine accuracy). This interfacial shape may not itself be exact (i.e. an exact circle/sphere) and the evolution of the velocity around the interface is expected to reflect the physical evolution (through damped capillary waves) from the initial "perturbed" condition toward the numerical equilibrium solution. One of course expects this numerical equilibrium interface shape to converge toward the exact equilibrium (circular) shape as spatial resolution is increased. See [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] for a demonstration of this convergence in the case of the VOF and HF-CSF surface tension method. Note that such a convergence is not trivial however, since it requires that the scheme guarantees evolution towards a constant curvature.

An important extension of this test was proposed in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] where a constant background velocity field ensures uniform translation of the droplet across the grid in a spatially periodic domain. Laplace's equilibrium solution is of course still valid in the frame of reference of the droplet. This test is more relevant to practical applications, par-ticularly when considering low-velocity/high-surface tension cases such as microfluidics or multiphase flow through porous media. This test was used extensively in the interesting comparative study of [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks[END_REF] who underlined the detrimental effect of interface (and curvature) perturbations induced by either interfacial transport (for VOF) or redistancing/interface reconstruction (for Levelset and coupled Levelset/VOF). See also in this article the demonstration of consistent well-balancing for levelset methods (without redistancing).

Capillary Oscillations

Capillary oscillations around equilibrium solutions are the next logical step. Analytical solutions can be obtained through classical linear stability analysis in the limit of vanishing amplitude and viscosity both for planar and circular/spherical interfaces [START_REF] Lamb | Hydrodynamics[END_REF]. [START_REF] Fyfe | Surface tension and viscosity with Lagrangian hydrodynamics on a triangular mesh[END_REF] considered the oscillation of a two-dimensional elliptical droplet in an inviscid fluid, for which the oscillation frequency is given by

ωn = (n 3 -n) σ (ρ1 + ρ2)a 3 (18)
where the droplet shape is given in polar coordinates by r = a + cos(nθ). This is an extension to two phases of a result by [START_REF] Rayleigh | On the capillary phenomena of jets[END_REF] who considered the stability of the cross-section of a jet. Although many variants of this test case exist, one of the most challenging is for large density ratios (1/1000) without viscosity. Details can be found in e.g. [START_REF] Torres | The point-set method: front-tracking without connectivity[END_REF][START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF][START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF]). The total energy (surface plus kinetic) should remain constant and any decay is the sign of numerical dissipation which should be minimised. Conversely, increasing total energy is a clear signature of surface tension imbalance. This setup is a stringent test of the accuracy of surface tension representation since physical or numerical viscosity cannot intervene to limit spurious currents. In addition to minimising dissipation, good numerical schemes can give second-order spatial convergence in the estimated oscillation frequency compared to (18) (see [START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF] for results with different schemes). Analytical solutions can also be obtained when viscosity is taken into account. The simplest analysis leads to exponential damping of oscillating modes, however, as studied in detail by [START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF], this gives significant deviations (several percents) compared to initial-value solutions taking into account the time-dependence of vorticity diffusion into the medium. Prosperetti derived closed-form solutions for the Laplace transform of shape evolution both for planar and spherical interfaces [START_REF] Prosperetti | Free oscillations of drops and bubbles: the initial-value problem[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF]. These solutions are the basis for a now classical test case, first proposed in [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF], which considers the oscillations of a linearly perturbed, planar interface. Although less stringent than the inviscid case, due to a simpler geometry which is less affected by imbalance and spurious currents, this test evaluates the quality of the full coupling between interfacial motion, surface tension, viscosity and inertia. Again, good schemes can demonstrate secondorder convergence toward the analytical solution with a small prefactor. See [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] for a comparison of different schemes.

More Complex Test Cases

Simple-looking test cases, for which analytical solutions exist, are often the most challenging, as illustrated by the history of spurious currents. More complex test cases are also useful however, in particular for assessing practical applicability of numerical schemes, including speed, robustness etc. An important issue for these tests is the availability of reference solutions: analytical solutions are usually not available, or have restrictions (e.g. on amplitudes, Reynolds numbers etc.) which can be difficult to enforce in numerical simulations ; experimental reference data can be available but error bars can be large and the experiments often include physical effects (e.g. surfactants, temperature gradients, compressibility etc.) which complicate their comparison with simpler numerical models. A popular example of this class is the case of rising bubbles, often used for validation of surface tension models. Due to a lack of accurate reference solutions, the validation is often qualitative, with a "visual" comparison of the shapes obtained experimentally or numerically. While this was useful when methods were inaccurate enough to cause obvious departure from the expected solutions (for example the extreme case of bubbles bursting due to spurious currents), this is insufficient to assess the relative accuracies of modern numerical methods. A useful approach, which requires substantial effort, is to provide accurate, converged, numerical reference solutions for non-trivial problems. For example [START_REF] Hysing | Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF], Featflow 2008) give reference solutions, using different numerical methods, for rising bubbles which can be reproduced accurately by other methods (Basilisk 2013, rising.c). A similar effort is made for Taylor bubble solutions by [START_REF] Marschall | Validation of interface capturing and tracking techniques with different surface tension treatments against a Taylor bubble benchmark problem[END_REF] and [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks[END_REF].

SELECTED APPLICATIONS

Numerical simulations are particularly useful in combination with laboratory experiments. Their advantages and drawbacks are often complementary so that simultaneous design of laboratory and numerical experiments can lead to deeper insight into complex physical phenomena. Figure 2 illustrates an example of this approach. A millimetre-size water droplet impacts on a pool and creates a complex splash structure. The top view is a zoom on the impact zone, seen as a vertical cut through the center of the drop. The axis of revolution is aligned with the left border of the image. The pool is coloured in blue and the droplet in red for visualization but there are only two fluids: water and air (light green). This configuration was studied in detail both numerically and experimentally by [START_REF] Thoraval | von Kármán vortex street within an impacting drop[END_REF]. The experiment and the numerical simulation are both very challenging due to the wide range of spatial scales and short duration of the phenomenon. Care was taken to ensure converged axisymmetric simulations. This required spatial resolutions of order 10 4 grid points per drop diameter i.e. resolved structures of order one micrometre. Besides an accurate surface tension model (well-balanced, height-function VOF-CSF), several numerical ingredients were necessary: adaptive mesh refinement, efficient multigrid pressure solver and parallelism [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Agbaglah | Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method[END_REF]. The numerical results were very consistent with experimental observations for the whole range of impact regimes (see also [START_REF] Agbaglah | Drop impact into a deep pool: vortex shedding and jet formation[END_REF] for an impressive comparison with high-speed X-ray imaging), but predicted a regime characterised by the unexpected von Kármán vortex street of Figure 2. This is associated with complex dynamics of the ejecta sheet which periodically entraps toroidal air bubbles. This regime was not observable using the side-view camera of the original experimental setup. This numerical result lead to the re-design of the experiment to use bottom-view cameras with the goal of observing the bubble rings predicted by the numerics [START_REF] Thoraval | Drop impact entrapment of bubble rings[END_REF]. A sample of the images obtained is given in Figure 2 bottom, for different pool depths. Although the subsequent three-dimensional breakup of the toroidal bubbles (some of them are still intact on the bottom-right frame) cannot be predicted by the axisymmetric simulations, the experimental results spectacularly confirmed the numerical discovery. The motion of gas bubbles in a liquid is a canonical example of the subtle balance between surface tension, viscous and buoyancy forces. The transitions between various regimes (straight, zig-zag or spiralling ascent) are particularly difficult to capture, either experimentally or numerically. They have been investigated numerically in a recent series of articles by [START_REF] Cano-Lozano | The use of volume of fluid technique to analyze multiphase flows: Specific case of bubble rising in still liquids[END_REF][START_REF] Ghidaglia | Capillary forces: A volume formulation[END_REF]. In contrast with the previous example, full three-dimensional simulations are necessary. The boundaries between regimes are controlled by the coupled interaction of the shape of the deformable bubble and the associated vorticity generation and wake formation. An example of the resulting trajectory, wake structure and bubble shapes is given in Figure 3. Accurate modelling of surface tension is vital to minimise spurious vorticity generation at the interface. As in the previous study, Cano-Lozano et al. were careful to check the numerical convergence of their results. This required a resolution of 128 grid points per bubble diameter. A very large tank of 8 × 8 × 128 diameters is necessary to be able to follow the bubble for a long time. This leads to formidable resolution requirements: 2 34 17 billion grid points on a regular grid! Adaptive mesh refinement brings this down to around 10 million grid points and make the simulations possible but still expensive (see http://basilisk.fr/src/examples/bubble.c for a full example). A large number of timesteps is necessary to capture the transition to established regime, in particular because of the explicit timestep restriction discussed previously. Note also that the parameters chosen correspond to those for a millimetric air bubble rising in a liquid roughly ten times more viscous than water. Further refinement would be necessary to properly capture the boundary layers for an air/water bubble. These two examples illustrate the capabilities, as well as limitations, of state-of-the-art models of surface tension. Obtaining numerically-converged results clearly requires considerable computing power. This is feasible, but challenging, for complex two-dimensional (or axisymmetric) configurations. Provided care is taken, very valuable insight can be gained from such simulations (see e.g. [START_REF] Samanta | A falling film down a slippery inclined plane[END_REF][START_REF] Fuster | Instability regimes in the primary breakup region of planar coflowing sheets[END_REF], Hoepffner et al. 2013[START_REF] Deike | Capillary effects on wave breaking[END_REF] for a small representative sample).

In three dimensions, only relatively simple configurations can be studied with confidence that results are fully independent from the numerics. That said, the situation was similar for two-dimensional simulations fifteen years ago, with the added limitation of less accurate surface tension models. Under-resolved three-dimensional simulations can still give very useful qualitative results for flows which are challenging to study experimentally, provided one controls the effect of resolution and checks consistency with the available experimental data and theoretical models. Representative examples of this approach include studies of atomisation [START_REF] Herrmann | Detailed numerical simulations of the primary atomization of a turbulent liquid jet in crossflow[END_REF], Desjardin et al. 2013[START_REF] Chen | High-fidelity simulations of impinging jet atomization[END_REF][START_REF] Jain | Secondary breakup of a drop at moderate Weber numbers[END_REF][START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF], industrial processes [START_REF] Mencinger | Numerical simulation of ligament-growth on a spinning wheel[END_REF] and waves [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]. A balanced trio of experimental, theoretical and numerical approaches can be extremely effective and I expect numerical models of surface tension to play an important role in future advances in our understanding of complex multiphase flows.

FUTURE ISSUES

1. None of the methods presented in this review satisfy both well-balancing and momentum conservation: properties which are required for consistency and robustness. Integral formulations, which have been relatively neglected, could be a promising research direction. 2. Although high-order height-function schemes have been demonstrated for curvature estimation, current volumetric formulations are formally first-order accurate.

This follows from detailed analysis of Peskin's scheme by [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF]. Although immersed interface schemes have been extended to second-order [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF][START_REF] Leveque | Immersed interface methods for Stokes flow with elastic boundaries or surface tension[END_REF][START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF], they have only been applied so far to Lagrangian interface discretisations and fluid-structure interactions. Their generalisation to generic two-phase flows, with an implicit interface representation, remains an open question. 3. Robust time-implicit schemes have already been formulated for Lagrangian immersed interface methods [START_REF] Mayo | An implicit numerical method for fluid dynamics problems with immersed elastic boundaries[END_REF][START_REF] Newren | Unconditionally stable discretizations of the immersed boundary equations[END_REF]) and there is no reason to believe that these results cannot be generalised to surface tension, however recent efforts in this direction have yielded somewhat confusing results [START_REF] Denner | Numerical time-step restrictions as a result of capillary waves[END_REF]. 4. Extensions to more than two phases, including the consistent treatment of triple points or lines has only been considered recently (see e.g. [START_REF] Li | Incompressible multiphase flow and encapsulation simulations using the moment-of-fluid method[END_REF]). 5. The schemes described in this review, and especially the height-function method, are most easily implemented on regular Cartesian grids (or their adaptive versions).

Their generalisation to unstructured grids, usually favoured for industrial applications, has so far concentrated on kinematics rather than dynamics. Surface tension schemes on these grids, often based on diffuse algebraic VOF formulations, are currently quite limited compared to the state-of-the-art on regular grids. 6. Finally, efforts should be pursued to provide standardised, benchmark cases showing numerical convergence (even if only first-order) for relevant, non-trivial physical configurations.

DISCLOSURE STATEMENT

The author is not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.

  Figure 1Approximations of the interface Heaviside and (normalised) Dirac functions (top and bottom figures respectively) for different volumetric methods.

Figure 2

 2 Figure 2 Top: Side view of an axisymmetric numerical simulation of the von Kármán vortex street created by the impact of a millimetric water droplet impacting on a pool. Bottom: high-speed experimental imaging of the bubble rings created by the vortex street for different pool depths. Adapted from Figures 10 and 11 of Thoraval et al. (2013).

Figure 3

 3 Figure 3 Numerical simulation of (a) trajectory, (b) vorticity distribution and (c) shapes of a gas bubble rising in a liquid in the "spiralling regime". The density ratio is close to air and water. The Galileo and Bond numbers are 100.25 and 10 respectively. Adapted from Figures 10 and 11 of Cano-Lozano et al. (2016).

Table 1

 1 Approximations of the interface Heaviside function for different volumetric methods.

S. Popinet

www.annualreviews.org • Numerical Models of Surface Tension

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

S. Popinet

ACKNOWLEDGMENTS