Computational interfacial rheology

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
Describing multiphase materials or multiphase material flows invariably involves making constitutive assumptions about the propensity of the interface to transmit stress [1][2][3]. Interfaces can be ''simple'' and be described by a constant interfacial tension. A first complexity arises when the concentration of surface-active species, which directly control interfacial tension, is changed as a function of time, necessitating a description of the transport phenomena, coupled to the deformation of the interface and the bulk and/or disturbance flows [4]. Further and intrinsic complexity arises when the interface becomes structured, due to an increased dense packing of those objects at the interface or due to lateral interactions between the surface-active moieties which can lead to (transient) network structures [5][6][7][8][9][10]. Such structured, or complex, soft-matter interfaces respond with extra stresses to a deformation, and linear and nonlinear viscoelastic or viscoplastic rheological properties then emerge which lead to the field of interfacial rheology [3,11,12].
Application areas are broad and diverse, ranging from aerated or multiphase food products [13][14][15][16][17], consumer care products [18], materials in biomedical applications as for example in lung surfactants [19][20][21], bacterial biofilms [22][23][24] or in drug delivery applications with monoclonal antibodies [25], and materials in the energy sector such as asphaltenes [26,27] and CO 2 storage [28] and finally polymeric materials (foams and blends) [29][30][31]. Goals of interfacial design are to stabilize products against coalescence [32] or Ostwald ripening [33], to control the morphology of disperse systems and polymeric blends N.O. Jaensson et al. COVID-19 pandemic. A prime example is the stability of the lungs and the lung collapse where a Laplace instability is suggested to occur, which is associated with acute respiratory distress syndrome [21]. Better understanding of these phenomena under dynamic conditions, crucial for improved treatments, would benefit from a computational modeling approach. The same holds true for better understanding of aerosol formation in breathing or sneezing, where the surfactant film in the lungs plays an important role in the production of submicron droplets [43]. Computational interfacial rheology could aid in predicting the size distribution of these droplets and help in designing mitigation strategies. There are many more problems where computational interfacial rheology could have significant impact; to whet the reader's appetite, some examples are given in Fig. 1.
The measurement of the equations of state and rheological material functions has a long history, going back to iconic scientists such as Ben Franklin, Lord Rayleigh, Agnes Pockels, Joseph Plateau and Irving Langmuir (see the book by Tanford for a historical overview [44]). Methods to measure the linear and non-linear rheological functions at constant surface area have been developed for shear and extensional deformations, as will be reviewed below [45,46], and experiments with clean kinematics and a rigorous separation of bulk and interfacial contributions have been developed, benchmarked and there is clarity on the operating windows where reliable and accurate data can be obtained [47]. The study of the rheological properties in compression or dilation is more challenging as the interplay between the changes in thermodynamic properties, due to changes in surface concentration, and intrinsic compressional viscoelasticity is difficult to deconvolute [48]. Many commercial instruments, such as the popular pendant drop instrument, analyze the response to deformations as if only a surface or interfacial tension is present. As a consequence, the intrinsic material functions are not readily obtained. This calls for computational interfacial rheometrical methods, combining such experiments with computational efforts to deconvolute the effects.
In our opinion, a ''perfect storm'' is brewing for generating even more broad advances using computational interfacial rheology and fluid mechanics. This review aims to provide firm footing for such progress, starting by laying out the basic assumptions commonly used in interfacial rheology in Section 2, followed by the theoretical description concerning how to set up the boundary conditions for rheologically complex interfaces in Section 3. Whereas the development of constitutive interfaces still allows for further progress, frame invariant descriptions of Newtonian, non-Newtonian, elastic and viscoelastic as well as basic viscoplastic equations are now available; the state of the art for these interfacial constitutive models is reviewed in Section 4. This is followed by a brief introduction to differential geometry in Section 5, which is crucial for the proper mathematical description of interfaces. Then, methods for interfacial rheometry, where numerical methods play an important role in extraction the rheological parameters of the interface, are reviewed in Section 6. Initial results on the simulation of deforming droplets and bubbles as well as drainage and thinning flows are discussed in Section 7, and we end this review with some open problems and an outlook in Section 8.

Basic assumptions
In interfacial rheology, the treatment of interfaces is typically done using a continuum approach, where the flow in the bulk of each liquid is described by conservation equations for mass, momentum and energy, and appropriate coupling conditions are employed at the liquid-liquid or liquid-fluid interface. In this ''sharp-interface'' framework, pioneered by Gibbs, the interface is considered a 2D dividing where ( , ) is the interfacial or surface tension, which is a state variable which depends only on the excess concentration and temperature , s = − is the surface unit tensor, where is the normal of unit length to the interface, and is the surface extra stress. For simplicity, in this review we only consider isothermal problems and spatial variations of the surface tension are thus only due to variations of the excess concentration. The surface stress tensor is a 2D secondorder tensor embedded in 3D space, and is symmetric and tangential, i.e., it maps tangential vectors to tangential vectors and it maps normal vectors to the zero vector [2]. Due to its symmetry, the surface stress tensor can be uniquely described by two orthogonal principal directions and two principal stresses, where the principal directions are tangential to the surface. Interfaces that have non-zero, possibly deviatoric, extra stresses will be referred to as ''complex''. In the literature, ''fully'' tangential tensors are sometimes defined separately from tangential tensors [55], which we will not do here, since the two definitions coincide for symmetric surface tensors.
To describe the surface extra stress, interfacial constitutive models are needed. Paraphrasing Pozrikidis [56], there are generally three approaches for deriving models that link surface stresses to deformations: (1) the interface is described as a 3D continuum, often assumed incompressible [57], and constitutive equations are derived in the limit of zero thickness, (2) the interface is described as a 3D continuum, and assumptions are made concerning the deformation and stresses in normal direction, possibly yielding non-tangential surface stresses, and (3) the thickness direction is disregarded from the outset, and the interface is described as a 2D surface embedded in 3D space. For many of the interfaces considered in interfacial rheology, it does not seem appropriate to describe interfaces as 3D continua since they are only a few molecules or particles in thickness. We therefore prefer the third approach, which also fits nicely within the framework of Gibbs. Moreover, we mainly treat models that split the constitutive response in an area-preserving part and a distortional part, since these offer much freedom for, e.g. fitting experimental data [58], and using experimentally observable material functions to predict the behavior during flow or deformation.
Using non-equilibrium thermodynamics, tangential slip has been identified as an entropy-producing mechanism for multi-phase systems [52], and may occur for large molecules [59,60]. However, the authors are not aware of reports of interfacial slip for small-molecular systems, suggesting that its magnitude is small enough to be neglected for practical purposes. Hence in most cases it can be assumed that the tangential components of the velocity are continuous across the interface. In addition, it assumed that there is no phase change across the interface, thus only ''material interfaces'' are considered. The latter assumption implies that the normal velocities are continuous as well, i.e. the velocity vector is continuous across the interface, and we therefore do not introduce a separate interfacial velocity. The 3D position vectors with respect to a chosen origin are denoted by , and the subset s contains those position vectors that point to the interface. A more formal definition of s will be given in Section 7.1.

Momentum balance for complex interfaces
We start by investigating the momentum balance for interfaces, by considering the conservation of the ''excess momentum density''. To describe the transport of the surface excess quantities in the sharpinterface framework, an expression for the interfacial conservation law is needed. However, a straightforward generalization from the bulk is complicated due to the possible motion of the interface through 3D space, and the idea of ''fixed point on the interface'' relying on the non-unique manner of how we choose to describe the location of the interface [51]. Following Slattery et al. we first introduce the concept of ''surface particles'' as infinitesimal particles that are bound to the interface, following its motion in 3D space [2]. It can be seen that, by definition, the velocity at the interface is the time derivative of the spatial location, keeping a surface particle constant. With computational interfacial rheology in mind, we follow an approach similar to the arbitrary Lagrangian Eulerian approach [61,62], where the interface is described by a grid that does not necessarily move with the surface particles. An interfacial grid velocity g is introduced, which is the time derivative of the spatial location, keeping a grid point constant, for which the only requirement is that the grid follows the shape of the interface. This requirement can be expressed as where we note that there are many choices possible for g , and a particular choice is generally motivated by the problem at hand [51].
We are now in the position to define an interfacial material derivative s ( )∕ as the time derivative from the ''point of view'' of an observer on a fixed surface particle, which yields [2,63]: N.O. Jaensson et al. Fig. 3. A surface element with contour . The unit vector normal to the surface is denoted by and the unit vector that is both tangential to the surface and normal to the contour is denoted by . The fluid below the interface, or in − direction, is called fluid 1, and is endowed with a bulk stress tensor 1 . The fluid above the interface, or in direction, is called fluid 2, and is endowed with a bulk stress tensor 2 . Inside the surface, the existence of a surface stress tensor s is assumed.
where the partial derivative s ( )∕ is for a constant grid point and where ∇ s = s ⋅ ∇ is the surface gradient operator (where ∇ is the gradient operator). The velocity − g , appearing in Eq. (3), is called the intrinsic velocity [2] or deformational velocity [52], and is tangential to the interface, ''unaware'' of the motion in 3D space.
To derive the momentum balance for interfaces, we proceed along the lines of Leal [1], but we retain the possibility of the interface having dynamical properties of its own and we consider changes in excess momentum. A sketch of an interfacial element is given in Fig. 3, where we define the normal vector and the tangential vector . Note that there are two different tangential vectors, but for brevity will denote them both by in this section. We assume that there is a bulk fluid present both below and above the interface, with bulk Cauchy stress tensor 1 and 2 , respectively. Moreover, as discussed above, we assume that there exists a tangential Cauchy surface stress tensor inside the interface, which we denote by s , and which is defined by s = s ⋅ , where s is the interfacial traction vector.
Collecting all the forces that act on the interfacial element and identifying them as sources of momentum, we obtain the following excess momentum balance equation [2]: Note that Eq. (4) contains a time derivative over the, possibly timedependent, surface element . Moreover, the second term on the right hand side is an integral over the contour . To rewrite Eq. (4) in a more useful form, we invoke the surface Reynolds transport theorem and the surface divergence theorem for tangential surface tensors to obtain [2,63]: where use was made of s s ∕ + s (∇ s ⋅ ) = 0 [53]. The term on the left hand side of Eq. (5) is associated with the inertia of the interface itself, which is typically small, or this term is not even present in the gauge fixed by s = 0, and it is generally neglected [47]. Moreover, since Eq. (5) is valid for an arbitrary surface element , we can remove the integrals. Substitution of the general expression for the surface stress tensor, as given in Eq. (1), yields after rewriting: where the first term on the right hand side represents Marangoni stresses, the second term represents the jump in normal stress due to the curvature of the interface (capillarity), and the last term is due to the extra surface stress (interfacial rheology). To explicitly show the dependence of Marangoni stresses on the equation of state and on spatial variations of the surface excess concentrations, we note that the first term on the right hand side of Eq. (6) can be rewritten for the isothermal case as ∇ s = ( ∕ )∇ s . Inspired by the nomenclature as introduced by Venerus and Öttinger [52], we will refer to Eq. (6) as the passive momentum balance, whereas Eq. (5) is the integral form of the active momentum balance, with the difference being that the active case involves time-derivatives, which arise due to the left hand side of Eq. (4), whereas the passive case is purely a jump-condition. Similar nomenclatures can be introduced for, e.g., the balance of excess mass or the balance of excess energy.
For deforming curved interfaces, the components of Eq. (6) are most conveniently described in curvilinear coordinates that are not necessarily orthogonal, as will be shown in Section 5. However, it is instructive to investigate Eq. (6) further for a planar interface using the standard Cartesian coordinates , and . It is assumed that the interface is located in the plane, yielding the tangential unit vectors and and the normal = . The velocity vector is given by = + + , where we note that is zero at the interface due to the planar assumption. Moreover, we assume that the fluids are Newtonian, thus for the th fluid the stress tensor reads = − +2 , where is the bulk pressure, is the bulk viscosity and = (∇ + (∇ ) )∕2 is the rate-of-deformation tensor. The stress balance in the tangential direction given by reads: where the notation | means that the term is evaluated at the interface, but on the side of the th fluid. The stress balance in direction is obtained by exchanging and in Eq. (7). The tangential stress balance shows that the discontinuity of the slope of the tangential velocity across the interface is due to a viscosity ratio, Marangoni stresses or interfacial rheological properties, and thus expresses how ''stresscarrying'' the interface is. In order to describe the interfacial stress tensor , constitutive models are needed that link (rate-of-)deformation of the interface to the surface stress, which will be topic of the next section.

Interfacial constitutive models
In this section we will review constitutive models for interfaces. The extra surface stresses described by each model are identified with the extra surface stress in Eq. (6).

Boussinesq-Scriven model
One of the earliest constitutive equations employed for interfaces was obtained by assuming a linear coupling between the interfacial rate-of-deformation and the surface stress [64,65]. The resulting ''Boussinesq-Scriven model'', which can be regarded as the surface equivalent of the compressible Newtonian model for bulk fluids, reads in coordinate-free form [2,66]: where s is the surface (or interfacial) shear viscosity, s is the surface (or interfacial) dilatation viscosity, is the velocity vector on the interface and d s is the deviatoric part of the rate-of-deformation tensor s = (∇ s ⋅ s + s ⋅ (∇ s ) )∕2, thus d s = s − tr( s ) s ∕2. The Boussinesq-Scriven model can also be derived by considering the fact that the excess entropy production must be positive, and then seeking for the simplest linear law which ensures this [52]. Note that for some interfaces, shape changes are much easier than area changes, e.g. biological membranes, which makes them effectively incompressible [67]. These systems can also be described by including the surface incompressibility condition ∇ s ⋅ = 0 as a constraint, in which case s becomes insignificant and assumes the role of a Lagrange multiplier [68]. Moreover, the presence of a surfactant can also strongly inhibit compression or dilation of the interface since spatial variations in surface excess concentration yield interfacial stress gradients, through the equation of state ( , ), that can oppose area changes [4]. Thus the surface incompressibility constraint has been used for surfactant-laden interfaces as well [69]. In this review, we generally allow for changes in area, unless otherwise noted.
In the literature, a number of confusions have arisen concerning the proper form of Eq. (8), which we would like to clarify here. The first is that the velocity vector as used in Eq. (8) should be the full 3D vector and not only the tangential projection as was done in e.g. [70][71][72][73][74]. It is easy to show that using only the tangential projection cannot be correct: if a sphere or cylinder is blown up by a purely normal velocity at the interface, the interface is clearly deforming and the surface stress should be non-zero. Note, that for planar interfaces, e.g. as studied by Elfring et al., the normal component of the velocity vector at the interface is assumed zero, so using the tangential velocity leads to the correct result [71].
The second point of confusion in literature concerns the projection operation which is applied on ∇ s and (∇ s ) in the definition of s . This operation is crucial since we define as a tangential tensor, thus the normal components of ∇ s and (∇ s ) must be removed before making a direct connection between and s [75]. However, this projection is absent in various works [42,57,[69][70][71][72][73]76,77]. Whereas for planar interfaces the projection is not necessary due to the velocity being a tangential vector (and it can thus be omitted) [69,71,76,77], for non-planar interfaces this needs to be considered. For certain geometries and situations including the projection might not always lead to significantly different results [42,78,79], but it is recommended to remove the normal components of ∇ s and (∇ s ) .
Again, it is instructive to investigate Eq. (8) further for a planar interface, as was done in Section 3 to emphasize the meaning of the different terms. Substitution of the Boussinesq-Scriven model in the tangential stress balance given by Eq. (7), yields: where we assumed 2 = 0 for brevity. Introducing the velocity scale , length scale , and characteristic surface excess concentration c , Eq.
where Ma = − c ( ∕ )∕( 1 ) is the Marangoni number where ( ∕ ) will generally depend on , but we treat it here as a constant for simplicity, Bq = s ∕( 1 ) is the Boussinesq number, = s ∕ s is the interfacial viscosity ratio. Asterisks denote non dimensional variables. Eq. (10) elegantly shows the intricate coupling between interfacial rheology, bulk rheology and transport phenomena. In addition, it shows how the bulk can be regarded as a sink for momentum, with the importance determined by Bq.

Generalized Newtonian interfaces
In the Boussinesq-Scriven model, as introduced in Eq. (8), it was assumed that the surface viscosities were constant, which has limited applicability for complex interfaces, where the viscosities are generally a function of deformation-rate [80][81][82][83][84] and/or surface excess concentration [4]. A straightforward improvement to the model can be made by assuming that the surface viscosities are a function of the invariants of rate-of-deformation tensor and of [66]. To maintain the split into an area-preserving part and a distortional part, we will make use of the first principal invariant of s , defined as s = tr( s ) = ∇ s ⋅ , and the second principal invariant of the deviatoric part of s , defined as d s = det s ( d s ), where det s is the surface determinant [2]. It is then assumed that the surface viscosity functions can be written as The viscosity functions in Eq. (11), are similar to the bulk case [85], such as a power-law or Carreau model, with shear-thinning or shearthickening behavior determined by the sign of the power-law exponent [84]. Similarly they could be expressed as functions of the invariants of the surface stress tensor. For a planar interface, a simple shear flow with ratėyields s = (̇∕2)( + ), where and are the orthonormal tangential unit tensors as introduced in Section 3. An isotropic flow with area rate of changėyields s = (̇∕2)( + ), where > 0 for surface dilatation and < 0 for surface compression. It can thus easily be verified that the first invariant of s is the area rate of change, whereas the second invariant of d s is related to the rate of shear. Using the invariants of s and d s , other non-Newtonian behavior, typically inspired by models for bulk rheology, can be carried over to interfaces in a rather straightforward manner [66,82]. To illustrate this, we present here, as an ansatz, the compressible Bingham model for viscoplastic interfaces, which is based on the bulk counterpart of this model [86][87][88], but maintaining the deviatoric-isotropic split. The deviatoric surface stress is given by: where is the magnitude of the deviatoric surface extra stress tensor and where sy is the shear yield stress. Note that | d s | and | d vp | can also be expressed as functions of the principal invariants introduced above. For the isotropic part of the surface stress we obtain: where iy is the isotropic yield stress, i.e. the yield stress in compression/dilatation. Using the notation from the previous paragraph, the model predicts in simple shear: where ( vp ) = ⋅ vp ⋅ , and similarly for the other components of vp . In dilatation/compression the model predicts: Similar to the bulk counterpart of the model, the yield stress can be a function of, e.g., a structural parameter [87] or the surface concentration [88]. Moreover, the yield criterion as used in Eq. (12) is a ''von Mises-like'' criterion for interfaces, whereas the yield criterion in Eq. (13) depends on isotropic contributions and thus differs from the von-Mises criterion [89]. Although the Bingham model for interfaces has been successfully used for area-preserving flows [82,90,91], we are not aware of any other experimental/numerical work that employs the Bingham model for interfaces with area-changing flows in the form given here. Although conceptually quite different from, e.g. the power law model, we categorized the interfacial Bingham model as a generalized Newtonian model since the effective viscosity is only a function of the instantaneous rate of deformation, albeit in a discontinuous manner. Thorough investigations into thermodynamic consistency, experimental validation and numerical stability, especially for area-changing flows, will be crucial for it to be as successful as its bulk counterpart [92].
Besides the deformation rate, the surface excess concentration can affect the value of the viscosity [93][94][95]. This dependency is typically described by an exponential relation between the viscosity and the surface pressure, where the latter is defined by = 0 − ( , ), with 0 the interfacial tension of the clean interface. Depending on the system under consideration, both -thinning and -thickening interfaces exist [4]. To describe the distribution of at the interface, an additional interfacial balance equation is needed, which can be readily derived in a similar manner as the momentum balance in Section 3. In the absence of reactions, and assuming Fick's law for diffusion at the interface, conservation of mass yields the following surface-convection-diffusion equation: where s is the surface diffusion coefficient and n is a source term that describes adsorption/desorption from and to the bulk. Often, the grid velocity is chosen equal to the normal velocity of surface particles, which yields [96]: where t = s ⋅ is the tangential interface velocity and n is the magnitude of the normal velocity n = ⋅ . Adsorption/desorption is mainly important for small-molecular surface active agents, that typically do not exhibit interfacial rheological behavior [95]. However, because changes in lead to changes in surface tension through the equation of state ( , ), a time-dependent interfacial tension can exist for these systems. Historically, this effect has often been described using a complex interfacial modulus [97] which is readily obtained using an oscillatory pendant drop instrument. We will not do that here, however, since the resulting parameters are not rheological material functions, but they are related to transport phenomena, as is clear from the fact that they are geometry dependent [98,99]. From a computational perspective, the material parameters that control this time dependence, and thus should be used as inputs in the modeling, are the bulk diffusivities and energy barriers for adsorption (interaction potentials). Rheological material functions should only relate stress to only strain or strain rate, and these are easily identified for structured, yet passive interfaces. Reporting apparent interfacial moduli, which arise due to time dependent properties and the diffusion to/from an active interface, in terms of properties of a passive interface have obfuscated the field of interfacial rheology and how to use any experimental data in a predictive modeling efforts. Further discussion on this topic can be found elsewhere [3,46,100].

Elastic interfaces
Elastic interfacial constitutive models can be constructed that couple the surface stress to the deformation with respect to some reference state. Although interfacial structure could in principle lead to anisotropic interfacial elasticity, for simplicity only isotropic elasticity is considered here, implying that the surface stress and surface strain are coaxial [57]. We introduce the stretch ratio , which describes the change in length of a line s , tangential to the surface in the deformed state, relative to a line s , tangential to the surface in the reference state: = ( s ⋅ s ) 1∕2 ∕( s ⋅ s ) 1∕2 , where s and s are the position vectors to the interface in the deformed and reference configuration, respectively. For small deformations, the linear surface strain tensor can be used to construct constitutive equations, which is defined by s = (∇ s ⋅ s + s ⋅ (∇ s ) )∕2, where is the 3D displacement vector on the interface. For small deformations, the principal values of s are 1 = 1 −1 and 2 = 2 −1, where 1 and 2 are the two in-plane principal linear strains and 1 and 2 are the two in-plane principal stretch ratios. For linear elastic behavior, we can define a strain-energy area density function of the form [57] where s is the interfacial shear modulus and s is the interfacial dilatational modulus. By taking the derivative of H with respect to the linear strain, we obtain the stresses acting on an infinitesimal patch of the interface, similar to the approach in bulk [101]. The resulting stresses can be written in tensor form to yield Hooke's law for linear elastic interfaces [76]: where d s is the deviatoric part of linear strain tensor: d s = s − tr( s ) s ∕2. Alternatively, Eq. (19) can be written in terms of the interfacial Lamé constants, Young's modulus or Poisson ratio [48,75,102]. Note, that Eq. (19) is sometimes written with additional 1∕ terms, which arise due to the area change between the deformed and reference state [103]. However, since Eq. (19) derives from the theory of linear elasticity, it is only valid for infinitesimal deformations and, to first order, no distinction is made between the deformed and reference state. To go beyond the small deformation limit, finite strain tensors are needed.
In the literature, isochoric ''Mooney-Rivlin'' or ''neo-Hookean'' models are often derived by assuming that the membrane has a finite thickness and adding the assumption that the membrane deformation is volume-preserving [104]. The strain-energy volume density function, together with the isochoric condition 1 2 3 = 1, where 3 is the stretch ratio in the normal direction, yields hyperelastic interfacial models. However, these models generally have limited usefulness, mainly because they lack an independent parameter to describe the constitutive response to dilatational/compression of the interface [57]. Moreover, the assumption of a continuum in the thickness direction is questionable for systems that only consist of one to a few molecules or particles in thickness direction [56].
Alternatively, the thickness direction can be neglected completely, and the membrane elasticity can be formulated in terms of a strainenergy area density function. Here, we give two of such models which are particularly useful since they both employ a split in the isotropic and shear response. Both models can be written in terms of the left Cauchy-Green interfacial strain tensor s = s s , where s is the 2D interfacial deformation gradient tensor given by s = s ∕ s . Note that s is symmetric and tangential to the surface in the deformed state. Furthermore, is the relative area deformation, which can be written as = det s ( s ) = 1 2 . The first is the Evans-Skalak model, which is often used for the modeling of bio-membranes, and for which the strain-energy area density function with respect to the reference state is given in terms of the principal stretch ratios as [105]: where = 1 2 is the relative area deformation. More recently, a new constitutive equation was derived by Pepicelli et al. [58], with the notable difference that a Hencky strain measure is used for the isotropic deformation, which is known to be able to describe moderate non-linearities and has the additional advantages that it is symmetric in dilation/compression and additive [48]: From the strain area density functions, the surface stresses can readily be calculated [58]. For the strain area density function as given by Eq. (21) this yields the neo-Hookean model for compressible interfaces, which reads in coordinate-free form [58]: wherēd s =̄s − tr(̄s) s ∕2 is the deviatoric part of the left-Cauchy-Green interfacial strain tensor at constant area, which is given bȳs = s ∕ . For the Evans-Skalak model as given by Eq. (20), the deviatoric part of the surface stress is similar to Eq. (22), but the isotropic part is given by s ( − 1) s . A comparison of different strain measures for describing the isotropic deformation of the interface is shown in Fig. 4.

Viscoelastic interfaces
Similar to bulk rheology, linear viscoelasticity of interfaces can be introduced by integrating the rate of deformation over all past times multiplied with decaying relaxation moduli [106]. The deviatoric and isotropic parts of the linear viscoelastic surface stress can then be written as [107]: where d lve denotes the deviatoric part of the linear viscoelastic surface stress and where ( ) and ( ) are the relaxation moduli for shear and dilation/compression, respectively. The constitutive relations given by Eq. (23) are the general linear viscoelastic model for interfaces, which can be used to, e.g., define the complex interfacial viscosity of the interface, in a similar way as their bulk counterparts [108]. The relaxation moduli can be estimated using the Maxwell model as exponentially decaying functions: where 0 and 0 are the instantaneous shear and bulk moduli, respectively, and s and d are the relaxation times in shear and dilation/compression respectively. Extension of this model to multiple relaxation processes is straightforward, and can be done by including multiple modes [106]. Similar to bulk rheology, this general linear viscoelastic model can be used to describe small-strain experiments such as creep recovery and stress relaxation after step strains. The latter approach was used to verify the Boltzmann superposition principle, one of the pillars of linear viscoelasticity, for a polymer multilayer at the air-water interface under isotropic compression, where the extra stress was more important compared to the thermodynamic one [109] (see Fig. 5). Other constitutive models, based on linear viscoelasticity, such as the linear Kelvin-Voigt and Jeffreys models, can be readily derived [107], and also simple thixotropy responses have been described in this framework [110]. However, it is important to emphasize that the model as given by Eq. (23) is not material frame indifferent and should only be used at small strains. A straightforward way to obtain non-linear viscoelastic solid interfacial constitutive models, which are valid for large deformations, is the superposition of the Boussinesq-Scriven model as given by Eq. (8), and a hyperelastic model for the interface, such as the one given by Eq. (22), yielding a quasi-linear Kelvin-Voigt type model for viscoelastic solid interfaces. This type of modeling was successfully applied in the description of biomembranes [111][112][113], and more recently to the surface viscoelastic behavior of graphene oxidesheets at the airwater interface [77,114]. An advantage of this approach is that it can be readily implemented in simulations [42,72,115,116]. Alternatively, non-linear viscoelastic liquid interfacial constitutive models can be derived following typical approaches from bulk-rheology, leading, e.g., to the surface upper-convected Maxwell (SUCM) model, which describes the evolution of the deviatoric viscoelastic stress, denoted by d ve , by the following evolution equation [53,107,117]: where the surface upper convected derivative is defined as Note that the surface upper convected derivative, as we define it here, is different from the one found often in the literature, which includes additional projection operations [53,107,117,118]. In Appendix, we show that the form including the additional projections is not frame invariant, whereas the form given in Eq. (26) is frame invariant, thus we believe the latter to be correct. Moreover, the third term on the left hand side of Eq. (25) is absent in the original presentation of the SUCM model [107,117], but it is added to ensure that d ve remains traceless [53,[119][120][121]. The evolution of the isotropic part of ve can be written as [107]: where s is the surface (or interfacial) dilatation viscosity. The SUCM model presented here is in a convenient form due to the split between the deviatoric and isotropic parts. However, similar to the case of bulk viscoelasticity [120], other variants of the SUCM could be constructed as well, that treat the isotropic part in different ways. Moreover, extensions can be made to this model to include, e.g., shear-thinning, yielding the interfacial Giesekus model [53,107]. Despite their success for bulk rheology, there is, as far as the authors know, no experimental literature using these non-linear viscoelastic models, and numerical work is scarce and often uses the incorrect form of the SUCM model [118].

Interfaces with resistance to bending
Bending stresses are important for biological membranes [122], but recently have also been shown to play a role for particle-laden interfaces, since they can have non-zero bending moduli, as observed in e.g. graphene oxide sheets [114], microparticles [123] or nanoparticles [124]. To describe the resistance against bending, the starting point is typically the Helfrich Hamiltonian, which assumes that the bending energy is quadratic in the mean curvature minus the spontaneous curvature and linear in the Gaussian curvature [125]: where and G are the bending moduli, m = ∇ s ⋅ ∕2 = ( 1 + 2 )∕2 is the mean curvature, where 1 and 2 are the principal curvatures, 0 is the mean spontaneous curvature and G = 1 2 is the Gaussian curvature. Variational procedures lead to a force density of the form [126,127] where it was assumed that no topological changes occur, thus, due to the Gauss-Bonnet theorem, the second term in Eq. (28) is a constant which does not lead to a force density. The tangential term in Eq. (29) is often neglected, but, for e.g. red blood cells in flow, its magnitude is small compared to the tangential elastic stress, or it can be absorbed in the surface pressure if the interface is incompressible [126]. The force density as given by Eq. (29) can be readily incorporated in numerical simulations, i.e. it is added to Eq. (6) as a jump in traction [126,128,129]. Moreover, b can be written as the surface divergence of a surface stress tensor b , which can be added to the surface stress [57,122,130]. We note, however, that care should be taken when using the latter approach within our framework, since the bending stress tensor is not tangential (i.e. b ⋅ ≠ ), which has consequences for, e.g., the surface divergence theorem [2], and thus the force balance given in Eq. (6). Since bending stresses act normal to the interface, this model can be readily combined with one of the other proposed constitutive models for tangential interfacial stress.

Differential geometry basics
For interfaces that do not deform, often a coordinate system can be chosen where iso-surfaces of one of the coordinates correspond to the shape of the interface, i.e. Cartesian coordinates for planar interfaces or spherical coordinates for spherical drops/bubbles. However, for interfaces that deform, a general description in terms of curvilinear coordinates is crucial, which will be introduced in this section. A more extensive introduction to the broad field of differential geometry can be found elsewhere [2,66,131,132]. Unless otherwise noted, we will make use of the Einstein summation convention, where Latin symbols ( , , etc.) are used to indicate that the index runs until two, whereas Greek symbols ( , , etc.) are used to indicate that the index runs until three.
A general description of interfaces can be obtained by a pair of coordinates, 1 and 2 , which provide a unique mapping of each point on the interface to its position in 3D space (see Fig. 6(left)): where s is the position vector for each point on the interface. In Fig. 6(right), we show that the curvilinear coordinates can be identified with the parametric coordinates, as used as in e.g. FEM, which will be shown useful in a later section.
The natural basis vectors are obtained by varying only one coordinate, while keeping the other constant: where we should note that 1 and 2 form a basis that is not necessarily orthogonal, nor are their lengths necessarily normalized. The dual basis vectors are defined by where is the Kronecker delta. Any vector or tensor that is tangential to the surface, denoted by and respectively, can thus be written as where and are the contravariant and covariant components of the vector , respectively, and and are the contravariant and covariant components of the tensor , respectively. Moreover, the surface gradient operator can be written as And finally, we can write for the surface unit tensor where the metric matrix = ⋅ and the dual metric matrix = ⋅ were introduced. To describe the velocity and position of the interface s , we make use of the standard Cartesian basis, thus = and s = , where the index is identified with the Cartesian coordinates , and . For example the surface gradient of the velocity can be written as [133]: and for the surface divergence of the velocity: The expressions in Eqs. (37) and (38) are particularly useful since they show how the surface operations can be written in terms of the Cartesian components of the velocity, natural/dual basis vectors and position vector, and can thus be incorporated in any standard numerical method for 3D bulk flow.

Orthogonal curvilinear coordinates
We conclude our brief treatment of differential geometry by investigating the special, but very important case of orthogonal curvilinear coordinates, where we have by definition: It is often more convenient to normalize the natural basis vectors: where ℎ = ( ⋅ ) −1∕2 (no summation) are the scale factors, also known as the metrical coefficients [66]. Note that the normalized natural basis vectors coincide with the normalized dual basis vectors, and we thus do not distinguish between contravariant and covariant components when using orthonormal basis vectors. The set of unit tangential vectors 1 and 2 and unit normal vector forms a useful basis for interfacial rheology, and the resulting components of vectors/tensors are referred to as the physical components [2]. Moreover, the surface gradient operator can be written as: which allows one to write the surface derivative operations in a more convenient form [66]. The scale factors for standard coordinate systems such as cylindrical-and spherical coordinates are well-documented, see e.g. [108,131,134], and are useful for cases where the interface does not deform and coincides with one of the coordinate surfaces, as we will show in the next section. Clearly, for the Cartesian coordinate system, all the scale factors are one, and the basis vectors are constant in space.   6. The curvilinear coordinates 1 and 2 , e.g. as used in FEM, provide a one-to-one mapping from real space (left) to parametric space (right).

Methods for interfacial rheometry
Due to the inherent hydrodynamic coupling of the interface with the bulk, as shown by Eq. (6), it is generally not possible to obtain directly a universal rheometric flow at the interface, i.e. a flow with a known surface (rate-of-)deformation tensor. The bulk fluids act as a sink for the interfacial momentum, with the ratio between the bulk and interfacial transport depending exactly on the a priori unknown interfacial rheological properties. Analytical and numerical methods are therefore crucial for calculating the flow field in a given geometry, and thus the velocity field at the interface. With the deformation rates at the interface known, rheological properties of the interface can be inferred by investigating the response of the probe, e.g. displacement or force/torque [46]. A large range of experimental setups is available for interfacial rheometry, a complete overview of which is beyond the scope of the current review. Extensive reviews are available in the literature [3,46,135]. The methods in this section have in common that they assume that the interface retains a fixed shape, allowing for a description in orthonormal coordinates, as introduced in Section 5. Moreover, by benefitting from symmetries occurring naturally in the problem, the governing equations and numerical methods for solving them are simplified considerably.

Interfacial shear flows
Several experimental setups are available for the generation of an interfacial shear flow, without the effects of compression/dilatation; here we will discuss a number of them. Among the first proposed geometries, were ones where an interfacial shear flow was generated by lowering a disk or bicone geometry into an interface, and rotating the cup that held the fluids [136]. Torques on the disk could be measured using, e.g., a torsion wire. For these geometries, Oh and Slattery [136] applied the Boussinesq-Scriven model for the surface stress and the Stokes equation for the bulk flow, neglecting effects from inertia. An implicit solution was obtained for the velocity field, and an iterative numerical scheme was applied to calculate the interfacial viscosity based on the measured torque on the probe in steady shear. This approach was later extended for sinusoidal deformations [137], and linear surface viscoelasticity can be investigated in the same geometries by [138] and [139], by assuming that the viscosity in the Boussinesq-Scriven model was complex. Erni et al. extended the analysis of Oh and Slattery, and proposed a bicone fixture that could be used in standard rotational rheometers [140]. Moreover, they introduced an improved numerical scheme that allowed for the calculation of the velocity field for a given torque, which can be used to infer interfacial properties, but fluid inertia was still neglected. A recent analysis of the bicone geometry was provided by Tajuelo et al., who used a finitedifference scheme to numerically calculate the flow in the bulk and at the interface and which does take inertia into account [141]. The same group later made the code, called BiconeDrag, freely available [142].
However, accurate measurements using the disk or bicone geometries are limited to interfaces with high interfacial viscosities as compared to the bulk viscosity [47]. This becomes immediately clear when defining a macroscopic Boussinesq number as [46,143,144]: where is a characteristic velocity, s and b are the characteristic length scales at which the velocity decays at the interface and in the bulk phase, respectively, s is the contact perimeter between the rheological probe and the interface, and b is the contact area between the probe and the bulk phase. The parameter = s b ∕( b s ) has units of length, and shows how the Boussinesq number can be effectively increased by, e.g., increasing the perimeter of the probe or decreasing its surface area in contact with the bulk. A geometry with a smaller value of , thus an increased sensitivity, is the double wall ring (DWR) setup, which was proposed by Vandebril et al. [145]. In their analysis, the Navier-Stokes equation was solved in the bulk and coupled to the flow at the interface, using finite-differencing schemes on 2D meshes, which was made possible by the radial symmetry of the problem. This approach included effects of inertia, and the response of the interface was described using the Boussinesq-Scriven model with real or complex viscosities. Interfacial viscosities could be inferred by matching boundary conditions between the simulations and the experiments in an iterative manner. The numerical codes were made freely available [145]. Another approach for decreasing , thus increasing the sensitivity, is using a ''knife edge geometry'', and similar numerical approaches to take the bulk flow into account for this geometry, typically using finite differencing methods, are published in the literature [83,146,147].
Another geometry for shear interfacial rheology consists of a (micro)needle at an interface, which is moved by a known magnetic force, as proposed by Brooks et al. [143]. With this setup, which is commonly called the interfacial shear rheometer (ISR), the interfacial viscosity can be calculated from the displacement of the needle, for example for sinusoidal oscillations of the force. Brooks et al. assumed that Bq ≫ 1, such that a linear velocity profile is generated at the surface, i.e. all the dissipation occurs in the interface. To allow for measurements at lower values of Bq, Reynaert et al. proposed a numerical calculation of the flow field for linear-elastic interfaces, which takes the full flow in the bulk, including effects of inertia, into account [148]. In this manner, measurements at lower values of Bq became accessible, also for the ISR. This approach was later verified by Verwijlen et al. using a point force solution at the interface [144]. The accuracy of the ISR can be improved by decreasing the needle size, thus decreasing the value of [149]. Finally, we note that more advanced numerical approaches have become available for the ISR, such as a boundary element method [150], who reported results in good agreement with previously mentioned works (see Fig. 7).
The high aspect ratio of the (micro)needle used in the ISR ensures that end effects play a negligible role. Disk-shaped interfacial probes  can be used to generate an interfacial shear flow by rotation inside the interface, for example by magnetically applying a torque to a ferromagnetic ''microbutton'' [93,[151][152][153][154]. Translation of probes with an aspect ratio of (1) at planar interfaces induces a mixed flow field, which makes the analysis significantly more difficult; these will be discussed in Section 6.3.
Rheometers that probe different deformation modes, typically inspired by bulk rheometers, have also been proposed. For example, the Cambridge Interfacial Tensiometer [155,156] is a 2D equivalent of a filament stretching rheometer. Although this device would ideally operate at constant surface area, the numerical/experimental analysis by Verwijlen et al. showed that the interfacial flow field contains dilational, shear and extensional components, depending on the ratio of the material parameters [157]. They conclude that these type of ''mixed flows'' contain a wealth of rheological information, but are typically difficult to analyze, limiting their accuracy. In a follow-up work, similar conclusions were drawn for the interfacial rheometer consisting of a ring-like fixture for the DWR with undulations that produces a more controlled mixed flow field ''by design'' [158]. A typical feature of such a mixed flows, is that the flow field is dominated by the ''cheapest'' deformation mode, i.e. the one for which the viscosity is lowest, which inherently leads to a low sensitivity [158].

Area-changing flows
Whereas the experimental and numerical methods for flows at constant surface area are well-established, the development of dilational/compressional rheometry is not as far. The main difficulty lies in obtaining a truly isotropic deformation of the interface, without any effects of shear [46]. Historically, the drop-shape analysis of a pendant drop is the technique that is most often used to investigate the response of the interface to dilation/compression. Note that the same analysis also applies to rising bubbles, but for brevity we will only talk about pendant drops in this section. With the pendant drop setup, the surface tension is obtained by fitting the Young-Laplace equation to the shape of the drop, which can be easily obtained from Eq. (6) by assuming only hydrostatic contributions to the bulk stress and in the absence of deviatoric surface stresses. Therefore, the Young-Laplace equation assumes an isotropic and constant surface stress, and small capillary-and inertial relaxation times such that the shape is in equilibrium. Moreover, the deformation of the drop due to gravity must be large enough in order to obtain an accurate value of the surface tension [159]. By oscillating the volume of the drop, the interface deforms in a complex, mixed manner which includes shear and dilation/compression, as shown in the FEM simulations by Balemans et al. [42] and reproduced in Fig. 8. The surface deformation yields a change in the surface excess concentration , and thus the surface tension through the equation of state ( , ). This effect is commonly called Gibbs elasticity, but its origin is generally different then, e.g., the origin of hyper-and viscoelasticity as introduced in Sections Section 4, which arise due to lateral interactions between surface-active moieties [48]. The response of the interfacial tension for systems with small molecular surfactants is often time-dependent, which is due to the relaxation processes as visible in Eq. (16): in plane convection/diffusion and transport normal to the interface [160]. Since these processes are geometry-dependent [98,99], the time-dependency measured for these type of systems is due to transport phenomena, as opposed to being a true material parameter [3].
To make drop-shape analysis suitable for complex interfaces, one should (1) ensure that the measured quantity is indeed a material property, e.g. by repeating the measurement with different drop sizes, and (2) allow for anisotropic surface stresses, the occurrence of which is now well-established by different groups [48,109,[161][162][163][164][165][166][167]. Estimating the surface stresses for complex interfaces can be aided greatly by including a measurement of the pressure inside the drop and obtaining the curvature locally, instead of a global fit of the Young-Laplace equation [48,166]. With the surface stresses known, iterative numerical procedures can be used to find the constitutive parameters of the interface, an approach generally called elastometry [48]. Alternatively, the constitutive equations can be applied directly to obtain the material parameters from the shape [103,164,167], which is a more straightforward, but seems to be a less accurate approach, as discussed in detail by Nagel et al. [48].
A closely related setup for measuring dilational properties of interfaces is the capillary pressure tensiometer (CPT), which combines a measurement of the radius of small drops, typically 10 to 100 μm, with a pressure measurement [168,169]. With the pressure inside the drop known, accurate measurements of the isotropic surface stress can be performed without the need for deformation due to gravity. In fact, smaller drops are preferred in the CPT, due to the higher signal to noise ratio of the pressure measurement, and the shorter capillaryand inertial relaxation times, allowing for measurements at higher frequencies [170]. In addition, the effect of shear becomes smaller for smaller drops, which have a shape closer to spherical [171], which can be exploited in the CPT to obtain interfacial deformations that are close to isotropic. As shown by the theoretical analysis by Kotula and Anna, the CPT can indeed be applied to complex interfaces in the presence of extra stresses [172]. Although the CPT is a large step forward in obtaining more isotropic deformations of interfaces, it can only be applied to interfaces of high curvature, due to the constraint of a spherical droplet shape.
Recently, a new experimental setup was proposed, called the radial through, which can apply truly isotropic deformations to a planar interface by moving an elastic band that is held in place by 12 ''fingers'' [58]. Due to the isotropic deformation, the stress state will be isotropic even for complex interfaces. This setup was later successfully applied to investigate surface viscoelasticity of a polymer multilayer at an airwater interface, and values were obtained that were consistent with  measurements using the rising bubble setup [109]. Moreover, it was shown for the first time that for these type of complex interfaces, concepts from bulk linear rheology, such as the Boltzmann superposition principle, indeed carry over to interfacial rheology (see Fig. 5). A miniaturized, 3D printable, version of the radial trough was developed by Kale et al., which allows for simultaneous dilatational deformation and interfacial microscopy [173].

Particle translating at a planar interface
Particles at complex interfaces are important for a variety of reasons. For example, they can be used to improve the stability of emulsions/foams [33], and the diffusion of particles at interfaces can give insights into the rheological properties of the interface [135,174]. Moreover, they can be used as idealized model system for the understanding of, e.g. the diffusion of proteins in biomembranes [175]. In the limit of Bq → ∞, the hydrodynamic problem of a particle translating at a viscous interface in an unbounded domain is closely related to the Stokes' paradox, which states that there is no creeping flow solution for a cylinder of infinite length translating perpendicular to its axis, which is due to the logarithmic decay of the disturbance velocity field [175,176]. Indeed, the slow decay of the velocity autocorrelation function leads to the inability of defining a surface viscosity, or diffusion coefficient, for such systems [177]. In the seminal work of Saffman and Delbrück, the logarithmic decay is cut short by including the subphases as ''sinks for momentum'', and the mobility of a disk of radius at an incompressible, viscous planar interface of finite thickness is derived [178]. The solution given by Saffman and Delbrück is valid for = ∕ s ≪ 1, where is the bulk viscosity and s is the surface shear viscosity [179]. This model was later extended to all values of [180], and for finite depths of the subphase [179].
Although the assumption of surface incompressibility seems to hold for many types of interfaces, e.g. biomembranes, this assumption does not hold in general for interfaces with other types of surface-active moieties, see e.g. [58]. Moreover, particles protruding into the bulk phases might not be accurately approximated by a flat disk embedded in the interface. Danov et al. presented a numerical model for a spherical particle embedded in a planar interface with both effects from shear and dilation/compression present, but in the absence of Marangoni forces [181]. In a later work, effects from Marangoni forces were included up to first order, but without the coupling to flow in the bulk [182]. In the work of Fischer et al., a spheroid at an incompressible interface between two infinite bulk domains is considered [69]. They rationalize the surface incompressibility of the interface by assuming that the Marangoni effects are large. This, assumption, however, might break down for interfaces at low surface coverage, or with a weak slope with respect to in the equation of state ( , ).
A particle translating through an interface will under most circumstances generate a complex disturbance velocity field. This was assessed experimentally by using needles of different aspect ratios and a using broad range of relevant non-dimensional numbers (Ma, Bq and the ratio between shear and dilatational s ∕ s ) [183]. Drag coefficients for a probe in an incompressible interface are not expected to depend on probe aspect ratio, suggesting effects of compressibility or Marangoni stress related effects play a role. Elfring et al. investigated theoretically this interplay between interfacial rheology, Marangoni stresses and coupling with the bulk assuming a thin subphase, allowing for a lubrication approach [71], see Fig. 9. They reached the same conclusion that a particle translating at an interface is not always ideal for characterizing the shear rheological properties interfaces, due to the possible occurrence of a mixed flow field (dilation/compression and shear) and the tangled effects of Marangoni stresses and interfacial rheology. Methods with clean kinematics are recommended to determine interfacial shear and dilational properties, as outlined in the previous sections.
Since a full overview of the field of the hydrodynamics of particles at interfaces is beyond the scope of the current review, we will conclude this section by referring to two recent reviews for further reading [4,135].

Wrinkle analysis
We conclude this section with methods aimed at recovering the bending modulus of the interface, which are based typically based on the analysis of wrinkles occurring during compression. Wrinkling (also called ''buckling'') occurs since, for low enough bending moduli, out-of-plane deformations can be a cheaper mode of deformation than in-plane compression. The resulting wave-length of the wrinkles can be used as a direct measure of the bending resistance of the interface. This analysis is typically done by assuming uniaxial compression of the interface, which allows for a 1D treatment, i.e. assuming invariance in the direction perpendicular to the compression, which is superimposed on a planar or curved shape. Starting point is the beam equation, which is given by [123]: where is the bending modulus, ℎ is the vertical deflection of the interface, is the compressive force per unit length, is the density difference between the subphases and is the acceleration due to  gravity. Eq. (43) can also be derived as a linearized form using the Helfrich bending energy as introduced in Section 4.5 [184]. Moreover, the surface stress can be calculated using the in-plane elastic constitutive models as introduced in Section 4.3 [164,184]. The beam equation has a sinusoidal solution for ℎ, with the wave length given by = ( 16 4 ∕( )) 1∕4 [123]. The bending modulus can be further expressed as function of bulk linear elastic parameters and a third-order dependency on the thickness of the interface [123], but we will not do that here due to the concerns of regarding the interface as a continuum in thickness direction, as raised in Section 2 of this review. From the wave length of wrinkles as observed in experiments, the bending modulus can be calculated. This approach has been successfully applied to graphene oxide sheets [114], and a number of complex interfaces in the pendant drop/rising bubble setup [103,164]. Non-linear calculations, based on the Helfrich model, are presented by [184], yielding a range of interesting deformation modes of the interface.
Finally we note that thermal fluctuations of vesicles are often used to probe the bending resistance of membranes [185][186][187]. Moreover, active particles can be included inside vesicles to act as local forces, inducing ''active fluctuations'' and yielding a large range of shapes that do not exist in equilibrium systems, as shown recently by Vutukuri et al. [188]. The large deformations, induced by the local forces, allow for an even more complete probing of the membrane mechanics, including the role of membrane (visco)elasticity and local changes in lipid concentration.

Flows with deforming complex interfaces
The approach of computational interfacial rheology follows the same lines as computational bulk rheology [189]: constitutive models are tested and calibrated in simple flows, as shown in the previous section, which can then be used to predict and benchmark the flow in increasingly complex flow scenarios. The latter is the topic of the current section of this review. For computational interfacial rheology, the complexity mainly arises due to the deformation of the interface, possibly leading to highly curved interfaces, thin fluid films and breakup and coalescence events. Since our focus is on continuum models, many of the numerical techniques available in the literature are based on the ''classical'' methods for bulk computational rheology, such as the finite difference method (FDM) [72], finite volume method (FVM) [190] and the finite element method (FEM) [42]. Moreover, the boundary element method (BEM), coupled to e.g. FDM or FEM for the surface computations, has been successfully applied as well [133,191].

Numerical approaches
Three main issues need to be addressed for the simulation of multiphase flows with deforming, rheologically complex interfaces: (1) describing the location of the interface, (2) evaluating the surface differential operators as introduced in Section 5 and (3) coupling the interfacial constitutive equation to the bulk flow. These issues are generally addressed by either interface capturing (IC) or interface tracking (IT) methods (see Fig. 10), which we will briefly review here.
In IC, the location is described in an implicit manner, typically by introducing a continuous function = ( , ), which is defined in the entire 3D domain, and the location of the interface is, by definition, the zero levelset of this function: ( s ) = 0. Assuming that remains constant for a fixed material point yields in the Eulerian description: which can evaluated in the entire domain using standard methods. The normal to the surface can be obtained from = ∇ ∕|∇ |, which can in turn be used to evaluate the surface operations from the definitions s = − and ∇ s = s ⋅ ∇. Coupling to the bulk solution of, e.g. the velocity, can be done by approximating the interface with a small thickness, effectively converting the surface operations to volume operations, and using the model in the sharp-interface limit [192,193]. Alternatively, the interface can be reconstructed by a triangulation of the zero-levelset, and approaches similar to the extended finite element method (XFEM), where the solution space is enriched by allowing for discontinuities within elements, can be applied [194][195][196]. A major advantage of IC methods is that all equations can be solved on stationary Eulerian meshes, and that break-up and coalescence events are handled relatively easily, but they are generally less accurate than IT methods [197]. Variants of the IC method that have been used for interfaces with surface viscosity are the volume-of-fluid (VOF) method, where the phases are described by volume fraction in the cells [190], and phasefield/diffuse-interface methods, that typically have a sound physical basis to describe the interface with a finite thickness [193,198].
In the IT method, a surface mesh, i.e. a 2D mesh embedded in 3D space, is introduced that describes the interface. Both structured and unstructured 2D meshes, of varying interpolation order, have been used in the literature; a comprehensive overview is given by Barthès-Biesel [57]. The surface mesh is identified with the grid as introduced in Section 3, and the mesh velocity is thus identified with g . The normal velocity of the mesh is determined by Eq. (2), but the choice of the tangential velocity is arbitrary, which can be used to, e.g., make sure that the elements of the surface mesh remain evenly spaces and/or do not deform significantly [199]. Within the IT framework, evaluation of the surface differential operations is relatively straightforward due to the availability of the surface mesh, and can be done by e.g., identifying the curvilinear coordinates with the reference coordinates of the element as shown in Eqs. (37) and (38), and in Fig. 6. The coupling to the bulk depends on how the volume is discretized. Often, an ''immersed boundary'' approach is used, where the surface mesh is embedded in a stationary Eulerian grid, and forces are calculated in the nodes close to the interface, e.g. by using a smoothed delta function [72].  This approach has the advantage that no boundary-fitted meshes have to be generated, but the solution close to the interface, i.e. within one volume element, is not accurate. Alternatively, meshes can be generated that are boundary fitted to the interface, such as the one shown in Fig. 10(left), yielding very accurate solution close to the interface, but requiring regular remeshing [42]. Lastly, a combination between the FEM for the solution on the surface mesh and the BEM for the solution in the bulk have been proposed, with the advantage that the full flow field in the bulk is obtained, without the need for a volume mesh [133].

Drops and bubbles in shear flow
An important goal of computational interfacial rheology is the simulation of drops and bubbles, endowed with complex interfaces, in flow. Due to the coupling between the flow field, interfacial properties and transport equations, the governing equations allow for analytical solutions only under the simplest circumstances. For full 3D flows, large deformations and/or interactions between multiple drops and bubbles, numerical approaches are invaluable. Although simulations exist for drops and bubbles with complex interfaces in several flow scenarios, such as the oscillating pendant drop [42,116], freely rising bubbles [195] and Poiseuille flow [195,200], we will focus here on simple shear flow as it serves as a stepping stone for more complex flows.
The dynamics of drops and bubbles with rheological interfaces in simple shear flow can give insight into the bulk rheology and stability of emulsions and foams, and is therefore relevant from both fundamental and practical point of view. In addition, understanding the flow conditions for which breakup of drops/bubbles occurs is important for predicting, e.g. the droplet size distribution after processing of the material. Finally, the rheometric bulk flow conditions allow for comparisons to experiments on suspended drops with soft-matter interfaces [201,202], e.g. by investigating their equilibrium shape.
The transient dynamics of a single, isolated droplet depend on the shear rate and the material parameters of the drop, the suspending fluid and the interface. Generally, three scenarios are possible: the drop can (1) elongate and break up into smaller drops, (2) exhibit periodic deformations that can persist for long times or even indefinitely, or (3) attain a steady-state shape with the interface continuously moving around the drop in a process called ''tank-treading'' [203]. If a steadystate is reached, the drop typically takes on an ellipsoidal shape for smaller deformations, with the major axis tilted with respect to the shear flow direction, described by the inclination angle . The deformation of the drop is quantified by the Taylor deformation parameter = ( − )∕( + ), where and are the major and minor axes of the ellipse, respectively. For larger deformations, or for example in the case of confined shear flow, the steady-state shape of the drop significantly deviates from an ellipsoidal shape and it is often useful to use the length of the drop to quantify the amount of deformation [204].
As early as the 1950s, Oldroyd derived an expression for the viscosity of a dilute emulsion of spherical drops endowed with viscous interfaces [205]. This analysis was later extended by Danov to include the effects of Marangoni stresses and Gibbs elasticity [206]. The influence of surface viscosity on the deformation of an isolated drop in shear flow was investigated by Flumerfelt using a first-order perturbation analysis, valid for small deformations [207], and the first full numerical simulations were presented by Pozrikidis, which were, however, limited to surface viscosity ratios of unity [191]. A complete 3D numerical analysis of the steady state of drops with viscous interfaces, including the separate effects of shear and dilational viscosity, was presented by Gounley et al. who used the BEM, coupled to the FEM, to investigate the behavior of drops endowed with Boussinesq-Scriven interfaces in shear flow [133], see Fig. 11. Their analysis revealed that increases with surface shear viscosity, whereas it decreases with surface dilational viscosity. Since an increase in can lead to droplet breakup, this analysis can have implications for the stability of foams and emulsions during flow, e.g. material processing. Similar conclusion were reached by Luo et al., who investigated the coupled effect of a -dependent surface viscosity and Marangoni stresses, assuming an insoluble surfactant [8]. They conclude that the coupling between surfactant transport, interfacial rheology and Marangoni stresses yield non-trivial and non-negligible alternations to the droplet dynamics.
Numerical simulations of droplets encapsulated by a (visco)elastic membrane in simple shear have been mainly motivated by research into biological systems, with the goal of mimicking e.g. red blood cells during flow. The systems under consideration can be broadly categorized as (1) vesicles, where the membrane consists of a ''liquid-like'' incompressible lipid bilayer with bending elasticity, (2) red blood cells, which are vesicles with an additional protein network attached to the bilayer to give it surface shear elasticity and (3) capsules, which is a generic term often used for droplets encased in a synthetic (e.g. polymer) compressible membrane [208]. Although biological systems are outside the scope of this review, the overlap with computational interfacial rheology is significant, so we will give a very brief overview here. More extensive reviews for the flow of vesicles and capsules are available in the literature [57,208].
In a series of papers [209][210][211], Pozrikidis and coworkers investigated the simulated capsules and vesicles in shear flow using a boundary element method (BEM) [212]. Following this work, different groups have successfully performed simulations of vesicles and capsules in shear flow using the BEM [68,[213][214][215][216][217][218]. Moreover, other numerical methods have been proposed as well, such as the phase-field method [219], the level-set method [220], spectral BEM [221], boundary fitted FEM [222] and immersed boundary FEM [223,224]. The majority of the cited literature above focuses on either purely viscous or purely elastic membranes, whereas membrane viscoelasticity has received much less attention. Viscoelasticity of the membrane has been incorporated in simulations mainly by using the quasi-linear Kelvin-Voigt model as introduced in Section 4.4, thus by assuming that the surface   [133]. Source: Reproduced from [133] with permission from Cambridge University Press. stress is given by a superposition of a viscous surface stress, as given by the Boussinesq-Scriven model, and an appropriate hyperelastic surface stress. Diaz et al. used BEM simulations to investigate the behavior of a capsule with a viscoelastic membrane in an elongation flow [225]. The influence of membrane viscoelasticity of a capsule in shear flow was investigated by Yazdani and Bagchi [115]. A similar approach was taken by Gounley and Peng, using an immersed lattice Boltzmann method [226], which was also recently used for the simulation of red blood cells [227]. Benchmark solutions for a drop, endowed with a Boussinesq-Scriven or Kelvin-Voigt interface, in shear flow were presented by Carrozza et al. [228]. Carrozza et al. extensively validated their boundary-fitted FEM implementation using the method of manufactured solutions, and showed higher order spatial and temporal convergence.
Immersed FDM approaches that include viscoelastic membranes have been introduced recently by Li and Zhang [72,73]. In many of the above-mentioned works, it has been observed that a direct implementation of the Kelvin-Voigt model is unstable and therefore a standard linear solid (SLS) model is often used for the viscoelastic surface stress [72,115,226,227]. The SLS model can be interpreted as an elastic element in parallel with a viscoelastic Maxwell element, and thus by choosing the modulus in the Maxwell element large, the Kelvin-Voigt model is recovered. Finally, we note that particle-based method have been used as well for the simulation of viscoelastic membranes, where the membrane is modeled as a network of viscoelastic springs [229].
We conclude this section by briefly mentioning that capsules and vesicles in simple shear can exhibit buckling/wrinkling instabilities due to compressive stresses, similar to the buckling as introduced in Section 6.4. This phenomenon has been as observed experimentally [230] and in simulations [115,214], and more information on this topic can be found elsewhere [231].

Drainage and thinning flows
Drainage and thinning flows are important since they play a major role in the stability of emulsions and foams [46,232] and coating operations [233,234]. Due to the coupling between capillarity, transport phenomena and interfacial rheology, it is often difficult to predict drainage times for all but the simplest systems. Moreover, due to an inconsistent separation of material-and transport parameters, especially for area-changing flows, a clear link between stability and interfacial rheology has been established for only a limited number of systems, e.g. [33,235]. Numerical and analytical approaches are thus crucial in elucidating this complex problem. Often, drainage problems are understood in terms of mobility of the interface, which is a measure for the magnitude of the tangential surface velocity. An alternative interpretation, which is particularly useful for rheologically complex interfaces, is to consider the ability of the interface to ''carry stress'', as expressed by Eq. (7). Interfaces with solely surface tension effects can also be ''rigidified'' by Marangoni stresses, as shown by Champougny et al. [233], and reproduced in Fig. 12. Since the films that are considered are generally much thinner than their lateral length scale, lubrication approaches can be applied which simplify the numerical The dashed line is the limiting case of a rigid interface, which is approached as Ma is increased, i.e. the interface is ''rigidified''. See [233] for definitions of the symbols used here. Used with permission of Royal Society of Chemistry, from [233]; permission conveyed through Copyright Clearance Center, Inc. methods considerably [1]. This approach allows for rewriting the governing equations in terms of the unknown height of the film, and by assuming radial symmetry, solutions can be sought on 1D meshes, often using simple FDM or FEM approaches.
The effect of interfacial viscosity and Marangoni stresses on the drainage between two drops was studied by Zapryanov et al., assuming the interface remained planar, thus the problem could be solved in a single point [236]. They concluded that, within the range of parameters included in their work, Marangoni stresses play a more important role in immobilizing the interface than interfacial rheology. Li studied the coalescence between two spherical bubbles being forced together, where the surface stress is given by the Boussinesq-Scriven model and a first-order effect of the Marangoni stresses is included [237]. Deformation of the interfaces was allowed, and dimples were observed, significantly altering the drainage behavior. Results were obtained that could explain experimental observation, but the analysis assumed no deviatoric stress inside the bubbles, and can thus not be extended to the drainage between drops. Chatzigiannakis and Vermant recently defined criteria how film rupture depends on the hydrodynamic conditions where a competition between the characteristic timescales of drainage and rupture has been experimentally determined. [238,239]. Comprehensive numerical modeling is needed here that accounts for drop deformation dynamics, film drainage until the point of rupture including the random evolution of propagating fluctuations. Already for systems with a constant interfacial tension or surfactant transport this has been shown to be extremely challenging [240,241]. In a series of recent papers, the drainage between drops was studied by Ozan and Jakobsen, whose boundary integral approach does solve for the flow inside the drops [242]. They considered the coupled effect of surfactant transport, Marangoni stresses and a surface excess concentration-dependent surface viscosity [243]. Finally, they included non-linear surface viscoelasticity based on the SUCM model [118], which is, to the best of our knowledge, the first numerical work that includes a full solution of the SUCM model for drainage between drops. However, the simulations by Ozan and Jakobsen are based on the incorrect SUCM model, which is not frame invariant and does not preserve the zero trace of the deviatoric tensor, as explained in Section 4.4.
Numerical analyses of other problems have been performed in a similar fashion. In the work of Hermans et al., a hemispherical dome was risen against into an air-water interface laden with different types of lung-surfactant replacements and the drainage of the film was investigated [54]. From numerical solutions of the drainage equation, which include the effects from both Marangoni stresses and surface rheology, experimental observations could be rationalized. Similar numerical approaches have been taken to describe the stability of antibubbles [244,245] and the classical Landau-Levich-Derjaguin (LLD) dip coating problem [234]

Open problems and outlook
From this review it is evident that quantitatively describing multiphase materials and in particular the role of a complex interface between the different components is a challenging task from both a computational and rheological point of view. Some of the challenges are quite similar as in the bulk flow of complex fluids; predictive constitutive models are essential that are valid not only in well-defined deformations such as shear, extension and compression/dilatation, but also in mixed flows. In this review several prototypical flows were discussed with a strong interplay between bulk flow and interfacial properties and despite all progress in literature, many challenges still need to be addressed. In particular, in systems where all length and time scales of bulk and interfaces are coupled and need to be resolved simultaneously. For example in drop-drop interaction up to five to six orders of magnitude in length scales are fully coupled and coalescence time scales are determined by the drainage of the liquid in the film, which is strongly dependent on interfacial transport and rheology. Even for well-designed smooth problems this is a very challenging (but also very relevant) problem to solve. Advanced computational methods that are accurate and efficient, such as full monolithic solvers, need to be developed for solving stiff numerical problems where all degrees of freedom are coupled. Also from a mathematical perspectives there are many open questions. Coupling different bulk and interfacial discretizations might give rise to (finite element) compatibility constraints and call might for stabilization techniques similar to the ones used to resolve the high Weissenberg number problem in bulk, such as the log-conformation approach [246,247]. Computational challenges also arise when interfacial phenomena or films reach a length scale well below those where continuum mechanics is still assumed to be valid. One could think of linking continuum bulk descriptions with molecular dynamics models for interfaces.
In addition to computational challenges also from a rheological point of view the scientific field is still wide open. For example in many industrial processes temperature plays an essential role which calls for experimental, theoretical and numerical methods to deal with temperature-dependent interfacial constitutive behavior. If coefficients in models, like a dilatational viscosity or modulus, are made temperature dependent one could argue if these type of constitutive relations are thermodynamically consistent. Other challenges are found when electrostatic charges become relevant or systems with a phase change across and at the interface. Also polymer crystallization at liquid-liquid interfaces is a process that simultaneously is influenced by interface chemistry and dynamics, and rheology.
To conclude, much progress has already been made in the field of computational interfacial rheology, but also many opportunities for further progress are present. An abundance in problems of both societal and industrial relevance could potentially benefit from this relatively young and exciting field, which offers a whole range of interesting numerical and theoretical challenges, waiting to be tackled.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
where s = (∇ s ) is the interfacial velocity gradient tensor and Ω s =̇⋅ s ⋅ describes the rate of rotation, where( ) has been used as a short-hand notation for the surface material derivative, i.e. the time derivative in a frame that translates with a surface particle. An alternative expression for the surface unit tensor in the rotated frame is * s = * − * * . Using the relations given in Eqs. (45) to (49), and using symmetry of the surface stress tensor, it can be shown thaṫ * thus the surface upper convected derivative given in Eq. (26) is frame invariant. Likewise, it can be shown thaṫ * which shows that the surface upper convected derivative appearing in, e.g. [53,107,117,118], is in fact not frame invariant.