Gradient flow formulation and second order numerical method for motion by mean curvature and contact line dynamics on rough surface

We study the dynamics of a droplet moving on an inclined rough surface in the absence of inertial and viscous stress effects. In this case, the dynamics of the droplet is a purely geometric motion in terms of the wetting domain and the capillary surface. Using a single graph representation, we interpret this geometric motion as a gradient flow on a Hilbert manifold. We propose unconditionally stable first/second order numerical schemes to simulate this geometric motion of the droplet, which is described using motion by mean curvature coupled with moving contact lines. The schemes are based on (i) explicit moving boundaries, which decouple the dynamic updates of the contact lines and the capillary surface, (ii) a semi-Lagrangian method on moving grids and (iii) a predictor-corrector method with a nonlinear elliptic solver upto second order accuracy. For the case of quasi-static dynamics with continuous spatial variable in the numerical schemes, we prove the stability and convergence of the first/second order numerical schemes. To demonstrate the accuracy and long-time validation of the proposed schemes, several challenging computational examples - including breathing droplets, droplets on inhomogeneous rough surfaces and quasi-static Kelvin pendant droplets - are constructed and compared with exact solutions to quasi-static dynamics obtained by desingularized differential-algebraic system of equations (DAEs).


Introduction
The dynamics and equilibrium of a droplet on a substrate are important problems with many practical applications such as coating, painting in industries and the adhesion of vesicles in biotechnology. The capillary effect, which contributes the leading behaviors of the geometric motion of a small droplet, is characterized by the surface tension on interfaces separating two different physical phases. Particularly, the capillary effect greatly affects the dynamics of the contact angle and the contact line of a droplet, where three phases (gas, liquid and solid) meet. Dating back to 1805, Young introduced the notion of mean curvature to study the contact angle of a capillary surface and proposed Young's equation for the contact angle, between the capillary surface and the solid substrate, of a static droplet. The geometric motion of a dynamic droplet is more challenging and extensively studied in the literature at the modeling [13,14,15,17,46], analysis [26,32,8,10,27,35,21,23], and computations level [58,48,45,49].
In this paper, we study the dynamics of a droplet placed on an inclined rough surface using vertical/horizontal graph representation. The contact angle hysteresis depends on the instantaneous contact angle of the droplet and the spatial heterogeneity of the substrate. Due to the roughness of the hydrophilic/hydrophobic substrate and moving contact lines, which leads to constant changes of local slope of the substrate, interesting phenomena such as pinning or capillary rise of the droplet will occur. After deriving the governing equations via gradient flows on a manifold, we propose first/second order unconditionally stable numerical schemes for the dynamics of the droplet and provide the first/second order convergence analysis for the quasi-static dynamics. Then we perform the accuracy check using quasi-static dynamics and conduct several challenging examples to accurately simulate the phenomena mentioned above.
Below, we give a more detailed outline of the results and strategies of the present paper. First, we give a kinematic description of a droplet using an incompressible potential flow, and derive the governing equations in a gradient flow formulation using a free energy including capillary effect and gravitational effect, and a Rayleigh dissipation function. After neglecting the inertial effect and the contribution of viscous stress inside the droplet, the dynamics of the droplet becomes a purely geometric motion, i.e., motion by mean curvature of the capillary surface coupled with motion of the contact lines; see Section 2. Using a single vertical/horizontal graph representation, we give a gradient flow formulation of the dynamic droplet by regarding the geometric motion of this droplet as a trajectory on a Hilbert manifold with boundary, where the obstacle occurs; see the resulting governing equations in Section 2.4 and derivations in Appendix A. Gradient flows on manifolds and the corresponding interpretation of minimizing movement with proper metrics have been the focus of recent research in both analytic and numerical aspects [5,4,41,1]. To completely describe the dynamics of the droplet, a free energy and a Riemannian metric (dissipation potential) in the gradient flow structure will be specific in different physical settings [13,46,3,17]. We also emphasize that there is a rich literature on droplets with different physical effects, such as viscosity stress inside the droplet, Marangoni effect, electromagnetic fields, electric fields or surfactant; see for instance [46,47,39,54,56,36,49,2].
The dynamics of the wetting domain for a 3D droplet is a challenging problem due to the geometric complexity. For example, cusp/corner singularity, topological changes such as merging and splitting. We refer to [32,22] for the studies of weak solutions to quasi-static contact line dynamics in the case that motion by mean curvature of the capillary surface is replaced by a Poisson equation. For the original mean curvature case, we also refer to [3,23] and the references therein for quantitative/qualitative stability theory of static droplets on rough surface. For the quasi-static dynamics where the capillary surface is determined by an elliptic equation, there are many analysis results on the global existence and homogenization problems; see [7,32,35,22] for capillary surface described by a harmonic equation and see [8,9,12,23] for capillary surface described by spatial-constant mean curvature equation.
Next, we propose unconditionally stable numerical schemes with second order accuracy for the 2D droplet dynamics described by the motion by mean curvature and the moving contact lines; see Section 3. The unconditionally stable schemes are based on an explicit boundary update, which decouples the computations for the dynamics of the contact lines and the capillary surface. In Section 3.1, we will give the stability and convergence analysis for the quasi-static dynamics based on the explicit boundary updates and some properties/dependence formulas of the quasi-static profile; see Proposition 3.1, Proposition 3.4 and Theorem 3.5, 3.6. For the full dynamic problem, the challenge of moving grids is handled by a semi-Lagrangian method with second order accuracy; see Section 3.3.3. To achieve a second order scheme efficiently, we also design a predictor-corrector scheme with a nonlinear elliptic solver upto second order accuracy; see Section 3.3.2.
Third, we construct several challenging and important computational examples to demonstrate the accuracy, validity and efficiency of our numerical schemes; see Section 4. (i) Using a quasistatic solution given by desingularized differential-algebraic system of equations (DAEs), we check the second order accuracy in space and time for our numerical schemes. (ii) We construct a breathing droplet with a closed formula solution, and we use it to show the long-time validity of our numerical schemes. (iii) We use complicated inclined rough surfaces such as the classical Utah teapot [43] to demonstrate the contact angle hysteresis of a droplet on inhomogeneous rough surface and the competition between capillary effect and gravitational effect. (iv) Using a horizontal graph representation and a desingularized formula for quasi-static dynamics, we also present simulations for Kelvin pendant drops with repeated bulges in the volume preserving case. For recent studies of steady solutions to Kelvin pendant droplet problem, we refer to [45]; see also [26,48]. Now we comment on existing references on different models and computations for droplets with hydrodynamic effects. Arguably, the contact line dynamics coupled with hydrodynamic effects is one of the mostly studied subjects in fluid dynamics. The contact lines experience an infinite driven force to overcome the no-slip boundary condition due to viscosity inside droplets. Various specific physical models describing slip boundary conditions and contact angle dynamics need to be imposed. For instance, the well-known Navier slip boundary conditions first used by [18] and the parameters in the Navier boundary condition can be determined by molecular dynamics [47]. A widely used moving contact line model coupled with fluids is derived in [47]; see (A.21) in Appendix A. In the coupled hydrodynamic model of droplets, the normal traction induced by the fluids is balanced with mean curvature F n = γ lg H, while in the purely geometric motion, the capillary surface is simply motion by mean curvature −ζv n = γ lg H with friction coefficient ζ. There are many numerical methods for the coupled hydrodynamic model of droplets. The level set method is first used in [39] and with reinitialization in [55] to simulate the moving contact lines. Various other related numerical methods are comprehensively reviewed in 2014 by [49], c.f. the front-tracking method [38,11], fixed domain method [42], the level set method [60, 16,50], the phase-field method [28,33,52] or the threshold dynamic method [20,19,51].
Instead, we focus on the purely geometric motion of droplets, in which the full dynamics is the motion by mean curvature of the capillary surface and the contact line dynamics; see Section 2. We mention particularly numerical methods that are closely related to the purely geometric motion case. The mean curvature flow with obstacles is theoretically studied in [4] in terms of weak solution constructed by a minimizing movement (implicit time-discretization). The threshold dynamics method based on characteristic functions are first used in [57,51] to simulate the contact line dynamics, which is particularly efficient and can be easily adapted to droplets with topological changes. They extended the original threshold method for mean curvature flows to the case with a solid substrate and a free energy with multi-phase surface tensions, in the form of obstacle problems. However, since they do not enforce the contact line mechanism [14,47], i.e., relation between contact line speed and unbalanced Young's force v cl = γ lg R (cos θ Y − cos θ cl ), thus their computations on contact angle dynamics are different with the present paper and only the equilibrium Young's angle θ Y is accurately recovered. Besides, the level set method developed in [39,55] can not be directly used and also can not deal with rough surfaces. Furthermore, to the best of our knowledge, the existing numerical methods for the contact line dynamics problem, including the level-set method, phase field method and threshold dynamics method, can not give the second order in time convergence as here proved in Theorem 3.6.
The rest of this paper is organized as follows. In Section 2, we describe the purely geometric motions of a dynamic/quasi-static dynamic droplet on inclined rough surfaces. With a specific free energy and a dissipation function, the derivation of governing equations using a gradient flow formulation is given in Appendix A. In Section 3, we propose the unconditionally stable 1st and 2nd order schemes for droplet dynamics on inclined rough surfaces. The stability and convergence analysis for the quasi-static equations are given in Section 3.1. The truncation error estimates and peusdo-codes for 1st/2nd order schemes are given in Appendix B and Appendix C respectively. In Section 4, we give some accuracy validations of our schemes compared with the quasi-static solution and demonstrate several challenging examples such as droplets in a teapot, breathing droplets and Kelvin pendant droplets.

Derivation of purely geometric dynamics of partially wetting droplets
In this section, we first revisit the kinematic equations of a droplet described by an incompressible potential flow and the dynamic mechanism driven by a free energy including capillary effect and gravitational effect, and a Rayleigh dissipation function; see Section 2.1, Section 2.2 and Section 2.3. After neglecting the inertial effect and the contribution of viscous stress inside the droplet, the dynamics of the droplet becomes a purely geometric motion, i.e., motion by mean curvature of the capillary surface and motion of the contact lines. In this case, using a vertical graph representation we describe the geometric motion of a droplet as a gradient flow on a Hilbert manifold; see the resulting governing equations for (quasi-static) dynamics of droplets with volume constraint in Section 2.4. Detailed derivations for the motion of a 3D droplet driven by a specific free energy and a Riemannian metric (dissipation function) are given in Appendix A. The governing equations of a 2D droplet on an inclined rough surfaces are presented in Section 2.5. When the vertical graph representation is broken, we replace it by a horizontal one; see Section 4.2.
2.1. Kinematic description of a droplet on a solid substrate. In this section, we first give a kinematic description of a droplet using a vertical graph function. More precisely, we study the motion of a 3D droplet placed on a substrate {(x, y, z); z = 0}. Assume the wetting domain is (x, y) ∈ D ⊂ R 2 with boundary Γ := ∂D. The droplet domain is identified by the area A := {(x, y, z); (x, y) ∈ D, 0 < z < u(x, y), u| ∂D = 0}.
Assume the fluid inside the droplet follows an incompressible potential flow with velocity v(x, y, x) = ∇φ(x, y, z) and constant density ρ. The motion of the shape of this droplet is described by the following two kinematic boundary conditions. The motion of the capillary surface on ∂A ∩ {z > 0} with the outer normal n is described by the normal speed and the motion of Γ (physically known as contact lines) with outer normal n cl in x-y plane is described by the contact line speed Notice the incompressible potential flow satisfies ∇ · v = 0. Hence together with kinematic boundary condition, we have For the dynamic wetting domain and droplets, all the quantities A(t), D(t), u(x, y, t), Γ(t) depend on time variable t. The compatibility condition for (2.3) is ∂A v n ds = 0, which turns out to be equivalent to the volume preserving constraint. Indeed, using the normal speed v n = ∂tu √ 1+|∇u| 2 in the graph representation, we have Then by u(x, y, t) = 0 on Γ(t) and the Reynolds transport theorem, we have where the last equality follows from the volume preserving constraint d dt D(t) u(x, y, t) dx dy = 0. Hence in the volume preserving case, the motion of the droplet can be completely described by the motion of capillary surface u(x, y, t) and the motion of contact domain D(t).

2.2.
Free energy for the droplet and Young's angle. To give a specific free energy, we will follow the same notations and terminologies as in the classical book of de Gennes [15]. To consider the interactions between the three phases of materials: gas, liquid, and solid, denote γ sl (γ sg , γ lg resp.) as the interfacial surface energy density between solid-liquid phases (solid-gas, liquid-gas resp.). γ sl , γ sg , γ lg > 0 are also known as the surface tension coefficients. Surface tension plays the most important role in the dynamics and equilibrium of the droplet. Surface energy consists of the contributions from the capillary surface (with surface tension γ lg ) and the wetting domain of the droplet (with the relative surface tension γ sl − γ sg ). Surface energy between liquid and gas is the excess energy due to the one half lower coordination number (in the mean field approximation) of molecules at the surface compared with those sitting in the liquid bulk (Doi [17]). We take the total free energy of the droplet as the summation of surface energy, gravitational energy and kinetic energy where ρ is the density of the liquid, g is the gravitational acceleration and φ satisfies (2.3). For droplets with small Weber number, the contribution of the fluid's inertia is small compared to the capillary effect. Thus for small droplets, we neglect the inertia effect and drop the last term in the free energy F has units of energy. We neglect other effects, such as viscosity stress inside the droplet, Marangoni effect (solutocapillary and thermocapillary effect), electromagnetic fields, evaporation and condensation, etc. The free energy (2.7) for small droplets in the current setup is particular useful for practical applications such as coating, painting in industries and the adhesion of vesicles in biotechnology.
The competition between the three surface tension coefficients will determine uniquely the steady shape of the droplet with a fixed volume V . Let the density of gas outside the droplet be ρ 0 = 0. We denote the capillary coefficient as ς := ρg γ lg and the capillary length as L c := 1 √ ς . For a droplet with volume V , its equivalent length (characteristic length) L is defined as V = 4π 3 L 3 in 3D and V = πL 2 in 2D. The Bond number Bo := ( L Lc ) 2 = ςL 2 shall be small enough to observe the capillary effect [15]. Notice for simplicity in presentations of the governing equations, we allow ς < 0 in the case of pendant droplet. Hence when ς < 0, the capillary length is L c = 1 √ |ς| and the Bond number is Bo = ( L Lc ) 2 = |ς|L 2 . Define σ as the relative adhesion coefficient between the liquid and the solid Adhesive forces between the liquid and the solid cause the liquid drop to spread across the surface (called a partially wetting liquid on a hydrophilic surface), while cohesive forces within the liquid cause the drop to ball up and avoid contact with the surface (called dewetting or non-wetting liquid on a hydrophobic surface). We remark that the spreading parameter S := γ lg γsg−γ sl γ lg − 1 could be positive in the so-called total wetting regime [14, Section 1.2.1]. In this case, the liquid spreads completely to a film of nanoscopic height, which can not be described using contact angle dynamics in the current paper. For the present contact angle dynamics setup, |σ| ≤ 1. By Young's equation [59], the equilibrium contact angle θ Y is determined by the Young's angle condition We call it the partially wetting liquid case (hydrophilic surface) if −1 < σ = − cos θ Y ≤ 0, 0 < θ Y ≤ π 2 , while we call it the non-wetting liquid case (hydrophobic surface) if 0 ≤ σ = − cos θ Y < 1, π 2 ≤ θ Y < π. The case θ Y = 0, σ = −1 is called completely wetting.

2.3.
Friction force for the motion of droplet and Rayleigh dissipation function. There are three types of friction and viscosity force on the droplet. First, the contact line friction force density is given by −Rv cl n cl = −R(n cl · ∇ x,y φ)n cl , where in 3D, R is the friction coefficient per unit length for the contact line with the units of mass/(length · time). Second, the friction force density on the capillary surface is given by −ζv n n, where in 3D, ζ is the coefficient per unit area for the capillary surface with the units of mass/(area · time). Third, the viscosity stress inside the droplet is µ(∇v + ∇v ⊥ ) where µ is the dynamic viscosity. We neglect the viscosity effect inside the droplet. Then the Rayleigh dissipation function (with the unit of work) is given by [31] (2.9) After neglecting the kinetic energy and viscosity dissipation inside the droplet, the dynamics of the droplet is purely a geometric motion driven by the free energy (2.7) and Rayleigh's dissipation function (2.9). Therefore, the motion of the droplet can be completely described by the geometric configurations: the boundary of wetting domain Γ(t) and capillary surface u(x, y, t), instead of by velocity potential φ.

2.4.
Dynamics of a droplet derived by gradient flow on manifold. In Appendix A, we will derive the gradient flow of F(η), η(t) = (Γ(t), u(t)) on a Hilbert manifold with respect to a Riemannian metric g η suggested by (2.9). For a 3D single droplet, this yields the following governing equations with initial data (Γ(0), u(x, y, 0)) and initial volume V . When there are topological changes, including merging and splitting of droplets, (2.10) becomes parabolic variational inequality; see Appendix A for more discussions.
In the following proposition, we summarize the dissipation relation, the quasi-static dynamics and the contact line speed mechanism. The proof of Proposition 2.1 is given in Appendix A.1 for completeness.
Proposition 2.1. Assume (Γ(t), u(x, y, t), λ(t)), t ∈ [0, T ] are solutions to (2.10) with regularities u ∈ H 2 0 (D(t)) and Γ(t) is a C 1 Jordan curve with a finite perimeter. Then we have (i) the velocity relation on the contact line (ii) the energy-dissipation relation (iii) if ζ = 0, the resulting quasi-static dynamics is a gradient flow for Γ(t) For the cases that singularity occurs on Γ(t), the solution (Γ, u) shall be understood in the viscosity sense with some geometric assumptions on the wetting domain D(t), which is beyond the scope of this paper. We will use the statement (iii) above, together with a DAEs solver, to check the accuracy of our numerical schemes proposed in next section.
The dimensional analysis for a 2D droplet is given in Appendix A. The resulting dimensionless equations in 2D are Here, with typical length scale L and time scale T , the coefficients κ = L 2 ς,λ = L γ lg λ, β = ζL 2 γ lg T and typical 2D volume [15] V = π are all dimensionless. The capillary number for the capillary surface β measures the ratio between the frictional force on capillary surface and surface tension, and indicates the capillary relaxation time on the capillary surface. The Bond number κ measures the ratio between the gravitational force and surface tension. In the remaining of this paper, we will use the dimensionless formulation (2.15) after dropping hat inλ.

2.5.
Governing equations for a single 2D droplet on a rough and inclined surface. In this section, with some modifications for the free energy, the Riemannian metric and the same derivation of the gradient flow formulation as in Appendix A, we summarize the governing equations for a single droplet on a rough and inclined surface. Given a rough solid surface, we follow the convention for studying droplets on an inclined surface and choose the Cartesian coordinate system built on an inclined plane with effective inclined angle θ 0 such that − π 2 < θ 0 < π 2 , i.e., (tan θ 0 )x is the new x-axis we choose; see Fig 1 (b). With this Cartesian coordinate system, the rough surface is described by the graph of a function w(x) and the droplet is then described by The motion of this droplet is described by the relative height function (capillary surface) u(x, t) ≥ 0 and partially wetting domain a(t) ≤ x ≤ b(t) with free boundaries a(t), b(t). Consider a manifold Given any q 1 = (α 1 , β 2 , v 1 ), q 2 = (α 2 , β 2 , v 2 ) ∈ T η M, define the Riemannian metric g η : T η M × T η M → R as The dynamics of a droplet on rough surface can be regarded as a trajectory η(t) = (a(t), b(t), u(x, t)) on M. The tangent direction of the curve η(t) is given by η (t) := (a (t), b (t), ∂ t u) ∈ T η(t) M. Now we consider the energy functional F : M → R associated with the rough surface (2.20) where ρ is the density of the liquid, g is the gravitational acceleration. In the inclined case, for a droplet with volume V in 2D, the effective Bond number is (2.21) Bo := L L c 2 cos θ 0 = ςL 2 cos θ 0 .
Then by same derivations as the gradient flow formulation in Appendix A, with h(x, t) := u(x, t) + w(x), the governing equations in the dimensionless form for a single droplet become where the two angles are defined as ∂ x w| a = tan θ 0a , ∂ x h| a = tan(θ 0a + θ a ) and ∂ x w| b = − tan θ 0b and ∂ x h| b = − tan(θ 0b + θ b ); see Fig 1 (b). Recall β, κ, V, λ are all dimensionless. It is easy to check the steady state a (t) = b (t) = 0 recovers Young's angle condition.
Remark 1. For w(x) = 0, i.e., the surface is a perfect inclined plane with angle θ 0 , the derivation above recovers the model for capillary droplets on an inclined surface [25,8,35].
is invariant with respect to the translation x 0 . As a consequence, the dynamics and the equilibrium profile do not depend on the coordinates we choose. More importantly, at the equilibrium, the right hand side in the first equation is exactly the hydrostatic balance [37, Section 61] due to the force balance on the capillary surface. With w ≡ 0, [25,Theorem 8.3] proved the nonexistence of the static profile for θ 0 = 0, π. With a sufficient substrate roughness, [8] proved the existence of the static profile.
Remark 3. In the case without volume constraint, we can calculate the rate of change of the volume The shrink estimate of the 2D droplet suggests the volume preserving constraint is necessary to observe interesting phenomena for long enough time. For this reason, all the numerical examples are with volume constraint.

numerical schemes for droplets dynamics on a rough and inclined surface
In this section, we consider a droplet (described by a vertical graph function) on a rough and inclined surface in the partially wetting case, i.e., the relative adhesion coefficient −1 < σ ≤ 0. Since a uniform estimate for the moving boundaries a(t), b(t) in (2.22) can be obtained, we have an unconditionally stable explicit scheme for the time stepping of the moving boundary, which in turn leads to the convergence of the numerical scheme. This explicit discretization decouples the computations for the moving boundary a(t), b(t) and capillary surface u(x, t). In Section 3.1, to first illustrate the idea, we give the stability and convergence analysis for the first/second order schemes for the quasi-static dynamics; see Proposition 3.1 and Theorems 3.5, 3.6. To design the numerical schemes for the full dynamics of droplets and achieve a second order scheme in time and space, we should particularly take care of the following issues. First, due to the moving grids at each time step, we need to map the capillary surface at different time to the same domain based on a semi-Lagrangian method upto second order accuracy. Second, to achieve a second order scheme efficiently, we design a predictor-corrector scheme with a nonlinear elliptic solver, which maintains the overall second order accuracy. Third, the effects from the inclined rough substrate and the volume constraint will be included. We will derive the first order scheme and give its truncation error in Section 3.2. Then we derive the second order scheme and give its truncation error in Section 3.3. The proofs for truncation error estimates will be left to Appendix B. Before we present the schemes, we list some key notations below in Table 1. 3.1. Stability analysis and convergence of numerical schemes for quasi-static dynamics. In this section, to illustrate the idea, we will first introduce the first order/second order scheme for w ≡ 0, θ 0 = 0 and the quasi-static case, i.e., β = 0 in (2.15) Then we give the stability analysis and convergence result in Proposition 3.1 and Theorems 3.5, 3.6 respectively. Based on the observation for the unconditional stability and convergence, in Section Table 1. Commonly used notations in this paper.

Symbols Meaning
Exact moving boundaries at t n a n , b n Numerical moving boundary at t ñ a n+1 ,b n+1 Predictor numerical moving boundary at t n+1 Numerical spatial variable at t n x n j = a n + jτ n , τ n = b n −a n N , j = 0, · · · , N Moving spatial grids at t ñ Predictor moving grids at t n+1 PDE solution at t n h n j , j = 0, · · · , N Numerical solution at time t n and spatial grid x n j h n (x n ) for x n ∈ [a n , b n ] Numerical solution at t n (with continuous spatial variable) Predictor numerical solution at t n+1 h n+1 j Predictor numerical solution at t n+1 and gridx n+1 Intermediate numerical rescaled solution from predictor 3.2 and Section 3.3, we will design the first/second order numerical schemes for the full dynamic problem, i.e., β > 0. We first present the first/second order schemes and then prove the stability and convergence. Let t n = n∆t, n = 0, 1, · · · with time step ∆t. Approximate a(t n ), b(t n ), u(t n ) by a n , b n , u n respectively. First order scheme: Step 1 . Explicit boundary updates. Compute the one-side approximated derivative of u n at b n and a n , denoted as (∂ x u n ) N and (∂ x u n ) 0 . Then by the dynamic boundary condition in (2.22), we update a n+1 , b n+1 using Step 2. Update u n+1 and λ n+1 implicitly.
where the independent variable for u n+1 is x n+1 ∈ (a n+1 , b n+1 ). Second order scheme: Step 1 . Repeat the Step 1 and Step 2 for the first order scheme. Denote the results as the predictor a n+1 ,b n+1 ,ũ n+1 .
Step 2. Explicit boundary updates based on a predictor-corrector method. Compute the oneside approximated derivative of u n at b n and a n , denoted as (∂ x u n ) N and (∂ x u n ) 0 . Then update a n+1 , b n+1 using Step 3. Update u n+1 and λ n+1 implicitly according to (3.3).
. Then for n < T ∆t , we have (i) the estimate for endpoints (iii) the estimate for u and u x where C is a constant depending only on the initial data a 0 , b 0 and T .
Proof. First, from (3.2) and (3.4), we know for both first and second order schemes, This leads to statement (i). Second, multiplying the first equation in (3.3) by u n+1 and integration by parts show immediately that λ n+1 > 0 since κ > 0. Then integrating the first equation in (3.3) from a n+1 to b n+1 , we have Then by (i) we have Third, multiplying the first equation in (3.3) by u n+1 and integration by parts show that This, together with the estimate for λ in (3.6) and (i), gives the estimate for ∂ x u and u and we conclude (3.7).
Before proving the convergence of the scheme, we first clarify the existence and properties in Proposition 3.4 for the quasi-static solution to It is easy to see (3.13) is translation invariant for x → x + , so without loss of generality we assume . Let θ be the contact angle such that tan θ = −∂ x u| b and thus (3.14) cos θ = 1 whose derivation via gradient flow is given in Appendix D for completeness. Equation (3.15) can be used to describe not only the quasi-static profiles with single vertical graph representation but also profiles with horizontal graph representation; see more details in Kelvin pendant drop problem (4.26). For instance, for the simple case κ = 0, the quasi-static profile is given by the spherical cap formula (4.11) with 2D volume formula (4.13). If 2V < πb 2 then the spherical cap has a single vertical graph representation. Given a contact angle θ, the existence and symmetry of a static droplet has been comprehensively studied in [25,24]. However, our numerical schemes require the existence of (3.13) for a given contact boundary b and crucially relies on the continuous dependence on b. Propostion 3.4 below obtains the unique critical wetting domain b c corresponding to θ c = π 2 such that the quasi-static solution to (3.13) has a single vertical graph representation. Moreover, it gives the estimate for the continuous dependence of the contact angle θ with respect to contact boundary b. This is also the key for the convergence analysis later. Before stating Proposition 3.4, we first give the following lemma for the relations between b, λ, θ, u m . Lemma 3.3. Suppose κ > 0 and the volume of the droplet is V . Then for any contact angle 0 < θ ≤ π 2 , the droplet profile obtained in Lemma 3.2 satisfies the following relations among θ, u m , b (3.16) Proof.
Step 1. Integrating once in the first equation of (3.15), we have Step 2. From (3.17), we know since κ > 0, the graph of J(u, θ) for fixed θ is a parabola open downwards. Hence its axis of symmetry is located to the right of u m and thus J(u, θ) is increasing w.r.t u for 0 ≤ u ≤ u m , which implies In particular, 0 < κu m ≤ λ and Now we derive the relations between λ and u m . Since J ≥ 0 Since X(u m ) = 0, we have the integral formula Then from the volume constraint, Combining (3.20), (3.24) and (3.25), we conclude Step 3. To further eliminate λ, we introduce the variable v ∈ [0, 1] such that u = u m v. Then by changing variables u = u m v in (3.18), we have where we used the first equation in (3.26) to eliminate λ. Thus the relations between θ, u m , b in (3.26) can be simplified as Proposition 3.4. Suppose κ > 0 and the volume of the droplet is V . For any 0 < θ ≤ π 2 , assume θ, u m , b satisfy (3.16). Then (i) the function u m (θ) obtained by Lemma 3.2 is strictly increasing w.r.t θ from u m (0 + ) = 0 to u m ( π 2 ) = u c m ; as a consequence, the inverse function θ = θ(u m ) maps u m ∈ (0, u c m ) onto θ ∈ (0, π 2 ) such that F 1 (u m , θ(u m )) = 0 and ∂(cos θ(um)) ∂um < 0; 2 ); and we have the estimate Proof.
Step 1. We first show statement (i) via inverse function theorem.
From the first equation in (3.16), taking partial derivative of F 1 with respect to cos θ, we have (3.20). Moreover, taking partial derivatives of F 1 w.r.t u m , we have From the inverse function theorem we know there exists a unique function θ = θ(u m ) such that From the definition of F 1 (u m , θ), it is easy to verify u m = 0 when θ = 0. Denote the value of u m corresponding to θ = π 2 as u c m . We conclude statement (i).
Step 2. Combining statement (i) and the second relation in (3.16), it is easy to see Step 3. We now show u m = u m (b) is a function of b via the inverse function theorem. Taking partial derivatives of F 2 w.r.t u m , we have From statement (i), we plug ∂ cos θ ∂um into (3.33) to see From Hölder's inequality, we have which shows I 2 in (3.34) is nonnegative. On the other hand, from (3.21), we have so the quasi-static profile X(u) is concave. Thanks to the concavity, for any α ∈ (0, 1), where we used J(v, θ) ≤ 1 similar to (3.30). Combining (3.35) and (3.36), we obtain Therefore b(u m ) is a strictly decreasing function w.r.t. u m . By the inverse function theorem, we conclude u m = u m (b) is a function of b. Moreover, from the definition of F 2 , it is easy to see when u m = 0, θ = 0, we have b = +∞. Denote the value of b corresponding to u c m as b c . We conclude (ii).
Step 4. From statement (i) and (ii), the composition θ(b) = θ(u m (b)) is a function of b and we conclude (iii).
Step 5. Finally, we give the estimate in statement (iv). Combining (3.32) and (3.37), we have Proposition 3.4 gives the continuous dependence of the contact angle θ = θ(b) with respect to the contact point b for symmetric contact points (a = −b), so we conclude the existence of solutions to (3.1) by the well-posedness of the ODE system where cos θ(·) is a function of b−a 2 . Now we state and prove the convergence result for the first order scheme (3.2)-(3.3).
are the exact solution to (3.1) and let a n , b n , u n (x n ), x n ∈ [a n , b n ] at t = t n be the numerical solution obtained from the first order scheme (3.2)-(3.3) with the same initial data (a 0 , b 0 , u 0 ). Then for n < T ∆t , we have the convergence where C m is the bound in (3.29).
Proof. Without loss of generality, we assume initially −a = b > 0 which is dynamically preserved for both exact solution and numerical scheme. So we have −a(t) = b(t), −a n = b n and we only prove the convergence for b. First, from the Taylor expansion for the exact solution and the boundary condition in (3.1), we have From the the estimate 0 ≤ ∂(cos θ(b)) ∂b ≤ C m in (3.29) and the boundary condition in (3.1), we have Now, subtract (3.42) from the boundary update (3.2) and denote ε n := |b(t n ) − b n |. From (3.29) and (3.43), we have ∂b ≤ C m since given b n the numerical profile u n (thus cos θ(b n )) solved in the first order scheme (3.3) is also quasi-static profile and (3.29) in Proposition 3.4 holds.
Now we state the convergence result for the second order scheme (3.4)-(3.3) and omit the proof.
are the exact solution to (3.1) and let a n , b n , u n (x n ), x n ∈ [a n , b n ] at t = t n be the numerical solution obtained from the second order scheme (3.4)-(3.3) with the same initial data (a 0 , b 0 , u 0 ). Then for n < T ∆t , we have the convergence where C depends only on the bound C m in (3.29).
3.2. First order unconditionally stable scheme based on explicit boundary update and semi-implicit motion by mean curvature. Based on the observation for the unconditional stability and convergence for the quasi-static dynamics of the droplet (in Section 3.1), in this section, we design the first order scheme and give the truncation error estimate.
3.2.1. First order scheme based on explicit boundary update and semi-implicit motion by mean curvature. Now we design a numerical algorithm for the motion by mean curvature case (2.22). With some proper spatial discretizations (such as finite difference, finite element, spectral approximation), we approximate h(x n , t n ) by h n (x n ) for x n ∈ (a n , b n ) and approximate λ(t n ) by λ n . We propose the following three-step algorithm for updating a n , b n , h n , λ n from time step t n to t n+1 .
Step 1. Compute the one-side approximated derivatives of h n at b n and a n , denoted as (∂ x h n ) N and (∂ x h n ) 0 . Then by the dynamic boundary condition in (2.22), we update a n+1 , b n+1 using Step 2. Rescale h n from [a n , b n ] to [a n+1 , b n+1 ] with O(∆t 2 ) accuracy using a semi-Lagrangian discretization. For x n+1 ∈ [a n+1 , b n+1 ], denote the map from moving grids at t n+1 to t n as (3.50) x Define the rescaled solution for h n as It is easy to verify by the Taylor expansion h n * ( Step 3. Update u n+1 and λ n+1 semi-implicitly. where the independent variable is x n+1 ∈ (a n+1 , b n+1 ) and V is the initial volume. For convenience, we provide a pseudo-code for this scheme in Appendix C.1. Similar to (3.5) in Proposition 3.1, from (3.49), we know for n∆t < T , since σ < 0 (3.53) so the explicit scheme for the moving boundaries is unconditionally stable.

3.2.2.
Truncation analysis for the first order scheme. Here, we state the truncation error for the first order scheme.
Lemma 3.7. Assume a(t), b(t) ∈ C 2 ([t n , t n+1 ]) and h(x, t) ∈ C 4,2 x,t ([a(t), b(t)] × [t n , t n+1 ]) be the exact solution to (2.22) with initial data at t = t n , a n , b n , h n (x n ) for x n ∈ [a n , b n ]. Then we have the first order truncation error estimates where h n * is given by For simplicity, h(t n+1 ) represents h(·, t n+1 ) in the lemma above and the remaining contents. By mapping moving domain to a fixed domain , the proof of this lemma is standard so we put it in Appendix B.
3.3. Second order numeric algorithm based on a predictor-corrector method with an unconditionally stable explicit boundary update. In this section, we use the predictor-corrector method to obtain a second order scheme. With the notations in Table 1, we still approximate a(t n ), b(t n ) by a n , b n respectively and approximate h(x n , t n ) by h n (x n ) for x n ∈ (a n , b n ). However, it is more convenient to use a fixed domain variable which is equivalent to Denote U (Z, t) := h(x, t). We will first present the second order numeric algorithm in Section 3.3.1 and then we give the derivation of the second order scheme in Section 3.3.2 and Section 3.3.3 based on the relation Here t n+ 1 2 := (n + 1 2 )∆t and Z(x n+ 1 with the independent variable x n+ 1 2 .
3.3.1. Second order predictor-corrector scheme and unconditional stability for explicit boundary update. Now we present the second order scheme with continuous spatial variables.
Step 1. Predictor. Since we show in Section 3.3.2 that the nonlinear elliptic solver for motion by mean curvature requires second order accuracy, after updating a n+1 , b n+1 by the first order scheme in Section 3.2, we replace the semi-implicit elliptic solver by an implicit nonlinear elliptic solver.
Precisely, for the independent variable x n+1 ∈ (a n+1 , b n+1 ), (3.59) β where h n * (x n+1 ) is the first order intermediate profile given in (3.51) and V is the initial volume.
Denote the results as the predictorã To solve (3.59), one can use standard Newton's iterative method.
Step 2. Explicit boundary update. Compute the one-side approximated derivative of h n at b n and a n , denoted as (∂ x h n ) N and (∂ x h n ) 0 . Then update (3.60) a n+1 − a n ∆t = 1 2 Step . We will give detailed derivation for the choice of the second order intermediate solutionh n * (x n+1 ) in Section 3.3.3. For convenience, we provide a pseudo-code for this scheme in Appendix C.2.

3.3.2.
Derivation of a second order scheme based on the predictor-corrector method for DAEs with an algebraic solver upto second order accuracy. To design a second order scheme, we illustrate the idea using the predictor-corrector method for an analogous DAEs. Assume we have an exact DAEs where the second algebraic equation is equivalent to u = G(b) for some function G. However, in practice, one may not solve u = G(b) exactly, which in our case, is to solve a nonlinear elliptic equation (3.59). Therefore, we design a predictor-corrector method to solve a DAEs with an algebraic solver upto second order accuracy. Let b n , u n be given such that u n = G(b n ) + O(∆t 2 ).
Step 1. Solve the predictorb n+1 by forward Euler scheme Step 2. Obtain the predictorũ n+1 by solving algebraic equation up to a second order accuracy Step 3. Solve the corrector b n+1 by the trapezoidal method Step 4. Obtain the corrector u n+1 by solving the algebraic equation up to a second order accuracy Indeed, we show the second order error estimate of this scheme below. Denote function which gives the second order accuracy for b n+1 and thus u n+1 .

3.3.3.
Derivation of the second order accuracy for the semi-Lagrangian termh n * (x n+1 ). Now we derive the second order scheme based on (3.58). Notice the spatial grids are moving along time. We need to map grids at different time back to the same fixed domain Z ∈ [0, 1] based on (3.58). Furthermore, to achieve second order accuracy, we apply midpoint scheme and define We illustrate the second order accuracy for the following term involving time derivative, for inde- Below, we approximate I 1 , I 2 , I 3 upto second order accuracy.
First, notice Z(x, t) at different time is related by (3.58). Thus whenever we evaluate some quantity U at different time, for instance at t n+1 , we mean U (Z(x n+1 , t n+1 ), t n+1 ). Recall U (Z, t) = h(x, t). Therefore, by midpoint scheme, I 1 (x n+ 1 2 , t n+ 1 2 ) = ∂ t U (Z, t n+ 1 2 ) can be approximated by where we use the numerical solution h n+1 (x n+1 ) ≈ h(x n+1 , t n+1 ) = U (Z, t n+1 ) and h n (x n ) is similar. Here the equality holds in the sense of changing variables x n , x n+ 1 2 , x n+1 to the fixed domain Z = x−a(t) b(t)−a(t) ∈ [0, 1] and Z = Z(x n+1 , t n+1 ) = Z(x n , t n ) = Z(x n+ 1 2 , t n+ 1 2 ). Second, by midpoint scheme, I 2 = ∂ Z U (Z, t n+ 1 2 ) can be approximated by Recall the scheme in Section 3.3.1 useh n+1 (x n+1 ) forx n+1 ∈ [ã n+1 ,b n+1 ] as predictor instead of the nonlinear unknown h n+1 (x n+1 ). From (B.13), we have the |x n+1 − x n+1 | = O(∆t 2 ), which enables us to replace the unknown term by the predictor. Then changing the intermediate variable Z back to x gives Third, we turn to approximate I 3 . Still by the midpoint scheme, the last term I 3 = ∂ t Z(x n+ 1 2 , t n+ 1 2 ) in (3.69) can be approximated by Notice from (3.57) and (3.68), we have which is recast as Notice also the relation a n+1 − a n+ 1 2 = a n+ 1 2 − a n = a n+1 −a n 2 . Therefore, the last term I 3 = ∂ t Z(x n+ 1 2 , t n+ 1 2 ) in (3.69) can be approximated by where we used (3.72) in the third equality. Therefore, we now define the intermediate variableh n * (x n+1 ) such that Using the approximated formulas for I 1 , I 2 , I 3 above, we propose the semi-Lagrangian term (3.74) In summary, we have the second order approximation Since this is a key step in the numerical discretization, so we also give the second order spatial discretization ofh n * (x n+1 ) to see it has a similar form with (3.51). Denote spatial grid size τ n = b n −a n N and τ n+1 = b n+1 −a n+1 N . Notice the second order spatial discretizations for I 2 , I 3 are Define the second order spatial discretization (3.76)

3.3.4.
Second order truncation error estimates for the predictor-corrector method. The strategy of the second order truncation error estimates is the same as that of Lemma 3.7. Namely, we notice that in a fixed domain in terms of Z ∈ [0, 1] the predictor-corrector method gives us a second order scheme and then we prove the mapping from Z to x n+1 keeps the second order accuracy. For completeness, we put the proof of Lemma 3.8 in Appendix B.

Validations and computations
In this section, we will first use the DAEs solution for the quasi-static dynamics to check the second order accuracy of the numerical schemes proposed in the last section. Then we design several challenging and important examples: (i) a periodic breathing droplet with a closed formula solution and a long-time computational validation; (ii) dynamics of a droplet on an inclined rough surface and in a "Utah teapot"; (iii) Kelvin pendant droplet with repeated bulges and the corresponding desingularized DAEs solver for quasi-static dynamics based on horizontal graph representation.

Desingularized DAEs formula and accuracy check for the quasi-static dynamics.
In this section, we will first derive the DAEs for the quasi-static dynamics using a desingularized formula. Then we give an accuracy check for the case w(x) = 0 and θ 0 = 0 using the corresponding quasi-static solutions, which can be obtained using the desingularized DAEs solver.
with boundary condition u(a(t), t) = u(b(t), t) = 0. Multiplying (4.1) by u and integration by parts give immediately that κ > 0 implies λ > 0. In this subsection, we will derive the following DAEs for b(t), λ(t), u m (t) in three steps, which completely describes the quasi-static motion where u m (t) is the maximal point of u(x, t) . Here J(u) is the shorthand notation for the function J(u, θ) = − κ(u 2 −u 2 m ) 2 + λ(u − u m ) + 1, where cos θ is solved from (3.19). Recall θ is the contact angle such that tan θ = −∂ x u| b . Thus we have (3.14) and from (3.19) On the other hand, to desingularize X u in the numerical implementation, denote ψ := √ u m − u. Then we have Recall (3.21) and particularly in the case u xx (0) < 0, we have Thus one can check can be approximated by .
With the desingularized formula (4.6), there is no singularity in the DAEs (4.2) for b(t), λ(t), u m (t), so it can be solved efficiently and accurately by any ODE solver such as ode15s in Matlab, whose results will be used to check the accuracy for our PDE solvers. Furthermore, we can solve the capillary surface X(u, t) by the integral formula (3.24). Finally, we give the equilibrium solution for DAEs (4.2). Taking b (t) = 0 in (4.2), we have the algebraic equation and can solve uniquely the steady solution (b, u m , λ).

4.1.2.
Accuracy check between DAEs and 1st/2nd order PDE solvers. We first use the DAEs solver ode15s in Matlab to solve DAEs (4.2) with the initial data Compared with the DAEs solution, we show below the accuracy check for the first order scheme in Section 3.2.1 and the second order scheme in Section 3.3 in Table 2. The absolute error in ode15s is set to be 10 −13 , which is smaller than the absolute error in the accuracy check Table 2. The residual tolerance in Newton's iteration in the second order scheme is set to be 10 −12 . We use the same parameters β = 0, κ = 0.1, initial angle θ in = 1.   in the graph representation. Here u is the piecewise graph representation of capillary surface. When κ = 0, the governing equation for the quasi-static dynamics becomes H = λ, where λ is a function of t. This means the quasi-static profile has constant mean curvature λ(t) everywhere on the capillary surface. Assume the initial droplet has the wetting domain {x ∈ R d−1 ; |x| ≤ b(0)}. Due to the rotation invariance for both equation and initial wetting domain, the solution will remain axially symmetric, denoted as u(r, t). As a consequence, the quasi-static profile is a spherical cap, whose center may be above the ground. To describe this spherical cap solution, we denote the height of the center as u * (t) ∈ R. Furthermore, notice the mean curvature of a d-dimensional ball is H = d−1 R , where R is the radius of the spherical cap.
Consider the partially wetting case in 3D when the droplet is represented by the single graph we can solve Then u * (t) = 2 λ(t) − u m (t) and we have (4.11) (u(r, t) − u * (t)) 2 + r 2 = R(t) 2 .
For a droplet in the non-wetting case, the capillary surface can not be expressed uniquely by the graph of a function u(r). In the non-wetting case, in which the center u * (t) is above the ground, one shall use two graph representation (with same notations) for 0 ≤ u ≤ u * (t) and u * (t) ≤ u ≤ u m respectively; see also the horizontal graph representation in Section 4.5. The spherical cap formula (4.11) holds true for these two pieces.
Recall the contact angle θ such as tan θ = −∂ r u| r=b . Then by elementary calculation we obtain the classical spherical cap volume formula in 3D and in 2D, the formula becomes

4.2.1.
Construct an exact breathing droplet solution and compare with numerical simulations. Motivated by the spherical cap solution, to check the long-time validation, we construct a breathing spherical cap solution with a prescribed oscillating contact angle θ(t) satisfying (4.14) where the parameters κ(t), σ(t) are determined below. Now we proceed to derive the formulas κ(t), σ(t) for this breathing droplet. Given θ(t) with oscillations, we will first calculate u(x, t) and b(t) from the spherical cap formula and then find κ(t) and σ(t) such that the PDE (4.14) holds with the Lagrangian multiplier λ(t).
Step 2. Calculate u(x, t) and b(t). From the spherical cap formula (4.13) , This construction automatically preserves the volume V and by elementary calculations, we know the following relations Computed by the first order scheme in Section 3.2.1 with κ(t) and σ(t) in (4.18) and oscillating contact angle θ(t) = θ in +A sin t, with θ in = 3π 16 , A = 0.2. The parameters in the PDE solver are β = 0.1, T = 30π, ∆t = 0.0628, N = 1000 for moving grids, initial domain [−3, 3] and initial u(x, 0) calculated by (4.16). Each subfigure shows the breathing droplet at time snapshots [0, π 2 , π, 3π 2 ] and the recurrence after 15 periods. Exact solution at 29π is shown using dotted orange line as comparison. Step 3. We find κ(t), σ(t) and λ(t) such that (4.14) holds. From the (4.17), upto some elementary calculations, we have (4.18) Particularly, for the quasi-static case β = 0, we have The constructed breathing droplet solution can be easily extended to 3D case using (4.12). Let the oscillating contact angle be θ(t) = θ in + A sin t, with θ in = 3π 16 , A = 0.2. Now we show the evolution of breathing droplet and the periodic recurrence for [0, 30π]. The dynamics of the breathing droplet in Fig. 2 is computed by the first order scheme in Section 3.2.1 with κ(t), σ(t) in (4.18) and with initial wetting domain [−3, 3] and initial profile calculated by (4.16) for t = 0. The parameters in the PDE solver are β = 0.1, final time T = 30π, time step ∆t = T 1500 = 0.0628, N = 1000 for moving grids in (a(t), b(t)). The exact solution at T = 29π based on (4.16) is also shown with dotted orange line in the lower left part of Fig. 2 as a comparison. 4.3. Capillary motion of a droplet in a Utah teapot. The Utah teapot is an important object in computer graphics history, whose 2D cross section can be completely described by several cubic Bézier curves [6]. In this section, we will use the bottom and the mouth of the Utah teapot as the inclined substrate to demonstrate the competition between the gravitational effect and capillary effect for droplets with small Bond number. We use four points (x i , y i ), i = 1, · · · , 4 to construct a cubic Bézier curve (x( ), y( )) with parameter ∈ [0, 1]. Denote the Bernstein basis polynomials as (4.20) Then the cubic Bézier curve is uniquely given by Now we construct the bottom and the mouth of the Utah teapot using 10 points x i , y i , i = 1, · · · , 10 listed in Table 3. For the bottom of the teapot, we use (x i , y i ) for i = 1, · · · , 4 and (x i , y i ) for i = 4, · · · , 7. For the mouth of the teapot, we use (x i , y i ) for i = 7, · · · , 10. Notice the inclined rough substrate is now expressed by parametric curve (4.21). Let (x) be the inverse function of x( ), then w(x) = y( (x)) and θ 0 = 0 in (2.22), which means we do not rotate the Cartesian coordinate system. To evaluate function w at endpoint a in the numerical implementation, one can use linear interpolation a = (1 − α)x( i ) + αx( i+1 ) for some α ∈ [0, 1]. Now we take the physical parameters as κ = 5, β = 3 and the initial droplet as To see clearly the instantaneous contact angle dynamics, we also track the contact angle at a(t) and b(t) for both the rolling down case (left subfigure) and the rising up case (right subfigure) in Fig.  4 associated with droplets dynamics in the teapot. Those contact angles vary as the effective slope of the mouth of the teapot changes, and then tend to the equilibrium Young's angle θ Y . Particularly, two cusps occur for the rolling down case (left subfigure) in Fig. 4 when the advancing (resp. receding) contact line a(t) (resp. b(t)) pass through the corner between the mouth and the body of the Utah teapot. The red line is the dynamics of the contact angle θ a at a(t) while the blue line is the dynamics of the contact angle θ b at b(t). One can see at the late stage, the sign of both a (t) = 1 cos θ 0a (cos θ a − cos θ Y ) and b (t) = − 1 cos θ 0b (cos θ b − cos θ Y ) are negative (resp. positive) for the rolling down case (resp. rising up case).

4.4.
Dynamics of droplets on an inclined groove-textured surface. In this section, we show the contact angle hysteresis (CAH) for a droplet on an inclined rough surface. Gravity will pull the droplet down while CAH will resist its motion. Therefore, one will observe the top of the droplet becomes thin while the bottom of it becomes thick. Besides, the contact line speeds depend on both the instantaneous contact angle θ a , θ b and the local slope of the rough surface θ 0a , θ 0b (in (2.22)), which changes constantly due to the boundary motion. Consequently, one can observe the contact line speed will change accordingly.

4.5.
Quasi-static dynamics of Kelvin pendant drop with volume constraint. In 1886, Lord Kelvin proposed a geometric integration procedure that the quasistatic profile remains no long graph representation and becomes "repeated bulges" when the height of the pendant drop exceeds a critical height u c . In this section, we compute the Kelvin pendant drop problem with volume constraint for κ < 0, which is not covered in Proposition 3.4. For the cases without volume constraint, we refer to [45]; see also [26,48] and references therein. To simulate the "repeated bulges", which certainly break the vertical single graph representation setting, we will first describe the droplet using inverse function X(u) (in horizontal graph setting) and give the gradient flow formulation in terms of X(u). By solving the DAEs for X(u) with κ < 0, which describes the quasi-static dynamics of a pendant droplet, we will recover multiple interfacial shapes including lightbulb and hourglass shapes with different Bond numbers. We refer to [45] for simulations and stability analysis of a liquid drop in hydrostatic states. For the case the capillary surface can not be expressed uniquely by the graph function u(x), we use a horizontal graph setting X(u). From [53, Theorem 1.1], for the quasi-static dynamics of a pendant droplet, there is only one vertical axis of symmetry such that any nonempty intersection of u with a horizontal hyperplane are disks centered at that vertical axis. Thus the maximum u m of the droplet is attained uniquely at X(u m ) = 0 and there exists a unique horizontal graph representation X(u), 0 ≤ u ≤ u m . For simplicity, we also assume the full dynamics of a pendant droplet has a horizontal graph representation in the derivation of (4.26). Now we identify the droplet on the right of x = 0 as  Next we give the following governing equations for a 2D droplet in terms of X(u) (4.26) The derivation using a gradient flow on a Hilbert manifold is similar to (2.15); see Appendix D.
To compute the Kelvin pendant droplet problem with volume constraint, we consider the quasistatic dynamics by taking β = 0 in (4.26). After desingularization, the quasi-static dynamics can be recast as the following DAEs for (X(0, t), u m (t), θ(t), λ(t)) (4.27) ∂ t X(0, ·) = − cos θ − σ, J(u m , θ) = 1,  κu−λ = 0 provided κu m = λ. After solving (X(0, t), u m (t), θ(t), λ(t)) from the above DAEs, we can further compute the formula for X(u, t) We use the DAEs solver ode15s in Matlab to solve the DAEs (4.27) with different initial data θ in , u m (0) and physical parameters κ, Young's angle θ Y ; see Table 4. We start the DAEs by first solving the compatible initial data b(0) = X(0, 0) and λ(0) from (4.27). At the final time T = 4, the final bulge shape of the pendant droplet (blue line) is illustrated in (left) Fig. 6, while the final lightbulb shape of the pendant droplet is illustrated in (right) Fig. 6. The corresponding Bond numbers and the final contact points b(4) are also shown in Table 4.
It is natural to assume the motion of the droplet η(t) is modeled by a gradient flow on manifold M described above. (i) The dynamics is driven by the free energy F(η) on manifold M; (ii) The dissipation mechanism of the dynamics is described by a Riemannian metric g η on the tangent plane T η M, which is discussed in (A.8) below. We will use the vertical velocity v = ∂ t u and the contact line velocity Γ t = v cl n cl to describe the tangent plane T η M. Since the geometric motion has an obstacle condition u ≥ 0, manifold M has a boundary, i.e., {η ∈ M; u(x, y) = 0 for some (x, y) ∈ D} (when the droplet has a splittingtype topological change). On the boundary, the tangent plane is not a linear space and has the restriction described below. Notice where we used the fact that ∇u(Γ(t), t) · n cl = −|∇u(Γ(t), t)| in the graph representation. The tangent plane at η is given by n| Γ = (sin θ cl n cl , cos θ cl ), v n | Γ = sin θ cl v cl . Now we describe the dissipation mechanism of the dynamics. From Rayleigh's dissipation function (2.9), since 2Q is the rate of energy dissipation due to friction [31], it is natural to introduce the Riemannian metric g η on T η M × T η M below. For any For similar derivations of the dynamics of droplets using a variational approach with various free energies and the same Riemannian metrics (A.8), we also refer to Davis [13], [46], Doi [56] and [32]. We now derive the gradient flow of F(η) defined in (2.7) on manifold M with respect to the Riemannian metric g η . For an arbitrary trajectoryη(s) = {Γ(s),ũ(x, y, s)} (physically known as a virtual displacement) passingη(t) = η(t) at the tangent direction q :=η (t) = {ṽ cl ,ṽ} ∈ T η(t) M, we know (A.9)ṽ(Γ) = |∇u(Γ(t), t)|ṽ cl .
Hence by taking differentη ∈ T η(t) M, the governing equations for u(·, t) ∈ H 1 0 (D(t)) and λ(t) are (A.14) with initial data η(0) = {Γ(0), u(x, y, 0)} and initial volume V . The variational inequality formulas above are able to describe the merging and splitting of several drops by using numerical schemes for parabolic variaional inequalities (PVI), such as splitting method with projecting operators, c.f. [34,40,29]. However, after the splitting of two droplets, the original two-phase interface becomes three-phase triple point at the splitting point, thus the purely PVI dynamics can not describe the contact line dynamics for the emerged triple point, i.e., new contact lineΓ. Indeed, according to PVI, the motion of the emerged triple pointsΓ is only guided by the projection of the motion by mean curvature on the substrate, thus is different from the true physical case, i.e., Rv cl = − [G + |∇u|(n cl · ∂ ∇u G)] Γ . Hence we need to detect the splitting points, enforce this contact line dynamics for newΓ and then treat them separately as two droplets. Existence of global weak solutions including topological changes (splitting and merging) is studied in [32] using a continuum limit of a time discretization based on variational methods.
A.1. Gradient flow for a single droplet: without merging and splitting. We call a single droplet as a sessile drop if u(x, y, t) > 0 for (x, y) ∈ D(t) with the gravity downwards, i.e., g > 0.
Another scenario is when a light drop is in a heavy fluid, the drop experiences buoyancy due to gravity. In this case, we call a single droplet as a pendant drop if u(x, y, t) > 0 for (x, y) ∈ D(t) with the gravity upwards, i.e., g < 0. Another equivalent problem is a drop pendant on ceiling with u < 0 for (x, y) ∈ D(t) and g > 0. In this paper, we only consider nonnegative u and use negative g for a pendant droplet. For those single sessile/pendant drop cases, the variational inequalities become variational equalities and the weak formulations can be equivalently converted to a strong-form PDE. Therefore the governing equations with volume constraint (A.14) become where V is the initial volume of the droplet. Particularly, for the energy (2.7), we have Therefore the governing equations are (2.10). We remark that when G(u, ∇u) = 1 2 |∇u| 2 + σ, the kinematic boundary condition for the contact line speed v cl in (A.14) becomes R γ lg v cl = 1 2 |∇u| 2 − σ, which recovers the kinematic boundary condition used in [32,35,22].

Proof of Proposition 2.1. Statement (i) comes from (A.3) directly.
For Statement (ii), from the Reynolds transport theorem and similar to (A.12), we have d dt D(t) G(u(x, y, t), ∇u(x, y, t)) dx dy which together with (A.3) and (A.15), gives For (iii), we derive the the gradient flow for a single droplet with quasi-static dynamics. If we consider the gradient flow in the quasi-static setting, i.e., ζ = 0, we can regard u(x, y, t) as being driven by Γ(t). In other words, we consider {Γ(t), u(x, y, t)} with u as the solution to This gives a reduced manifold Γ(t). Correspondingly, we have the quasi-static trajectory η qs (t) := Γ(t), the quasi-static free energy F qs (Γ(t)) := F((Γ(t), u(x, y, t))). and the quasi-static tangent plane T ηqs . With the quasi-static metrics g ηr (η r ,η r ) = R Γ(t) v clṽcl ds, we have the gradient flow for quasi-static dynamics In 2D, the units of γ lg becomes energy/length, R becomes mass/time, ζ becomes mass/(length·time) and ς is 1/(length 2 ). Let t = Tt, x = Lx, u = Lû, a = Lâ, b = Lb, V = L 2V and λ = γ lg Lλ where L is the typical length of the droplet and T is the typical time scale to observe the motion of contact lines. In other words, we choose typical time T such that RL γ lg T = 1 and typical volume for unit disk V = π. We denote the capillary number for the capillary surface as β := ζL 2 γ lg T and set κ = L 2 ς, both being dimensionless. Then the dimensionless equations (after dropping hat) in 2D are given by (2.15).
A.2. Additional hydrodynamic effects inside the droplets. We also remark the relations between our pure geometric motion of droplets and the other contact line dynamics coupled with hydrodynamic effect of viscous fluid.
Consider further the hydrodynamic effect, which specifically is (i) adding kinetic energy 1 2 ρ A |v| 2 dx due to inertial effect in the free energy F, (ii) adding viscosity dissipation inside the droplet and (iii) adding energy dissipation on solid-liquid interface due to Navier slip boundary condition. Then the energy dissipation relation becomes [47, eq(39)] where α is the slip length. The corresponding governing equations are Rv cl = γ lg (cos θ Y − cos θ cl ), on Γ(t).
Here, v is the fluid velocity, p is the pressure, µ is viscosity and ρ is the fluid density. Notice in the third equation in (A.21), the normal traction (normal force) induced by the fluids F n = −p + n T · µ(∇v + ∇v T ) · n = γ lg H is balanced with the mean curvature γ lg H. Meanwhile, in the purely geometric motion, the normal velocity of the capillary surface ζv n = −γ lg H can be understood as velocity induced by an underlying normal frictional force F n = −ζv n . Roughly speaking, the purely geometric motion derived in Appendix A captures the same main feature (motion by mean curvature of capillary surface) as the original hydrodynamic one in [47], where the normal velocity of the capillary surface is induced by the fluid velocity following Navier-Stokes equation.
Recall the first equation in (2.10) in the absence of gravity, i.e., The energy dissipation relation (2.12) is exactly same as the dissipation relation in [47, eq(38)], where v is the velocity of the capillary surface. Here λ(t) does not contribute due to ∇ · v = 0 inside the droplet.
Appendix B. Proof for the truncation analysis of first and second order schemes In this section, we give some truncation error estimates for the first and second order schemes in the case w(x) = 0, θ 0 = 0 and thus h(x) = u(x). Now the governing equation (2.22) becomes (2.15).
Proof of Lemma 3.7 (first order truncation error estimates). Let a(t n+1 ), b(t n+1 ), u(x n+1 , t n+1 ) for x n+1 ∈ [a(t n+1 ), b(t n+1 )] be the exact solution to (2.15) evaluated at t = t n+1 with initial data at t = t n , a n , b n , u n (x n ) for x n ∈ [a n , b n ]. We outline the idea of proof below.
Proof of Lemma 3.8 (second order truncation error estimates). Let a(t n+1 ), b(t n+1 ), u(x n+1 , t n+1 ) for x n+1 ∈ [a(t n+1 ), b(t n+1 )] be the exact solution to (2.15) evaluated at t = t n+1 with initial data at t = t n , a n , b n , u n (x n ) for x n ∈ [a n , b n ]. We will prove truncation error estimates (3.77) and (3.78) separately in Step 1 and Step 2. We outline the idea of proof below.
Step 1. Second order truncation error for the moving boundary (3.77).
We first illustrate the idea of the usual truncation error estimates for the predictor-corrector ODE solver for v = f (v) with v = f (v)f (v). By Taylor expansion, which is equivalent to Moreover, for any smooth function W (v), we have the second order estimate (B.18) also gives another method for second order truncation error estimate by evaluating the equation v = f (v) at t n+ 1 2 , which is (B.17). The truncation error in Step 2 will rely on this method.
Finally, we proveã n+1 −a n ∆t − (a ) n and ∂ ZŨ n+1 −∂ Z U n ∆t − (∂ Zt U ) n | Z=0 have O(∆t) accuracy. Since the predictorã n+1 is given by the first order scheme in Section 3.2, we knowã n+1 −a n ∆t − (a ) n has O(∆t) accuracy and we obtain (3.54). To estimate ∂ ZŨ n+1 −∂ Z U n ∆t , we give the following claim. Claim 1: Assume we have the error estimates Then we have the second order accuracy The proof of Claim 1 is based on changing moving domain to fixed domain by Z =x , which is similar to (B.20) and will be omitted. Notice the first order accuracy of predictor a n+1 ,b n+1 and we used implicit elliptic solver with second order accuracy in (3.59) for predictor u n+1 , so the assumptions in claim 1 are satisfied automatically. Thus from the Taylor expansion and claim 1 we know Therefore, we complete the second order truncation error estimates for the moving boundary (3.77).
Step 2. Second order truncation error estimates (3.78) for u n+1 . First from the similar argument for (B.18), we have the following generalized claim Claim 2: For any smooth function W (v(x, t), v x (x, t), v xx (x, t), x, t), we have where the equality holds in the sense of changing variables to fixed domain Z = x−a(t) b(t)−a(t) with the relation (B.5).
Using further Claim 1 and Claim 2, we obtain the second order accuracy for (3.78). Therefore, we complete the second order truncation error estimates for (3.78).
6. Update the moving grids x n+1 j = a n+1 + jτ n+1 , τ n+1 = b n+1 −a n+1 N , j = 0, 1, · · · , N. 7. Calculate h n * j based on (3.76), for j = 1, · · · , N − 1. 8. Solve h n+1 implicitly. For j = 1, · · · , N − 1, with h n+1 Appendix D. Gradient Flow for single sessile drops in non-wetting case using horizontal graph representation X(u) Recall the description of the non-wetting droplet in terms of X(u) in (4.25). We consider the manifold based on X(u) Similar to Appendix A, we calculate the gradient flow on manifold M. Below, we directly use dimensionless quantities for simplicity. LetX(u, s) withũ m (s) be trajectory on M starting from X(u, t) with u m (t). Based on the relations in Lemma 3.3, the tangent plane T η M at η can be described by the contact line speed v cl = X t (0) and the horizontal velocity v = X t of the capillary surface. Notice from X(u m (t), t) = 0, on the tangent plane T η M we have (D.2) ∂ tX = −X u ∂ tũm on u = u m .
For any q 1 = (v cl1 , v 1 ), q 2 = (v cl2 , v 2 ) ∈ T η M, define the Riemannian metric g η : Consider the free energy where we used (D.5). Then the gradient flow of E on manifold M with respect to the metrics g η gives the governing equations (4.26).
Remark 4. Using the similar derivations above, we have the governing equations for a 3D axisymmetric droplet in terms of R(u, t) (D.6) where θ is defined as tan θ = −R u | u=0 . Here , is the mean curvature in terms of R(u).