A linear, second-order, energy stable, fully adaptive finite-element method for phase-field modeling of wetting phenomena

We propose a new numerical method to solve the Cahn-Hilliard equation coupled with non-linear wetting boundary conditions. We show that the method is mass-conservative and that the discrete solution satisfies a discrete energy law similar to the one satisfied by the exact solution. We perform several tests inspired by realistic situations to verify the accuracy and performance of the method: wetting of a chemically heterogeneous substrate in three dimensions, wetting-driven nucleation in a complex two-dimensional domain and three-dimensional diffusion through a porous medium.


Introduction
Capillarity and wetting phenomena, driven primarily by interfacial forces, are ubiquitous in a wide spectrum of natural phenomena and technological applications. Examples range from the wetting of plant leaves by rainwater and insects walking on water to coating processes, inkjet printing, oil recovery and microfluidic devices; for reviews, see e.g. [17,5]. From an historical point of view, two of the concepts essential to the understanding of capillarity and wetting were introduced and studied already in 1805: these are the Laplace pressure [40] and the Young-Dupré contact angle [65]. Later, following the work of Plateau on soap films [45], Poincaré [46] linked interfacial phenomena with the theory of minimal surfaces.
Wetting phenomena typically involve a fluid-fluid interface advancing or receding on a solid substrate and a contact line formed at the intersection between the interface and the substrate.
The wetting properties of the substrate determine to a large extent the behaviour of the fluids in the contact-line region, and in particular the contact angle at the three-phase conjunction, defined as the angle between the fluid-fluid interface and the tangent plane at the substrate. At equilibrium, this is precisely the Young-Dupré angle. When one of the two fluids moves against the other, the contact angle becomes a dynamic quantity, and when the problem is formulated in the framework of conventional hydrodynamics, the contact line motion relatively to the solid boundary results in heterogeneities, or both.
Like the OD2 scheme on which it was based, the time-stepping scheme we propose is semiimplicit and linear. We show that it is mass-conservative and satisfies a discrete free-energy law with a numerical dissipation term of order 2 in time. Space discretization is achieved using a finite-element method, leading to an unsymmetrical sparse linear system to solve at each iteration.
We use a mesh refinement strategy to capture interfaces precisely, and an adaptive time step to limit the variation of free energy at each step, with the aim of increasing the resolution in time during fast phenomena.
To test the efficiency of the proposed numerical scheme we consider several wetting problems as test cases. We first study relaxation towards equilibrium in two situations: the spreading of a sessile droplet and the coalescence of two sessile droplets on a flat, chemically homogeneous substrate. We then consider two-component systems in complex geometries delimited by chemically heterogeneous substrates in both 2D and 3D.
In Section 2, we introduce the CH system and the non-linear wetting boundary condition. In Section 3, we outline our numerical scheme and prove the associated conservation properties. In Section 4, we present the results of several numerical experiments. Conclusions and perspectives for future work are offered in Section 5.
2 Phase-field model for wetting phenomena Throughout this study, Ω ⊂ R d corresponds to a d-dimensional domain, ∂Ω denotes its boundary with outward unit normal vector n, Γ S is the solid substrate and Γ G = ∂Ω \ Γ S . The CH system we use to describe the dynamics of two immiscible fluids in contact with a solid substrate, is a freeenergy-based model. The starting point is the introduction of a locally conserved field, denoted by φ : Ω → R, that plays the role of an order-parameter: two equilibrium values, say +1 and −1, represent the pure phases, and the interface is conventionally located at the points where φ = 0 [12,11]. We consider systems with a free energy given by where the two terms, E m and E w , represent the mixing and wall components of the free energy, respectively. Here F m (φ) = 1 4 (φ 2 − 1) 2 and F w is taken to be a cubic polynomial, following e.g. Refs [55,54]: where θ = θ(x) is the equilibrium contact angle, which can depend on the spatial position x. From the expression of the free energy, we calculate that, for a sufficiently smooth function ψ : Ω → R: with f m = F m and f w = F w , so the chemical potential is equal to and the natural boundary condition associated with the surface energy is We assume that the dynamics of the system is governed by the CH equation, where b(x) is a mobility parameter, assumed to be uniform hereafter. This leads to the following mass-conservation property: so the mass flux at the boundary can be specified using the condition b∇µ · n =ṁ(x), wherė m(x) is the desired mass flux. In particular, we will setṁ(x) = 0 at the solid boundary, Γ S . In summary, the equations we are solving in this study are the following: In addition to the conservation of mass, Eqs. (9a) to (9d) imply the following energy-conservation law, involving the phase field and the chemical potential: An advantage of the cubic surface energy (3) over other surface energy formulations (see [33] for a review of wetting boundary conditions for binary fluids) is that the well-known hyperbolic tangent profile is an equilibrium solution in more than 1 dimension. Specifically, the function is solution to the CH equation posed in the half plane {y ≥ 0} with the boundary condition (3) at {y = 0} and constant θ(x) = θ. A schematic representation of this solution and the corresponding fluid-fluid interface is given in Fig. 1.
A drawback of the cubic wall energy (3) is that the conservation of energy no longer seems to imply stability bounds for the solution, making it impossible to use the tools traditionally employed (see e.g. [21]) to prove the well-posedness of the system. Indeed, an application of the I n t e r f a c e Figure 1: Schematic of the profile geometry of a fluid-fluid interface intersecting a solid boundary and illustration of the stationary solution (11).
trace inequality gives only that, under appropriate regularity assumptions on φ: where we used Hölder's inequality and Young's inequality with a parameter. Therefore, the wall energy cannot be controlled by the mixing energy for arbitrary domains. This issue can be remedied by a simple modification of the wall energy outside of the physical range [−1; 1]; instead of (3), we consider the following wall energy: This function is such that , and (F * w ) is absolutely continuous, which makes it possible to prove the second order convergence of our time-stepping scheme, see Section 3. Another possibility would have been to choose constant values for F * w outside of the interval [−1; 1], but this would have lead to F * w being only C 1 (R), making it more difficult to show second order convergence theoretically. The weak formulation of Eqs. (9a) to (9d) with the modified wall energy (13) is as follows: find (φ, µ) such that φ ∈ L ∞ (0, T ; H 1 (Ω)), ∂φ ∂t ∈ L 2 (0, T ; (H 1 (Ω)) ), µ ∈ L 2 (0, T ; H 1 (Ω)), (14) and the following variational formulation is satisfied: with f * w := (F * w ) and where ·, · , (·, ·) and (·, ·) ∂Ω denote respectively the duality pairing between (H 1 (Ω)) and H 1 (Ω), the standard inner product in L 2 (Ω), and the standard inner product in L 2 (∂Ω). For simplicity of notations, the symbols F w , f w and E will refer in the rest of this paper to F * w , f * w , and E m + ∂Ω F * w dσ, respectively.

Well-posedness of the model
We can show the following existence result for the weak formulation of the Cahn-Hilliard system with the modified boundary condition presented above, under appropriate regularity assumptions for the initial condition and the mass fluxṁ.
Proof. See Appendix A.

Numerical method
In this section we introduce a new time-stepping scheme to solve the CH equation (7) with the non-linear wetting boundary condition (6), which is a generalisation of the optimal dissipation scheme of order 2, OD2, developed in [27]. We decided to extend this particular scheme because the authors of [27] showed that, among all the linear schemes they proposed, it is the most accurate and the least dissipative. And in selected test cases, they showed that for a large enough time step, it is the only scheme that leads to the correct equilibrium solution. We refer to our scheme as OD2-W, with W denoting wetting, and show that it leads to a consistent discrete energy law.
We also develop a new adaptive time-stepping strategy which, combined with adaptation in space, leads to a fully adaptive finite-element method. An excellent introduction to the finiteelement method and corresponding mixed formulations can be found in Ref. [13] and to mesh generation and adaptive refinement in Ref. [24].

OD2-W scheme
In this section, we assume for simplicity thatṁ = 0 and that θ is uniform on ∂Ω. We denote by ∆t the time step, and by φ n and µ n+1/α the numerical approximations of φ and µ at times t n and t n + 1 α ∆t, respectively. To define a discretization in time of the CH system appropriate for wetting phenomena, we follow the approach proposed in [27] to design an optimal dissipation scheme, and consider the following generic implicit-explicit numerical scheme: given φ n ∈ H 1 (Ω), In these expressions,f m ,f w are functions to be specified, linear in their second argument. The parameter α ∈ {1, 2} determines the accuracy of the numerical scheme, and the parameter β ∈ [0, 1−1/α] controls the numerical diffusion. The function φ n+ 1 α +β is defined by linear interpolation between φ n and φ n+1 , and δ t φ n+1 is the approximation of the time derivative of φ given by In most numerical experiments presented in this paper, we consider the case (α, β) = (2, 0) (OD2-W), but we note that other usual choices include (α, β) = (1, 0) (OD1-W) and (α, β) = (2, O(∆t)) (OD2mod-W). By taking ψ = µ n+ 1 α and ν = δ t φ n+1 in (16), we obtain where N D(φ n , φ n+1 ), representing the non-physical numerical dissipation introduced by the timestepping scheme, can be broken down in three parts: with Notice that the philic dissipation is always nonnegative, with N D philic (·, ·) = 0 if (α, β) = (2, 0) if (α, β) = (1, 0) (OD1-W). The two other terms can be expanded using Taylor's formula, taking into account that F m is a polynomial of degree 4 and using the integral form of the remainder: This suggests the following choices for the functionsf m andf w : where the last expression is convenient for programming purposes. We note that this methodology to derive a second-order scheme can be applied mutatis mutandis when using the unmodified wall energy (3), although we haven't been able to prove that the weak formulation is well-posed in that case. Doing so leads tof w (φ n , φ n+1 ) = −( √ 2/2) cos θ (1 − φ n φ n+1 ), which coincides with (23b) when φ n ∈ [−1, 1]. In either case, we have the following property: Property 3.1. Assume that α = 2 and β = 0. Then the numerical dissipation term in Eq. (19) is such that , provided that all the terms in the definition of C are well-defined.
Proof. In [27], the authors show that: For the wall term, we obtain from Eqs. (22b) and (23b): In addition to the energy law (19), the numerical scheme (16) satisfies a discrete version of the conservation law (8) presented in Section 2, which can be seen by choosing ψ = 1 in Eq. (16a).
Property 3.2. The numerical solution satisfies the following mass conservation law:

Space discretization and adaptive mesh refinement
Our approach for mesh adaptation is based on a method proposed in [30], and implemented through the FreeFem++ functions adaptmesh (in 2D) and mshmet (in 3D). The idea of the method is to define a metric on the computational domain based on the solution at the current time step, and to use for the next time step a mesh that is uniform in that metric. The metric we consider corresponds to the following metric tensor, depending only on the phase field φ: where (λ i (x)) d i=1 are the eigenvalues of the Hessian of φ at x, R(x) is the matrix containing the associated orthonormal eigenvectors, and γ > 0 is a parameter controlling the interpolation error. A standard algorithm of Delaunay type is used to generate a mesh that is equilateral and uniform with characteristic length 1 in that metric. This mesh definition ensures that the interpolation error of the phase field is roughly equi-distributed over the parts of the domain where In most of the simulations presented in the next section, we set h min to a value lower than or equal to ε/5, to ensure that enough mesh points are available for the discretization of the interface region in its normal direction, and h max to a value small enough that a good approximation of the chemical potential is possible. For 3D simulations, however, choosing h min ≤ ε/5 when ε is of the order of 0.01 leads to a prohibitive computational cost; in these cases we have thus used a less precise value, as specified in the relevant sections.
For a given mesh T = with P ρ the space of polynomials of degree ρ. In the numerical experiments below, we used both quadratic elements (ρ = 2) and linear ones (ρ = 1). Space discretization is achieved by replacing (16), leading to a sparse unsymmetric linear system at each iteration.

Time step adaptation
Here we assume thatṁ = 0 in the boundary condition (9d). From Eqs. (8) and (10), this implies that M (φ) is constant in time and E(φ) decreases. Numerical exploration suggests that large freeenergy variations are usually caused by topological changes of interfaces, corresponding to physical phenomena such as the coalescence of droplets. Since capturing such phenomena precisely is crucial to the accuracy of the solution, we propose an adaptive strategy aimed at limiting the variation of free energy at each time step. We adapt the time step based on the dissipation of free energy, parameters enter in our time-adaptation scheme: • ∆t min , ∆t max : the time steps below which we stop refining and beyond which we stop coarsening, respectively.
• ∆E min : the variation of free energy below which we increase the time step at the next iteration.
• ∆E max : the variation of free energy beyond which we refine the time step and recalculate the numerical solution.
• f > 1: the factor by which the time step is multiplied or divided at each adaptation.
The condition (E(φ * ) − E(φ n )) > ∆E max /100 serves to guarantee that the method does not blow up. The choice of a nonzero right-hand side is motivated by the fact that, when the system is close to equilibrium, it can happen that E(φ * ) > E(φ n ). This is because, in contrast with the sign of N D philic (φ n , φ n+1 ), which is always positive or zero according to Eq. (21), the signs of N D phobic (φ n , φ n+1 ) and N D wall (φ n , φ n+1 ) are in general unknown.
In the numerical experiments presented in Section 4, we chose ∆t min = 0. Since the numerical dissipation term scales as ∆t 2 , the inequality E(φ n+1 ) ≤ E(φ n ) + ∆E max /100 will always hold for ∆t small enough, so the refinement process is guaranteed to terminate at each iteration.

Numerical results
The new numerical method is applied on a number of test cases. For the implementation, we have used FreeFem++ [29] for the implementation of the finite-element method and 2D mesh adaptation, umfpack [16] for the linear solver, mshmet [23] and tetgen [51] for the mesh adaptation in 3D, and gmsh [26] for the description of the geometry, post-processing and 3D visualisation. In Section 4.1 we check that the numerical scheme leads to the correct equilibrium solution in the simple case of a droplet spreading on a philic or phobic substrate. In Section 4.2 we study the convergence of the method with respect to the time step and the mesh size, when a uniform mesh and a constant (a) θ = π/3, θ * /θ = 0.991.

Equilibrium contact angle
We consider a 2D sessile droplet on a flat substrate where we impose the no-flux condition and the wetting condition (6) with the modified wall energy (13) and uniform contact angle θ: Our aim in this section is to check that our method is able to accurately capture the imposed contact angle, θ. Figure 2 shows the equilibrium position of a droplet for different values of θ, for b = 1 and ε = 5 × 10 −3 . In all cases we used the scheme OD2-W with adaptation in space using the parameters h max = 10 h min = 0.01, and we computed the contact angle of the φ = 0 isoline at the substrate. A very good agreement is achieved between the imposed equilibrium contact angle and the observed numerical one.

Convergence of the method
Here, we study the convergence of the method when both time step and mesh size decrease. The problem we considered to that purpose is the coalescence of two adjacent sessile droplets as they spread on a flat substrate. For the simulation, we used the initial condition   mesh adaptation and for ε = 0.1, so that enough data points could be generated at a reasonable numerical cost. Since the exact solution to the CH equation in this case is not known analytically, we calculated the error by comparison of the numerical solutions to the solution obtained with the smallest value of h. Results are presented in Fig. 3. As we can see, the observed convergence rate is almost equal to 2, which is the optimal rate in the case of linear basis functions. Now we address the convergence with respect to the time step. For this case, we used the parameters ε = 0.1, b = 10 4 , and the minimum time step we considered was ∆t * := 0.00665. In Fig. 4, we present convergence curves for OD1-W, OD2-W, and OD2mod-W. We note that the convergence rates are close to the expected ones, and that the use of OD2 gives significantly more accurate results than the other two methods. In Fig. 5, the total numerical dissipation produced by the numerical schemes is presented. Here too, numerical results agree with the theoretical results of Section 3. (c) OD2mod-W with β = 10 ∆t Figure 5: Total numerical dissipation generated by the numerical schemes in the simulation used to produce Fig. 4. OD2-W is by far the scheme producing the least numerical dissipation, even for relatively large time steps. OD2mod-W, on the other hand, introduces significant numerical dissipation for large time steps, owing to the large value of β that was chosen for the simulation, but is much less dissipative than OD1-W for small time steps.

Time-adaptation scheme
In this section, we examine the performance of the adaptive time-stepping scheme in the case of two droplets evolving on a chemically homogeneous substrate. We start from the situation where The evolution of the time step, of the number of recalculations, and of the free energies is displayed in Fig. 8. In all three cases, the time step is refined several times at the first iteration, to accommodate for the discontinuity of the initial condition. Since the initial angle between the interface and the substrate is equal to π/2, the number of recalculations performed at the first iteration is higher for θ = π/4, 3π/4 than for θ = π/2. After the initial refinement, the time step steadily increases to its maximum allowed value for θ = π/2 and θ = 3π/4, but when θ = π/4 a second refinement occurs to capture the coalescence of the droplets.
In this latter case, we observe, simultaneous with the second refinement of the time step, an increase in the rate of dissipation of free energy. After the formation of a new stable interface, the total free energy continues to decrease, but more slowly, as a new droplet, formed by the merging of the two original droplets, moves towards its equilibrium position. We clearly identify the coalescence time by looking at the singularity in the curve corresponding to the mixing energy.
This energy increases before coalescence, as the interfaces are being stretched, and it decreases steadily after. The wall energy, on the other hand, decreases at first and increases in the later stage of the simulation. As prescribed by Algorithm 1, the time step detects the variations of free energy; it decreases when the rate of variation of the total free energy increases, and conversely.
For comparison purposes, we also included in Figs. 8d to 8f data corresponding to the case where a fixed time step is used for the simulations presented in this section. There does not  currently exist any result with conditions on the time step that ensure the stability of OD2, and we haven't been able to show stability results for OD2-W either. In practice, we observed that the time step required to ensure stability of OD2-W with the set of parameters we use in this test case would lead to a very high computational cost. We point out that, contrary to what we expected, the time step required to achieve stable integration in time with the modified wall energy (13), which we use in this paper, seems to be generally smaller than with the cubic formulation (3). To keep the computational cost at a reasonable level, we carried out the simulations with a fixed time step using the method OD1-W, the greater stability of which enabled us to choose ∆t = 0.02. In Figs. 8d to 8f, we see that, for the same contact angle, the curves corresponding to a fixed and an adaptive time step are almost undistinguishable. The agreement is also very good at the level of the phase field and chemical potential, although we do not present snapshots of the solutions obtained with a fixed time step.
The CPU times corresponding to the three contact angles considered are presented in Table 1.
Without adaptation, the simulations take significantly longer to run, which is consistent with the fact that more iterations (20000) were necessary to reach the final time. In addition, among the simulations that used an adaptive time-step, the difference between the CPU times is also significant, with the case θ = π/4 taking more than twice as long as the case θ = π/2.  Table 1: CPU times (hh:mm:ss) using an Intel i7-3770 processor for the simulations presented in Section 4.3 (two droplets on a substrate), with or without time-step adaptation. The method OD2-W was used for the simulations with an adaptive time step, and the method OD1-W was used for the simulations with a fixed time step. In both cases, an adaptive mesh was used, with the parameter h min equal to ε/10 = 0.001.

Wetting in complex geometries and with heterogeneous substrates
We now present the results of numerical experiments in more complicated and realistic settings, in both 2D and 3D systems.

3D droplet on a chemically heterogeneous substrate
We study the dynamics of a 3D sessile droplet on a flat substrate with chemical heterogeneities, i.e. the contact angle has a spatial dependence now, say θ = θ(x, y). This situation typically arises in electro-wetting settings [42]. It is widely accepted that the droplet shape can be controlled using patterned substrates, e.g. Ref [61,50], that may also be modelled efficiently using a space varying contact angle [61]. We consider chemical heterogeneities on the substrate of the form with θ 0 = π 2 the mean contact angle, a = π 6 the amplitude, and f x = f y = 4 the frequencies in x and y directions, respectively. As initial condition we take a droplet of base radius r 0 = 0.24   centered at x 0 = (0.5, 0.5, 0). The initial values of the phase field are given as Results are displayed in Fig. 9. The droplet, initially spherical, spreads on the hydrophilic regions of the substrate, and retracts from the hydrophobic patches. While we do not present any quantitative analysis of the error in this case, we note that the wetting behaviour agrees qualitatively with what one might expect intuitively from our understanding of wetting phenomena. While it progresses towards equilibrium, the droplet adopts a diamond-like shape.
For this test case, we used the method OD2-W with adaptation in space and time. The parameters used were the following: b = 10 4 , ε = 0.02, h max = 10 h min = 0.1, ∆t 0 = 0.0016, ∆t min = 0, ∆t max = 16 ∆t 0 , f = √ 2, ∆E max = 10 ∆E min = 0.0001. With these parameters, the time step was refined only at the beginning of the simulation, which is consistent with the absence of coalescence events in this case. There were 24 recalculations at the first time step, corresponding to a refinement of the time step by a factor f 24 = 4096.

Diffusion in a 3D porous medium
Here we consider a binary fluid in a model porous medium consisting of a cube filled with spheres.
The cube has edges of length 1, and the spheres have radius 0.1 and are located at positions Figure 10: Evolution of the isosurface φ = 0 of the phase field when a constant flux is imposed at the bottom boundary; The pictures correspond to iterations 0, 200, 400, 800 and 1000. Note that, because of the neutral boundary condition imposed at the spheres, the isosurface tends to stay normal to them as long as they are not completely covered.
(1.5, 1.5, 1.5) + 2∆(i, j, k) with ∆ = 1/7 and i, j, k ∈ {0, 1, 2}. We take all the substrates to be neutral, i.e. θ = π 2 , and the initial condition is the same as used before, defined by Eq. (33). In addition, we include an inflow boundary condition at the bottom of the cube to represent a pore where liquid can be pumped in. The radius of this pore is 0.1 and is located at (0.5, 0.5, 0). This boundary condition can be incorporated by imposing The imposed contact angle at the spheres is π/2, forcing the isosurface to stay normal to the spheres as long as these are not completely covered. Because of the boundary condition (34), the mass increases linearly, and the free energy increases, in agreement with Eqs. (8) and (10).
This case study demonstrates the ability of our method to easily tackle complex geometries. The parameters used for this test case are the same as in Section 4.4.1, except that we employed the fixed time step ∆t = 0.001.   (Figure 10). In this case, the mass increases linearly because we impose a constant mass inflow at the pore. The free energy increases as well, because the size of the interface grows, in agreement with both the mass and energy laws (8) and (10).

Nucleation processes with complex boundaries
The last problem we study is the process of phase separation in a domain with complex boundary characterized by different length scales. Specifically, we consider a domain defined by the coastline of the two islands that form the United Kingdom and Ireland. Starting from a satellite black and white picture, we extracted the isolines that define the contour of the different islands, which we passed to the FreeFem++ mesh generator to obtain a triangular mesh (for this, we based our code on a FreeFem++ example for the Leman lake). At the boundary we consider the contact angles θ = π/4, π/2, 3π/4, and we assume that the phase field is initially set to a random value at each grid point, drawn from a random normal distribution with variance 0.1. A fixed mesh was used for this simulation, and the parameters used were b = 1000, ε = 0.02, ∆E min = 0.02, ∆E max = 0.04, f = √ 2, ∆t min = 0, ∆t max = 1.
The evolution of the phase field and of the chemical potential in the case θ = π/4, obtained with the adaptive time-stepping scheme 1, is presented in Fig. 12. For each of the contact angles considered, we also ran a simulation with the fixed time step ∆t = 0.01, using the method OD1-W instead of OD2-W to benefit from the stabilizing effect introduced by the philic numerical dissipation of OD1-W. We note in particular that OD2-W is unstable for the selected value of ∆t, with oscillations appearing in the energy curves from the first iterations, and that the time step would have to be reduced significantly to ensure stability. The final configurations (time 500) are presented in Fig. 13 for the three contact angles considered. We observe that the final configurations are different depending on whether or not an adaptive time step is used, which can be attributed to the high sensitivity of the solution to perturbations of the initial condition chosen for this test case; the areas where separation of the phases first occurs is influenced by numerical errors in the early stages of the simulation.
Simulation data are presented in Fig. 14. With an adaptive time-stepping scheme, it appears from Fig. 14 (a) that, overall, the time step increases steadily as the frequency of coalescence events decreases. At specific times, the time step decreases slightly in order to accurately capture the evolution. As expected, the total free energy has a roughly constant negative slope when plotted Figure 12: Evolution of the phase field and chemical potential for the nucleation in a domain with complex boundaries, when starting from a random distribution. As before, blue corresponds to φ = 1 and green to φ = −1. The contact angle imposed at the boundaries is θ = π 4 .
against the iteration number. Here too, we observe a small discrepancy between the fixed and adaptive cases, which is consistent with differences observed at the final time in Fig. 13.
The CPU times corresponding to the simulations presented in this section are displayed in 2.
For the parameters selected, the adaptive time-stepping scheme leads to a lower computational cost. This test demonstrates the advantage of using a finite-element approach, as it would have  Table 2: CPU times (hh:mm:ss) using an Intel i7-3770 processor for the simulations presented in Section 4.4.3 (nucleation in a domain with complex boundaries), with or without time-step adaptation. The method OD2-W was used for the simulations with an adaptive time step, and the method OD1-W was used for the simulations with a fixed time step. In all cases, we used a fixed mesh with mesh size h = 0.01 (the size of the domain is roughly 5 by 5) and P 1 elements, leading to 181587 unknowns. The parameter ε was set to 0.02.
been very complicated to solve the CH equation in the geometries we consider here with e.g. a spectral method or finite differences.

Conclusions
We have proposed a new, fast and reliable numerical method to solve the CH equation with a wetting boundary condition. Our method is a generalisation of the OD2 scheme introduced in [27], which considered only the homogeneous condition ∇φ · n = 0. In addition, we have designed a new time-step adaptation algorithm, leading to a scheme that is adaptive both in space and time, and we have shown that this scheme is mass-conservative and satisfies a consistent discrete energy law.
We checked the validity of the proposed numerical scheme with several examples. First we (a) θ = π/4 (b) θ = π/2 (c) θ = 3π/4 Figure 13: Comparison of the solutions at the final time, with (left) and without (right) time-step adaptation. At the initial time, the phase field is set to a random value at each grid point, drawn from a random normal distribution with variance 0.1. Although one could expect the solutions for θ = π/4 and θ = 3π/4 to differ only by a sign, this is not the case. There is also a significant difference between the solutions obtained with and without time-step adaptation. These differences can be explained by the sensitivity of the evolution to perturbations of the initial condition and to numerical errors in the early stage of the simulation.
considered the relaxation towards equilibrium of a sessile droplet and the coalescence of two sessile droplets on flat, chemically homogeneous substrates; then we considered several multiphase systems in complex geometries or surrounded by chemically heterogeneous substrates.
Compared to finite differences or spectral approaches, the method introduced here has the advantage that it can be used without modification with complex geometries. Furthermore, the numerical scheme we have proposed can easily be extended to include at least two additional features. First, a linear, energy-stable, second-order scheme could be developed for the threecomponent CH model with wetting boundary conditions, building on the work of [7,8]. Second, we remark that in our work, we considered a regime in which contact line motion is controlled by diffusive interfacial fluxes, or in other words, we considered a large diffusivity limit, where any possible advection effects are neglected. To account for such effects the model must be appropriately modified to include an advection term coupled to the Navier-Stokes equations [1,49,35,36,37,63].
Such generalisations are indeed possible within the proposed numerical scheme and we hope to address these and related issues in future studies.

Acknowledgements
We are grateful to Dr.   . Overall, the time step increases steadily when the adaptive time-stepping scheme is used, which is consistent with the decreasing frequency of coalescence events. The time step is refined at times to ensure that the incremental decrease of free energy at each iteration is approximately constant.
We also recall two other well-known compactness results; see e.g. [41]. Let X, Y, Z be Banach spaces with a compact embedding X ⊂ Y and a continuous embedding Y ⊂ Z. Then the following embeddings are compact: Proof. Without loss of generality, we assume that the mobility, b, is equal to 1. In the spirit of [21, Theorem 2], we apply a Faedo-Galerkin approximation. Let {ϕ n } n∈N and {λ n } n∈N denote the eigenfunctions and eigenvalues of the Laplace operator with a homogeneous Neumann boundary condition, i.e. − ∆ϕ n = λ n ϕ n in Ω, normalized such that Ω ϕ n ϕ m dΩ = δ mn .
We assume without loss of generality that λ 1 = 0. To build an approximation of the solution to Eqs. (15a) and (15b) in the finite-dimensional space S N := span{ϕ 1 , . . . , ϕ N }, we consider the following ansatz, and the variational formulation To this formulation corresponds the following system of ordinary differential equations, with unknown functions a N n N n=1 and b N n N n=1 : for n = 1, . . . , N . Local existence and uniqueness of a solution to this system of equations is guaranteed by the fact that that the right-hand side of (41a) depends continuously on the coefficients a N n N n=1 . To show the existence of a global solution, we will use the a priori estimate presented in the following lemma.
Proof. Settingφ = µ N ,μ = ∂ t φ N in Eqs. (40a) and (40b) and subtracting leads to the equation Using a trace inequality, Hölder's inequality, and Young's inequality with a parameter, we have, for all u ∈ H 1 (Ω), Now we use the simple fact that, for any β > 0 and 0 ≤ s ≤ t, the inequality |x| s ≤ β s + β s−t |x| t holds true for all x ∈ R, to obtain for a constant C independent of u.
By integration by parts of the first term in Eq. (42), we obtain N n=1 λ n (a N n ) 2 < C. This result, together with the inequality implied by Equation (41a) and the fact that λ 1 = 0, show that the coefficients a N n N n=1 do not blow up, and by Eq. (41b) neither do the coefficients b N n N n=1 , implying global existence. In addition to (42), we have the usual estimate on ∂ t φ N : denoting by Π N the L 2 (Ω) projection on S N , for all ψ ∈ L 2 (0, T ; H 1 (Ω)) the following holds: ≤ C ψ L 2 (0,T ;H 1 (Ω)) .