A unified numerical model for wetting of soft substrates

The wetting of deformable elastic structures has been recently shown to comprise a rich variety of new physical phenomena (stick‐slip motion, durotaxis, Shuttleworth effect, etc.) whose fundamental understanding demands for numerical simulation tools. In this article, we develop a novel unified model and numerical method for this problem. As special features, the method includes exact incompressibility and a linear monolithic assembly of the Navier‐Stokes and Cahn‐Hilliard equations which stabilizes dominant surface tension at small interface lengths. Also, solid viscosity is included which offers the opportunity to simulate the surfing behavior of droplets on viscoelastic substrates for the first time. We show that the method is highly accurate and more robust than most previous approaches and illustrate its potential by numerical simulation examples.


INTRODUCTION
Wetting of flexible substrates (also: soft wetting, or elasto-capillarity) plays a major role in a broad variety of phenomena. Technical applications include the patterning of cells 1 or droplets 2 onto soft surfaces, the optimization of condensation processes, 3 or the deposition of droplets in ink-jet printing or additive manufacturing. Also in biology, capillary interactions with flexible materials play a major role in self-organization of cell tissues, 4 cell motility, 5 and cancer cell migration. 6 Still, our understanding of the dynamics of such soft wetting processes lags by far behind of what is known about rigid surfaces. 7 Despite significant progress in computational modeling of (de)wetting of rigid substrates and the interaction of single fluids with elastic solids separately, the continuum modeling and simulation of wetting of elastic substrates has remained essentially unexplored. 8 One reason for this discrepancy is certainly the challenging nature of the problem which is inherently multi-scale in space (processes at the micro-scale determine dynamics on the macro-scale) and time (different time scales ranging from microseconds to seconds). Additionally, the strong coupling of hydrodynamics and substrate dynamics near the contact line requires special attention from a numerical point of view. Another reason for the lack of computational results in this field lies in the diversity of the required numerical techniques: the simulation of wetting typically employs interface capturing techniques (like the phase-field method 9,10 ), while the interaction of a fluid with an elastic structure (fluid-structure interaction, FSI) is typically modeled by interface tracking techniques, in particular the arbitrary Lagrangian-Eulerian (ALE) method. 11 Accordingly, first simulations of wetting of elastic substrates appeared only very recently. 8,[12][13][14][15] These references comprise basically the whole available literature on numerical simulations of soft wetting.
The phase-field model (also diffuse-interface model) is probably the most popular continuum model to describe wetting phenomena with moving contact lines. The idea is based on an order parameter, the phase-field function , which is used to indicate the fluid phases. The phase-field function varies smoothly across the interface leading to a (thin) diffuse interface which is also physically motivated by the composition of real fluid-fluid interfaces. 16 This diffusive nature of the interface regularizes the stress singularity at the contact line, making it a very natural approach for moving contact lines. The rigorous thermodynamic substructure of the model allows for consistent modeling of topological transitions (e.g., Reference 17), and enables energy stable discrete formulations 10 and robust time discretizations. [18][19][20] This makes phase-field models physically more sound than sharp interface models, like level-set or ALE, which typically need to introduce experimental fitting parameters 7 while offering an almost similar accuracy for two-phase flow simulations. 21,22 Consistent boundary conditions for the phase-field modeling of wetting have been derived in References 23 and 24. In Reference 9, an improved phase-field model was derived from variational arguments and showed impressive quantitative agreement to molecular dynamics simulations.
Accordingly, the available approaches for simulation of soft wetting 8,[12][13][14][15] all use a combination of different phase field methods for two-phase flow coupled to an ALE discretization of the fluid-solid interaction. For the latter part, advanced standard FSI methods are used in all works. For example, computations on the elastic structure are carried out in a stationary reference domain and mappings are used to transform the corresponding quantities to real space. While this approach is useful to precompute the stiffness matrix in standard FSI, it looses this advantage in soft wetting, because adaptively changing grids are required to resolve the region near the contact line. Hence, one could avoid the use of mappings in the ALE discretization by actually moving the grids. Also, the solid incompressibility, which is given for basically any relevant soft material, is not exploited in any current numerical method.
In this article, we aim to exploit these properties to come up with a unified and simpler mathematical model for soft wetting. The model is amenable to a simpler numerical treatment, yet we will show that the unified approach makes it numerically more robust with respect to time step restricions than previous explicit models. 8,[12][13][14] As solid viscosity can be easily included, we will consider soft wetting of a viscoelastic material for the first time. We will illustrate the accuracy and efficiency in a numerical benchmark case and show the potential of the method by the first simulations of the surfing behavior of droplets on viscoelastic substrate.
Next, in Section 2, we introduce the model equations and reformulate them in a unified way. In Section 3, we present the numerical method. Benchmark tests and further illustrative simulations are carried out in Section 4. Conclusions are drawn in Section 5.

Configuration
We consider a two-phase fluid domain Ω f ⊂ R d , which consists of both, a liquid-and an ambient phase and is connected to ). An illustration is given in Figure 1, where the example of a drop on a substrate is shown. In order to indicate the fluid phases, an order parameter is introduced, referred to as phasefield. We define = 0 in the ambient phase (liquid or gas) and = 1 in the liquid phase. Across the ambient-liquid interface, the phasefield varies smoothly, following a tangent hyperbolic profile. This leads to a thin diffuse interface of width . 25 Furthermore, we denote Γ as the sharp fluid-solid interface. Three surface tensions are present along the three contact lines. A fluid-fluid tension contracts the interface described by the phase-field and two solid-fluid tensions, 0 and 1 , contract the surface of the solid material which is in contact with the ambient and liquid phase, respectively.

Binary fluid model
If variable (phase-dependent) densities are involved, different formulations of the phase field model equations have been proposed. [26][27][28][29][30] However, in real scenarios away from the critical point, the involved interface thickness is small such that F I G U R E 1 Illustration of the soft wetting scenario. A deformable solid domain Ω s borders on a two-phase fluid domain Ω f . The two fluids are indicated by the value of the phase field function . The fluid-fluid interface is diffuse with a thickness and carries a surface tension . The fluid-solid interface Γ has two distinct surface tensions 0 and 1 depending on the contacting fluid all these models yield comparable results. 21 Therefore, we use the simplest model from Reference 26 which considers a volume-averaged velocity formulation with fluid velocity v f and reads with viscous, pressure, and capillary stress Here, the following denotations were used: the phase-dependent fluid density f , constant mobility m, the chemical potential , the double well potential W, the (scaled) fluid-fluid surface tensioñ, and the phase-dependent viscosity f . In our case, the double well potential is chosen as W( ) = 2 (1 − ) 2 , which yields a scaling of a physically motivated surface tension according tõ= 3

√
2 , see, for example, Reference 21. Compared to the previous literature on soft wetting, 15 we omit the term I(̃∕2 |∇ | 2 +̃∕ W( ) − ) in the definition of the capillary stress S ca , as it only rescales the interfacial pressure. That is, the pressure p f used here is related to the real physical pressure p phys by p f = p phys −̃∕2 |∇ | 2 +̃∕ W( ) − . Omitting the left out terms also from the interfacial force balance (Equation (5) 2 below) ensures that this approach is consistent.

Solid elasticity
We intend to formulate and solve the equations of solid elasticity in the domain Ω s . That approach is different from typical approaches in FSI in which the equations are defined on a stationary reference domain and mappings are used to transfer operators to the real physical domain. Furthermore, we will exploit the fact that all soft elastic materials (soft wetting occurs typically at Young's moduli in the range of few kPa), are essentially gels, which are incompressible, due to high fluid content. In most work on soft wetting (a notable exception being Reference 15), this property is approximated by using a Poisson ratio close to 1/2. Instead we will here use the exact incompressiblity condition and formulate the governing equations accordingly very similar to the Navier-Stokes equations. Using a linear elastic model, the equations of momentum and mass conservation are with viscous, pressure, and elastic stress Here, s , s , and G denote the solid density, solid viscosity, and elastic shear modulus, respectively. Furthermore, u is a displacement field, that is, the difference between current coordinates and initial coordinates of material points. This displacement field can be tracked either by an additional PDE, t u + v s ⋅ ∇u = v in Ω s , or by moving the grid with the velocity v s while memorizing initial coordinatesx of each grid point, such that u = x −x. We use the latter approach here. Note, that we include the nonlinear Euler-Almansi term in the elastic stress S el in Equation (4) to make the model geometrically nonlinear. Contrary to previous approaches for numerical simulation of soft wetting, we include a viscosity in the solid material, which leads to an effective description as a Kelvin-Voigt viscoelastic material. The viscosity of the solid is relevant for many dynamic phenomena in soft wetting, and its influence is hardly understood. 31 Moreover, let us note that in general G can be chosen as an arbitrary time-and space-dependent function in Ω S , for example, in order to implement an elasticity gradient in the substrate (for example to describe durotaxis 2 ). However, throughout this work we assume G to be constant.

Coupling conditions
The following coupling conditions apply at the solid-fluid interface Γ: The first equation is the continuity of velocities across the interface which is equal to the usual no slip condition. This condition also implies a topological condition which ensures that Ω f and Ω S are contiguous at the interface (see Section 3.1). In typical sharp interface formulations such a no-slip condition results in a pinning of the fluid/fluid contact line to the solid boundary. This effect, also known as contact line singularity, is overcome in diffuse interface models. In these models, the fluid/fluid interface can move not only by advection, but also by the Cahn-Hilliard diffusion which is constructed from energetic arguments, such as to minimize the overall surface energy (i.e., to realize the correct contact angles). Hence, the use of a phase field model for the fluid/fluid interface is an elegant and physically motivated way to avoid the contact line singularity.
The second equation is the dynamic condition which describes the balance of forces across the interface (namely traction and capillary stress). Here, n is the outer normal to Ω f . The solid surface tension s ( ) represents the tension along Γ between the solid and the fluid indicated by the value of the phase field. We use the following differentiable function of the phasefield s ( ) = ( 1 − 0 ) 2 (3 − 2 ) + 0 , which implies s (0) = 0 and s (1) = 1 . Further, ∇ Γ ⋅ denotes the surface divergence operator and P = I − n ⊗ n is the surface projection operator. Note that the present formulation coincides with the boundary conditions given in References 8,12, and 15 as ∇ Γ ⋅ ( s ( )P) = s ( ) n + ∇ Γ s ( ) with = −∇ Γ ⋅ n being the total curvature.
The third equation, Equation (5) 3 , is responsible for the formation of a dynamic contact angle, where is an inverse relaxation time. This condition has been derived in a variational model for moving contact lines. 9 Note, that when tends to infinity, the equation reduces to the static contact angle conditioñn ⋅ ∇ = − ′ s ( ). Let us note that both the kinematic and the topological condition are comparable to those in conventional FSI problems. However, the dynamic condition as well as Equation (5) 3 differ from conventional FSI problems in terms of contribution of the capillary stress S ca and the fluid-solid surface tension s , leading to a discontinuity of the traction across the interface (see also sec. 2.3 of Reference 15).
Finally, Equation (5) 4 ensures mass conservation (no penetration) on Γ. With this, we have specified a closed system of equations (up to outer boundary conditions) for the unknown variables v f , p f , v s , p s , , .

Unified model
We now unify the previously described equations of binary-FSI similar to the approach presented in Reference 32. The first step is to exploit the continuity of velocity (5) 1 to introduce the common velocity field Furthermore we denote the material derivative of a quantity by • t = t + v ⋅ ∇. The fact that the domain Ω s is not a reference domain, but moves with the material, makes it possible to define a common connected computational domain Ω ∶= Ω s ∪ Ω f . In this domain, we introduce the characteristic functions f and s to distinguish the domain parts, that is, = 1 in Ω for ∈ {f , s} and zero anywhere else. The surface Γ is identified with the surface Dirac delta function Γ . With this notation, we can add up the variational forms of the equations to end up with Here, = (x, (x)) and = (x, (x)) are the space-dependent density and viscosity. For both variables, we use an interpolation of the form: where f ,1 and f ,0 denote the density in the liquid and in the ambient phase, respectively. Let us remark here that the last two case distinctions are necessary to limit to physically meaningful values. The reason for this is that the phase field does not satisfy a maximum principle, hence its values can slightly exceed the interval [0, 1], which can create unphysical density, especially for large differences between f ,0 and f ,1 . This also applies in an analogous manner to . Equation (6) consolidates Equations (1) 1 -(1) 2 , (2), (3), and (4). Additionally also the kinematic and dynamic coupling conditions equations (5) 1 -(5) 2 are included. Hence, we end up with a single Navier-Stokes equation which incorporates momentum and mass balance of fluid and solid phases together with the interfacial balance of forces in a monolithic way. It is therefore expected that this formulation leads to superior stability of the numerical method, while greatly simplifying the numerical treatment of the system, as we will illustrate in Section 4.1.
For completeness let us note here that these equations are still accompanied by the Cahn-Hilliard system in the fluid domain: equipped with the boundary conditions (5) 3 -(5) 4 .

ALE discretization
We use the ALE method to discretize the numerical grids. The two domains Ω s and Ω f are discretized on two separate but connected, moving numerical grids. Within the elastic structure as well as on Γ, grid points move with the material velocity v, that is, they are material points. In contrast, the grid points in the fluid structure move with a continuous harmonic extension of this velocity in order to keep a proper shape of the mesh. In this work, the grid velocity w is calculated in the fluid domain by solving the Bi-Laplace problem and w = v in Ω s . Alternatively, other methods can also be used to extend the grid velocity into the fluid domain, for example, by penalizing length changes of triangles/tetrahedra, see Reference 33 for more details. The calculated grid velocity w is then subtracted in all convective terms of the governing equations. Therefore, the material derivate • t is replaced by where t,x * defines the time derivative along a quantity on a moving grid point.

Time discretization
The problem is discretized with equidistant time steps of size . At each time step, the general workflow of the numerical solution procedure is as follows. First, the coupled system of momentum balance, mass balance, and phase-field evolution, Equations (6)- (7) is solved in one monolithic system. Afterward the grid velocity w is computed as explained in the previous section. Then in a last step, the grid is updated, that is, each grid point is moved by the corresponding value of w.
An IMEX (implicit/explicit) Euler method is used to formulate a time discretization of Equations (6)- (7) which is linear in the solution variables. We denote quantities on discrete time points by a superscript, where () n − 1 refers to the previous time step and () n refers to the current time step. Explicit discretizations (i.e., step n − 1 ) are used for the nonlinear elastic term and most occurrences of the phase field in the momentum balance. The linear elastic term is taken implicitly, G(∇ u n + (∇ u n ) T ), by using the identity u n = u n−1 + v n (see also Reference 32). This way, we obtain an explicit part G(∇ u n − 1 + (∇ u n − 1 ) T ) and an implicit part G(∇v n + (∇v n ) T ), which effectively leads to an increment of the solid viscosity by G.
Additionally, stability is increased by an implicit treatment of the fluid-fluid surface tension force. The capillary stress S ca is taken semi-implicitly by using̃∇ n ⊗ ∇ n−1 . Note that we define the involved tensor product such that . Let us remark that the alternative discetizatioñ∇ n−1 ⊗ ∇ n does not lead to a stabilizing effect, since only the first occurrence of ∇ effectively contains the interface curvature, see Reference 18 for a more detailed discussion of stabilization of large surface tension in phase field models. This treatment must be accompanied by using the new velocity field v n in the advective term of the phase field. Further, the nonlinear derivative of the double well potential W ′ ( ) is linearized by a first order Taylor expansion. We end up with the time discrete system: In each time step n, find v n , p n f , p n s , n , n such that Here v n−1 * and n−1 * denote the corresponding quantities from the last time step, but after the applied change of the grid point coordinates due to the mesh update. Further, n−1 = (x, n−1 (x)), n−1 = (x, n−1 (x)).

Space discretization
In our finite-element approach, we consider the separated triangulations T f of Ω f and T s of Ω s , whereby the membrane Γ that connects the two domains is triangulated by T Γ = T f ∩ T s . In particular, the connection is ensured by the fact that every grid point on the interface of Ω f has a corresponding grid point on the interface of Ω s , sharing the same point coordinates. This allows us to enforce continuity of velocity across Γ, as well as to implement the jump conditions of the stress there. We present an example for the numerical mesh in Figure 2. We now present the combined discrete system of Equations (1)-(5) in the weak form. Adopting the approach from Reference 33, we introduce the following finite element spaces: Here, V h is the finite element space for the components of the velocity v. It ensures continuity of the respective variables across T Γ . Note that in Equation (11) 2 , denotes a placeholder for the distinction between the solid and the fluid domain. The use of two separate spaces, P h, f and P h, s , is motivated by the discontinuity of pressure across Γ caused by the solid surface tension. The use of standard finite element spaces for the discretization of the discontinuous pressure leads to poor numerical properties, with an approximation order of only ( √ h) w.r.t. to the L 2 norm. 34 Accordingly, the usage of two separate spaces, P h, f , P h, s , extends the standard Taylor-Hood finite element space by additional degrees of freedom of the pressure at the interface, such that the discontinuity can be exactly resolved. The remaining finite element space C h refers to and in Equations (1) 3 and (1) 4 .

Navier-Stokes system
With the previous arguments, we can establish a uniform weak formulation of the momentum and mass conservation equation for the combined domain Ω. We assume the density in Ω f to be constant in the following. The weak form reads:

Cahn-Hilliard system
We assume a static contact angle condition in this work, Equation (5)  Find Let us note that the Cahn-Hilliard system is still directly coupled to the Navier-Stokes system by the presence of n in Equation (12) and the presence of v n in Equation (14). This coupling, which removes time step restrictions for small interface length, 18 is linear. Hence, the coupled system can be solved by assembling Equations (12)-(15) in a monolithic way, and no subiterations are needed.

Validation study
We commence the verification and validation of the proposed model with a popular wetting test case: the relaxation of an initially half-spherical drop on a substrate. To follow up this benchmark case, we aim to reproduce the numerical reference results in Reference 15 and the experimental data presented in Reference 35. The computational setup and parameters are chosen as in Reference 15, but repeated here for completeness. As indicated in Figure 2, we consider a cylindrical computational domain with radius 350 μm and height 300 μm. The bottom of the domains consists of a 50 μm thick soft elastic substrate. A liquid droplet is placed in the center of the substrate. The shape of the drop is initially a spherical cap with radius R C ≈ 177.8 μm, centered slightly above the substrate surface, such that the fluid-ambient contact line assumes the equilibrium contact angle on the flat substrate. Note that this alignment of the correct initial contact angle is not a necessity of the numerical algorithm but is used to comply with Reference 15. For an initially undeformed planar substrate, the contact angle is given by Young's relation which applies only to flat rigid substrates. In particular the time-dependent order parameter (t, x) initially reads where the center x C of the spherical cap is placed such that its distance to the surface of the substrate is |R C cos |. 36 As in Reference 15, we choose the following surface tensions: which leads to an equilibrium contact angle of approximately 96.24 • and a distance of x C to the substrate of R C | 0 − 1 | ≈ 19.3 μm.
The shear modulus of the soft substrate has been reported to be G = 1 kPa. The width of the fluid-fluid interface is set to = 2 μm, which is sufficiently smaller than the elasto-capillary length G . That is, we could expect a soft-solid behavior and the formation of a three-phase contact angle, which should roughly follow Neumann's triangle. 37 Further parameters are Since a solid viscosity is not included in Reference 15, we have chosen a negligible value for s here. In the associated experiment, the data are measured examining a glycerol droplet in air on a silicone-gel layer. Let us note that ambient parameters, f ,0 , f ,0 , in Equation (16) are chosen to comply with Reference 15. There the authors deviated from the physical values of air as used in the experiment, 35 to improve numerical stability. The experimental results are only compared in the stationary state, which is independent on ambient viscosity and density.
As in Reference 15, we use an axisymmetric implementation of the method to incorporate the rotational symmetry of the experimental setup (see Appendix). At the bottom boundary of the solid substrate, we prescribe v = 0, while we set a free-slip condition at all other boundaries.
To resolve the dynamics of the wetting ridge formation, the time step size is set to 2 μs. The distance between the grid nodes is set to 0.39 μm along the interface, which leads to approximately 16 degrees of freedom across the interface region ( ∈ [0.1, 0.9]). Grid size is increased away from the interface region to a maximum node distance of 6.25 μm. The constant mobility is set to m = 0.01 m 3 s kg . We first consider the equilibrium profile of the substrate, for which the choice of fluid and ambient viscosities and the material densities is not relevant. Figure 3A shows the resulting profile at t = 6 ms in the (almost) stationary state. Maximum strain is about 30% in a very small region near the tip of the wetting ridge. A comparison with the reference simulation shows a very good agreement in both, the substrate profile under the drop as well as the wetting ridge. The slightly smaller height of the wetting ridge could be explained by the fact that we model a linear elastic material in contrast to the Neo-Hookean material in the reference. This is also confirmed by comparing to Reference 8 where a linear elastic model lead to a height of slightly less than 7 μm for the wetting ridge, which is exactly the same as in our simulations. However the results of Reference 8 are not included in the comparison as they deviate largely in the indentation of the soft substrate below the droplet. Already in the original paper the authors conjectured that this discrepancy is related to compressibility effects due to their pseudo-incompressible solid model. Moreover, the comparison of the simulation results with the experimental data from Reference 35 also shows a good agreement, especially in the wetting-ridge region. Since the stationary wetting ridge profile can be considered as a standard benchmark case for soft wetting, we published all data to reproduce Figure 3A in Reference 38. Figure 3B shows the evolution of the profile line between t = 100 μs and t = 500 μs, in which a large part of the clearly observable dynamic behaviour takes place. At this point, it is necessary to choose viscosities and densities that match the reference. Using the above mentioned values, we refer to test case 2 in Reference 15. A good agreement in the comparison to the reference is also shown here.

Time step stability
The used time step size of 2 μs was chosen so that in particular a precise analysis of the evolution of the wetting ridge is possible. The small time step size ensures that kink formation (∼ μs) and ridge formation (∼ ms) are reasonably well resolved. Many interesting soft wetting effects happen on larger time scales of seconds or minutes (e.g., droplet motion, cheerios effect, soft nucleation). Obviously, much larger time steps are necessry to resolve these phenomena within reasonable compute time. In this case, one can sacrifice temporal accuracy during the formation of the wetting ridge, as one is more interested to resolve the (slower) motion of an already formed ridge moving over the substrate. However, the necessary larger time steps may pose problems to numerical algorithms as stability restrictions kick in. Such restrictions were investigated for the above benchmark problem in detail in Reference 8. It was reported that The largest stable relaxation factor was found ∼10 −2 which means that on the order of 10 2 subiterations were necessary in each time step. This reduces the reported time step size of 10 μs to an "effective" time step size of 0.1 μs. Similarly, in Reference 14, time step sizes of 1 μs were chosen to conduct the same benchmark study. Nothing is reported on stability requirements there, but it can be expected that much larger time steps also posed problems, as the computations took already 47 hours on 128 cores. The unified approach for the fluid and solid equations presented here seems to stabilize the coupling significantly. The numerical simulations of the benchmark example are found to be stable for time step sizes of up to 16 μs without any subiterations. For higher time step size, the explicit coupling of velocity evolution and geometry evolution (grid movement) becomes unstable. While this time step size is already much larger than the sizes reported in previous literature on soft wetting (except for the fully monolithic approach in Reference 15), the stability can be further increased by considering matched viscosities and densities, which is a valid assumption to obtain the equilibrium state. For example, with the choice f ,0 = f ,1 = 1260 kg m 3 , f ,0 = f ,1 = s = 1.41 Pa ⋅ s, the time step size may be increased to 128 μs. Also, coarsening of the finite element grid leads to weaker time step restrictions. In our case, we found that doubling the mesh size enabled approximately 4 times larger time steps (512 μs). A more detailed study on time step restrictions, stability criteria, and time stepping errors is left for future work.
Let us note that the method still contains an explicit coupling of geometry evolution and flow computation. It can therefore not cope with a completely implicit coupling as very recently proposed in Reference 15. There the authors showed that their fully monolithic solution procedure can provide robust results at time step size of 1000 μs, and probably even above. The fully implicit coupling however complicates the solution procedure, including automatic differentiation, Newton iterations, and a technique called -continuation. Despite being not as robust as this fully implicit approach, our proposed method is more compact and can provide results at low computation time of each time step (35 s per time step for the benchmark example on the finest grid with ≈ 280,000 DOFs on a single core).

Fluid substrates
We now illustrate that our approach is capable to describe even much softer substrates. To this end, we consider the extreme case where the substrate is a purely viscous fluid without elasticity. This approach is particularly interesting for multi-phase flow with free surfaces, such as gas-water-oil mixtures. Figure 4 shows the scenario of the above benchmark test case, but now with G = 0 Pa. Figure 4B illustrates the significantly different dynamics during the first 1000 microseconds that result from the omission of an elastic influence. While the elastic substrate is already close to the equilibrium state at this time, the purely viscous substrate keeps deforming. The Laplace pressure pushes the droplet deeper and deeper into the substrate. As the substrate here is very thin, this dynamics would eventually lead to a breakup of the substrate phase in the center such that the drop touches the bottom boundary of the computational domain. Obviously, such a topological change is not correctly described by a grid-based numerical model. Here the simulations crash before the topological change due to degenerate mesh elements.
The observation of a fluid substrate also allows for a comparison of the equilibrium contact angle in the simulations with the Neumann's triangle. In Figure 5, we reproduce the results of three-phase flow simulations from fig. 10 of Reference 39 using the same relation of the surface tensions. We chose uniform values for viscosity and density throughout the whole domain here, namely, 1 Pa ⋅s and 1000 kg m 3 . The initial droplet is described as a small spherical cap upon a planar fluid substrate. The obtained equilibrium shapes are displayed in Figure 5. The intersection between the fluid-fluid meniscus and the substrate reveals a good agreement of the developed equilibrium contact angle. In particular, Figure 5A indicates that for = 1 = 0 the respective angles in all three phases are equal to each other (≈120 • ). Also, the case 1 ∶ 0 ∶ = 1 ∶ 1.4 ∶ 0.8 shown in Figure 5B yields the expected shape (cf. fig. 10 in Reference 39).

Interaction between liquid drops on viscoelastic substrates
After the validation studies for elastic and fluid substrates in the previous sections, we next illustrate the potential of the method. To this end, we present the first numerical simulations of two experimentally observed soft wetting phenomena in Sections 4.4 and 4.5.
In this section, we consider the substrate-mediated interaction of droplets. When placed sufficiently close, liquid drops on soft solid substrates are found to attract or repel each other, due to a combination of substrate elasticity and capillary forces. 40 This mechanism is also referred to as inverted Cheerios effect. With the present model, we are able to simulate such scenarios including the topological transition during fusion of drops.
We consider a 2D scenario with two nearby drops in the following. Figure 6 shows an exemplary series of snapshots to describe the dynamics in more detail. Here, we set a square domain Ω with length 150 μm, divided in half by Ω S and Ω f . We place two drops centered symmetrically so that they have an initial distance of 10.75 μm. The drops are initially prescribed as half-spherical caps of 20 μm radius. Furthermore we chose = 2 μm, which is significantly smaller than the distance between the drops, thus preventing the diffuse interfaces to "feel" each other. To strengthen the substrate-mediated interaction, we choose a relatively soft substrate, G = 200 Pa, with higher solid surface tensions, 1 = 0 = 46 mN m , and = 31 mN m . In the whole domain, we choose constant viscosity and density, namely, 1 Pa ⋅ s and 1000 kg m 3 , respectively. At the beginning there is a relatively rapid deformation of the substrate in vertical direction due to the capillary forces (t = 160 μs). This dynamics decays after t = 3200 μs and the attraction force between the droplets becomes dominant, moving the two drops in horizontal direction toward each other (t = 44,800 μs), see also eq. 3 of Reference 40. After the two drops have merged (t = 45,760 μs), the resulting single drop develops an almost stationary state, as indicated in the last snapshot.

Surfing on a viscoelastic substrate
Being able to describe viscoelastic substrates, the present model can simulate another interesting soft wetting phenomenon: the observed surfing behavior of liquid droplets. In Reference 31, it was reported that a fluid-fluid interface on a viscoelastic substrate moves in front of the wetting ridge, which suggested the term "surfing on the ridge." Different rheological models were considered analytically for this purpose, 31 and provide a potential reference for results from numerical simulations. Figure 7A shows an exemplary experimental setup from Reference 31. Water is injected into a hollow circular cylinder of viscoelastic material. While water is pushed in, ambient gas drains out and the fluid-fluid interface moves across the substrate, dragging the wetting ridge behind it.
To simulate this scenario, we prescribe a cylinder of 2 mm inner radius, 6 mm length, and 2 mm thickness. The parameters are motivated by the experiments in Reference 31. We set where v Flow denotes the average inflow velocity, which also determines the velocity of the fluid-ambient interface. Again, we use a higher viscosity of air to improve numerical robustness which barely perturbs numerical results. The injection of water is realized using the Dirichlet boundary condition of a parabolic velocity profile. Note that gravity is not included in the simulation setup as the forces at these small length scales are dominated by capillary effects. The evolution of the substrate profile around the contact line is shown in Figure 7B. One can see in the initial stages that the fluid-ambient interface moves faster in flow direction than the peak of the wetting ridge. This behavior is caused by the substrate viscosity which slows down substrate dynamics significantly. Accordingly, the three-phase contact point is displaced from the peak of the wetting ridge. This gives rise to the "surfing" of the contact line along the front of the ridge.
The dynamics of the solid domain is further illustrated by velocity streamlines shown in Figure 8 at two different times. Note that x-and y-axis are scaled equally here, which clarifies the actual flatness of the wetting ridge. It can be seen that two vortices appear in the early stage of surfing, complemented by a third vortex in a later stage of the surfing process.
The here developed numerical method is the first one able to simulate the surfing behavior and provides a framework to study viscoelasto-capillary dynamics beyond the limited theoretical study in Reference 31. Therefore, it can be used to increase the fundamental understanding of the dynamics of soft wetting with many potential applications, such as the patterning of cells 1 or droplets 2 onto soft surfaces, or the optimization of condensation processes, 3 which we leave for a further study.

CONCLUSIONS
Soft wetting is an emerging young field of research with many potential applications (e.g., References 1-3) and a richness of not fully understood physical phenomena (e.g., stick-slip motion, 31 durotaxis, 2 Shuttleworth effect 41 ) which demand for numerical simulation tools. We have contributed here a novel method which differs from previous methods in exploiting some characteristical properties of soft wetting. The result is a certain simplicity of the model (a unified approach for fluid and solid equations) and of the numerical discretization (no mappings, no subiterations, no nonlinear solvers needed). Yet, we have shown that the method is highly accurate and more robust than previous methods based on explicit coupling of the involved subproblems, only surpassed by the fully monolithic method presented in Reference 15. As special features, the method includes exact incompressibility and a linear monolithic assembly of the Navier-Stokes and Cahn-Hilliard equations which stabilizes dominant surface tension at small interface lengths. 18 Solid viscosity is included and we illustrated that solid material properties can be continuously tuned from purely viscous (i.e., three-phase flow) to purely elastic. With this we reproduced the surfing behavior of droplets on viscoelastic substrates for the first time. The robustness and flexibility of the method will permit detailed investigations of many soft wetting phenomena in the future, which we will address in a set of follow-up articles.

Rotational symmetry
A part of the results in Section 4 used axisymmetric computations, for which we now want to outline some details. Including rotational symmetry boils down to conducting simulations in the two-dimensional z-r plane (see Figure 2). The same finite element spaces are used as described in Equation (11), with dimension d = 2. Axisymmetry is realized on the level of differential operators by using an axisymmetric divergence and Laplace operator, respectively, which replace the corresponding operators in Equations (9)-(10). According to Figure 2, x Z and x R refer to the coordinates in lateral and radial direction, respectively. To compute the axisymmetric divergence of the viscous and elastic stress and the projection matrix, additional terms have to be implemented to account for the missing tensor components in 2D. For a more detailed overview of the time-discrete axisymmetric Navier-Stokes and Cahn-Hilliard system, see sec. A.2 of Reference 42. Further, the solid surface tension force is reformulated as ∇ Γ ⋅ ( s ( )P) = s ( ) n + ∇ Γ s ( ).
Here, is the total curvature of the solid surface which is computed in the discrete axisymmetric setting as in sec. 4.2 of Reference 33. Additionally, symmetric boundary conditions are imposed at the axis of rotation: n ⋅ = 0, n ⋅ = 0, n ⋅ v = 0.