Phase field simulation of liquid filling on grooved surfaces for complete, partial, and pseudo-partial wetting cases

Wetting of solid surfaces by liquids is ubiquitous in nature and critically important for many technological and industrial applications ranging from printing, coating, and microfluidics to oil recovery and carbon capture.^1–4 Given the importance of surface wettability, along with rapid advances in surface engineering techniques, such as lithography, 3D printing, and surface self-assembly,^5–7 understanding the roles of surface topography on the wetting behavior of liquids has emerged as a prominent area of research.

Numerous works to date have investigated how surface structures can give rise to advantageous surface wettability,^8–11 including superhydrophobicity, self-cleaning, drag reduction, and directional spreading. However, the majority of these studies take a macroscopic view of wetting phenomena where the liquid–solid–gas interactions are represented by a single parameter describing the contact angle. At the same time, it is well established in the literature that the intermolecular interactions between the liquid and solid molecules can be highly complex.^1,12,13

Such intermolecular interactions include hydrogen bonds, van der Waals, dipole–dipole interactions, and others.¹⁴ However, their effect on wetting can be understood by looking at a thin liquid film of thickness e on a solid substrate, from which all intermolecular interactions can be incorporated in terms of an effective interface potential. This is defined as the cost of energy per unit area to maintain the thin film at a given thickness.¹³ From this effective interface potential, one can derive the effective force per unit area between the solid–liquid and liquid–gas interfaces, known as the disjoining pressure.¹ 

When considering phenomena at length scales smaller than or comparable to the range of the effective interface potential (typically, of the order of hundreds of nanometers^13,15), a contact angle description of the three-phase interaction in itself is not adequate. It is important to consider the distance-dependent interactions from the solid surface. Indeed, depending on the disjoining pressure profiles (which capture the aforementioned distance-dependent interactions), different wetting states can arise:^14,16 complete wetting, where liquid fully spreads on the solid surface; partial wetting, where a finite contact angle is formed at the three-phase solid–liquid–gas contact line; and pseudo-partial wetting, where the macroscopic liquid domain (i.e., droplet) is surrounded by a thin liquid film.

Previous approaches to computationally study nanoscale fluid phenomena incorporate atomistic details, such as using Molecular Dynamics (MD)¹⁷ and Density Functional Theory (DFT).^18–22 Here, we show how a mesoscale model, the phase field model, can be augmented to enable a wide variety of short- and long-range solid–fluid interactions described above. This is distinct from previous phase field wetting simulations, which typically treat wetting as a boundary condition at a solid surface^23,24 and neglect long-range forces. Since we are interested in static and quasi-static phenomena in this paper, we will directly minimize the free energy of the phase field models.

As phase field models are computationally less demanding than traditional nanoscale methods, the incorporation of long-range interactions should allow highly complicated structures to be studied. This is relevant not only because smaller and more complex features can be reliably manufactured²⁵ but also because they are key for the emergence of interfacial phase transitions, such as liquid adsorption and liquid filling,^20,26–28 which start at the smallest surface features. Such phase transitions are important for many applications, such as thin film condensation and evaporation^29,30 and heat transfer.³¹ 

To demonstrate the versatility of the phase field method, we apply it to study liquid filling and emptying on grooved surfaces as the liquid pressure is varied.^18,32,33 We will compare the results for short-range and long-range liquid–solid interactions. We will also contrast them for complete, partial, and pseudo-partial wetting scenarios. To the best of our knowledge, this is the first systematic liquid filling transition study for the pseudo-partial wetting case. Due to the competition between short-range attraction and long-range repulsion, it leads to several possible pathways and critical pressure dependences on geometry that are distinct from the complete and partial wetting cases.

This paper is structured as follows: In Sec. II, the theoretical basis of the model is introduced. We present our main results in Sec. III, which contains two subsections. The first part of Sec. III is devoted to simulation results on a flat surface, while the second part is for grooved surfaces. We summarize our work and discuss avenues for future work in Sec. IV.

In the simulation, the free energy functional and its derivative are discretized into three types of nodes: bulk fluid nodes, solid nodes, and surface nodes (at the solid boundary). The bulk fluid nodes comprise of a cubic lattice with every adjacent node separated by a lattice spacing, G, with G = 0.02 in simulation units. The interface width for the liquid–gas interface is typically chosen to be ɛ = 2G. The solid and surface nodes are also arranged in a cubic lattice; however, their separation is G/(2g_res − 1), where g_res = 1, 2, 3, …, n denotes the grid resolution. Typically, we use g_res = 2. This grid refinement is useful for increasing the accuracy of calculation of the interaction energy density, $Fs(r)$⁠. Every node in the bulk fluid (or solid) is fully occupied with a unit volume of fluid (or solid). At the surface nodes, however, a node is part fluid and part solid, where the corresponding fractions depend on the local surface structure.

When we compute the numerical integration of interaction energy density $Fs(r)$ as in Eq. (5), all effective liquid–solid interactions are taken into account. These include contributions from bulk fluid–solid, bulk fluid–surface, surface–solid, and surface–surface node interactions. The detailed scheme is provided in the supplementary material (Sec. S1). In addition, we also employ several periodic images of the solid domain to ensure the long-range liquid–solid interactions are sufficiently accounted for. Although the interaction energy density calculation is computationally expensive, particularly for a large domain and a high solid node resolution, it is only calculated once at the start of the simulation.

Upon the energy minimization routine, the discretized order parameter in the bulk fluid and surface nodes will evolve toward the minimum energy configurations. We employ the L-BFGS algorithm due to its efficiency for problems with a large number of degrees of freedom. For details on the energy minimization routine, see Refs. 34 and 35.

In our model, $γlg=8/9$⁠, and γ_ls can be calculated when the film thickness is sufficiently thick such that γ_ls = F(e → ∞) − γ_lg. P(e) can then be determined by evaluating F(e) at varying film thickness values. To get the variation of P(e), we simulate a liquid film with small interfacial area on a flat surface to avoid the coexistence between a liquid film and a dry solid. For convenience, here we employ the volume constraint as given in Eq. (13) and vary the film thickness by adjusting the target volume V₀ in the simulation. Different variations of P(e) representative of different wetting states are shown in Fig. 1 (for the disjoining pressure Π(e) profiles of these wetting states, see the supplementary material, Fig. S2). These capture the profiles previously proposed in the literature, such as those presented in the work of Brochard-Wyart et al.¹⁴ In the insets, we show an equilibrium state of a sessile droplet placed on a flat surface under the respective wetting states.

Figures 1(a) and 1(b) show the complete wetting case indicated by the positive value of P(0) (S > 0) and the formation of a liquid thin film (insets). In panel (a), we use a large negative α in the effective interaction $ulr1$⁠. Here, the functional P(e) decreases as the film thickness increases, which results in positive Π(e) and S. In panel (b), the value of α is less negative than that in panel (a). This makes P(e) non-monotonic, leading to negative Π(e) at small e. The resulting value of S is also smaller but remains positive. A similar profile of P(e) can be seen when using the effective interaction $ulr2$⁠. Here, |−u_o| ≫ |u_p| is used to give a strong attractive interaction near the solid surface to allow the droplet to spread across the surface.

Figures 1(c) and 1(d) show the partial wetting case indicated by the negative value of P(0) (S < 0), leading to droplets with finite contact angles [θ < 90° for panel (c) and θ > 90° for panel (d)], as depicted in the inset of the figures. In panel (c), we still use a negative α in effective interaction $ulr1$ but the value is smaller than that in the complete wetting case. As a result, the liquid–solid interaction is weaker. Similar to Fig. 1(b), P(e) is increasing at small e but decreasing at large e. However, the resulting spreading parameter S is negative. If we now switch to positive α, P(e) is monotonically increasing with a larger negative S, as shown in panel (d). For α = 0, θ = 90° and S < 0. We can also show the partial wetting case using the effective interaction $ulr2$⁠, as depicted in Figs. 1(c) and 1(d), where u_p is increased to make a stronger repulsive interaction preventing the droplet from completely spreading.

If we tune the variable parameters in $ulr2$ such that the short-range attractive interaction is strong enough to allow the liquid spreading and the long-range repulsive interaction is sufficient to stabilize a droplet, we will obtain a pseudo-partial wetting case, where a droplet is surrounded by a thin liquid film wetting the solid surface, as shown in the insets of Figs. 1(e) and 1(f). The spreading parameter S can be negative or positive,^14,37 and the P(e) profile is characterized by a minimum at a certain e. In panel (e), the attractive term is quite strong at short ranges (due to a large |−u_o|) that P(e) is decreasing. Since S > 0, the thin film extends indefinitely. At long ranges, the repulsive term becomes more dominant (due to moderate value of |u_p|) to change the direction of P(e) to be increasing. The droplet formed in this condition has a lower contact angle as seen in the inset. When the strength of the repulsive terms is increased but the attractive term is kept unchanged, the decreasing trend of P(e) at short ranges reduces and the increasing trend of P(e) at long ranges increases, as shown in panel (f). As can be seen, the contact angle of the droplet is larger (see inset). Moreover, since S < 0, the thin film does not extend indefinitely, as illustrated in the inset. The pseudo-partial wetting case cannot be obtained with effective interaction $ulr1$⁠.

When considering large-scale wetting phenomena, the long-range liquid–solid interactions discussed in the preceding subsection are often not directly relevant as they occur at much smaller length scales. The short-range surface energy density implemented in the free energy is directly related to contact angle at the surface via Eq. (11). Similar short-range energy densities have previously been demonstrated in the literature.³⁸ Here, the main difference is the quartic form of f(ϕ_s). In Fig. 2(a), we compare the measured contact angle of a sessile 2D drop from the simulation, labeled θ_sim, with the input contact angle, θ_theory. To measure the contact angle, we fit a circular arc to the drop profile. We found an excellent accuracy of the contact angle with the error of <1° [Fig. 2(b)]. Such accuracy is superior compared to a range of frequently used forms of f(ϕ_s)³⁹ including linear and cubic models. The comparisons between the different forms of f(ϕ_s) are provided in the supplementary material, Table S1 and Fig. S3.

Our next investigation is the wetting behavior of liquid on a structured surface. In this context, we consider a long periodic grooved surface with groove width of d, depth h, and wall barrier width w as shown in Fig. 3(a). To reduce the simulation cost, it is only necessary to simulate a single groove unit cell with periodic boundary condition being applied in the x and y directions to capture the periodicity of the grooves. In this work, the simulation domain size is chosen to be N_x = w + d, N_y = 6G, and N_z = h + 50G. The groove dimension is taken as w = h = 50G and d = 20G, unless stated otherwise. The typical liquid–solid interaction energy densities due to long-range interactions across the system are shown in Fig. 3(b) for the effective interactions $ulr1$ and $ulr2$⁠. Here, the energy density is scaled by γ_lg/G. For $ulr1$⁠, the parameter α is taken as a negative value, hence liquid and solid experience an attractive interaction that is higher at the surface and decays toward zero farther from the surface. This is depicted in Fig. 3(c) for complete and partial wetting cases at x = 0 and varying z. The decay rate of $ulr1$ depends on the parameter β. The higher β, the slower the decay and the longer the interaction tail. For $ulr2$⁠, the attractive interaction only occurs near the surface, and the interaction becomes repulsive in the bulk of the liquid, as depicted in Fig. 3(c) for the pseudo-partial wetting case. It is also worth noting that for both effective interactions, the liquid–solid interaction is stronger at the bottom corners and weaker at the top corners of the barrier wall, consistent with observations from MD⁴⁰ and DFT^19,20 simulations.

Figures 4(a)–4(c) show the filling and emptying transitions as the pressure Δp is varied, for the case of complete wetting. Here, we use $ulr1$ for the long-range interaction with α = −8 × 10⁻⁴ and β = 0.5. Similar results are obtained when $ulr2$ is used. Upon increasing the liquid pressure, liquid begins to fill the grooved surface, as shown in Fig. 4(a). We can categorize the liquid filling process into three stages,⁴¹ viz., (i) pre-filling, (ii) capillary filling, and (iii) post-filling, which occur after one and the other with increasing liquid pressure.

The pre-filling stage occurs at large negative Δp, which means the pressure in the liquid is much lower than that in the gas phase. Here, the liquid forms a thin film that follows the shape of the groove structure. The thickness of the film depends on the strength of the interaction (parameter α in $ulr1$⁠) and increases as the liquid pressure is increased. The dependency is well approximated by e ∝ (2α/Δp)^−1/3 at the bottom, top, and sides of the barrier wall, where α can be associated with the Hamaker constant.^42,43 Assuming the grooves have the dimension of order of a hundred nm, the values of α in our simulation translate 10⁻²⁰ to 10⁻¹⁹ J, which are the typical values of the Hamaker constant.¹² 

The key to understanding the capillary filling stage lies in changes to the film thickness in the bottom corners of the groove. As the menisci in the corners grow and approach d/2 in size, liquid from either side merges and rapidly fills the gap. The liquid interface then rises up from the bottom of the groove. This can be seen from the sudden increase in liquid film thickness as calculated at the middle gap, $ẽg$ [Fig. 4(b)] and from the sudden decrease of the liquid–gas interfacial area, $Ãg$ [Fig. 4(c)]. We define the critical pressure Δp_c as the pressure value with the largest gradient in the $ẽg$ and $Ãg$ plots. At this capillary filling stage, the growth of thin film at the side walls and top of the barrier wall still follows e ∝ (1/Δp)^−1/3.

The liquid, however, does not immediately fill the whole gap of the groove. In the post-filling stage, with increasing pressure, the liquid–gas interface between the barrier walls starts to smooth out until it becomes flat. The thin film thickness at the top of the barrier wall also increases more rapidly compared to the pre-filling and capillary filling stages. This occurs when Δp has small negative values, which means p_l → p_g. As p_l > p_g, the film thickness increases to infinity as liquid fills up the whole domain.

If Δp is reversed from positive to large negative values, the liquid will be emptied from the grooved surface. Upon decreasing Δp, the liquid–gas interface follows the reverse path as the liquid fills the groove surface [Fig. 4(c)]. Therefore, the liquid filling does not exhibit hysteresis behavior for the complete wetting case. Recently, similar filling transitions have been investigated via DFT.²⁰ It was also observed that the filling transition is mediated by the growth of the menisci in the bottom corners of the groove. However, here, we are also able to show the contribution of the films on the sidewalls and top of the barrier wall.

Next, we want to compare the effect of long-range and short-range liquid–solid interactions on the liquid filling transition. In this case, we use the effective interaction $ulr1$ in the interaction energy density $Fs(r)$ for the former and $Fs(rs)$ as in Eq. (10) for the latter. Although the filling behavior for both interactions is qualitatively the same, the filling transition occurs at different critical pressures [Fig. 4(b)]. This is because the liquid film thickness at the wall is different, which changes the effective separation between the walls. The critical pressure dependency can be inferred from the Kelvin equation, in which Δp_c is expected to be inversely proportional to the effective wall separation.

The liquid film formed due to long-range interactions is thicker than that due to short-range interactions. For the former, the contribution of the liquid–solid interaction is determined by how long the tail of the decaying interaction is until it becomes essentially zero. This is controlled by the parameter β. The higher the value of β, the longer the tail. As a result, the interaction with higher β forms a thicker liquid film at the wall. For the short-range interaction, the liquid–solid interaction is assumed to occur only at the surface of the solid. Therefore, there is no liquid–solid interaction contribution farther from the surface. Hence, the liquid film is thinner, and the filling transition occurs for larger Δp.

We now turn our attention to the partial wetting case (θ > 0°). The results presented here employ the short-range interaction. Equivalent results are obtained for the long-range interactions once the contact angles are mapped. The liquid filling behavior for the partial wetting case is illustrated in Figs. 4(d)–4(f). In the same manner as the complete wetting case, we can also group the filling process into (i) pre-filling, (ii) capillary filling, and (iii) post-filling stages. The capillary filling stage is also marked by a critical pressure at Δp_c. Here, we have to divide our discussion into two scenarios,¹⁸ which are for θ < 45° and for θ > 45°.

For θ < 45°, in the pre-filling stage, liquid condensation could be nucleated at the corner of the groove forming menisci at a large negative Δp [Fig. 4(d)]. In this case, the corner menisci grow as the liquid pressure increases until they merge as a single meniscus. Once this has happened, we enter the capillary filling stage, in which the liquid starts filling the gap while maintaining the shape of the meniscus. In sharp interface models, at a certain Δp_c, the filling transition occurs as signified by an abrupt increase in $ẽg$⁠. Δp_c can be predicted using the Kelvin equation, as will be discussed in Sec. III B 4. Using the diffuse interface model in the present study results in a rounding of this first-order phase transition. However, as is shown in the supplementary material (Sec. S4), this effect is marginal if there is a suitable separation of length scales (at least a factor of 10) between the diffuse interface width and the wall height. Thus, as is shown in Fig. 4(e), the partial wetting filling transition is still sharp compared to the complete wetting case. After the filling transition occurs, the liquid again does not completely fill the gap, as in the complete wetting case, but both ends of the liquid–gas interface are pinned in the top edge of the wall. In the post-filling stage, the meniscus starts to flatten as p_l → p_g. When Δp turns positive, the curvature of the meniscus also turns sign from negative to positive. As the liquid manages to overcome the contact line pinning, it fills up all the gas phase.

For θ > 45°, the pre-filling stage is marked by a gas-like phase with a completely dry solid [Fig. 4(e)]. The corner menisci do not form, and the filling transition immediately occurs when the pressure has reached Δp_c. The liquid will then be pinned at the top edge of the groove with smaller curvature due to higher θ. The post-filling stage is then similar to that for θ < 45° except that the positive curvature at positive Δp could grow larger in size before it overcomes the contact line pinning and fills all of the gas phase. This means that Δp at which the liquid fills the gas phase occurs at a larger value than that for θ < 45°.

Figure 4(f) shows liquid filling and emptying paths for increasing and decreasing Δp. The hysteresis behavior is clearly pronounced. Due to contact line pinning at the top corner of the walls, during the filling process, the meniscus curvature changes from negative to positive as Δp increases. During the emptying process, however, upon decreasing Δp the liquid continues to wet all the surface and maintains a flat liquid–gas interface until a significantly lower pressure difference. Once the top of the wall is fully dewetted, the liquid gets pinned at the top corners with negative curvature. The liquid emptying path then follows along the same path as the liquid filling (see the supplementary material, Fig. S5, for snapshots of configurations during the filling and emptying processes). This hysteresis behavior in the partial wetting case has also been reported elsewhere.^32,44

Figures 4(g)–4(i) show the liquid filling behavior on a grooved surface for the pseudo-partial wetting case. The effective interaction $ulr2$ is used in the interaction energy density $Fs(r)$⁠. The magnitude of parameter −u_o controls the strength of the attractive interaction. To obtain a pseudo-partial wetting state, a large enough −u_o is employed to get a liquid film near the surface. The parameter σ controls the thickness of the liquid film. Here, we use σ = 3, u_o = −5 × 10⁻⁷, and u_p = 2 × 10⁻⁷.

In the pre-filling stage (at large negative Δp), in contrast to the partial wetting case at the same Δp, the liquid wets the bottom surface of the groove and the side walls, forming a liquid film, but it leaves the top of the barrier wall dry as the liquid film is pinned at the top edges of the wall. At the bottom corners of the groove, the meniscus of liquid condensation is not as pronounced as it is for the full and partial wetting case. This is because the interaction at the surface near the corner slightly reduces due to the effect of the repulsive term in the effective interaction $ulr2$⁠. As the liquid pressure increases, the liquid overcomes the contact line pinning at the top edges and covers the top of the wall. This is shown in Fig. 4(i) (right-pointing triangle), in which the liquid–gas interface area increases abruptly at Δp/(γ_lg/hG) ≈ 0.3. The bottom corner menisci only slightly grow with increasing pressure, unlike for the complete and partial wetting cases where they grow and merge as their size approaches d/2.

In the capillary filling stage, the critical pressure for the filling transition occurs sharply at positive Δp_c [indicated in Fig. 4(h)]. The sharp transition applies for narrow and wide groove widths. Compared to the negative Δp_c observed for complete and partial wetting cases, this suggests the filling transition is more energetically expensive for the pseudo-partial wetting case. At the critical pressure Δp_c, liquid fills the gap, and it forms a droplet in the middle of the gap coexisting with the liquid film on top of the barrier wall [Fig. 4(g)]. Such coexistence is reminiscent of the morphology observed on a flat surface. However, the range of stability of the droplet is limited. With increasing pressure in the post-filling stage (with positive Δp), the droplet becomes unstable and the liquid fills the simulation domain.

The hysteresis behavior is also clearly observed in the pseudo-partial wetting case, as shown in Fig. 4(i) (see the supplementary material, Fig. S5, for snapshots of configurations during the filling and emptying processes). During the filling process, the contact line pinning at the top corners of the walls allows the liquid to form a droplet bulge in the middle of the gap as Δp increases. Upon decreasing Δp, however, the droplet bulge slowly flattens until the liquid filling the gap abruptly drains, leaving a liquid film that follows the shape of the groove structure. The top of the wall remains covered by a liquid film, hence we find higher $Ãlg$ than in the liquid filling path.

The pseudo-partial wetting case, however, cannot be captured by a relation akin to Eq. (18) or Eq. (19). We argue that this is because the corner menisci are not so apparent during the capillary filling stage and do not merge into a single meniscus before the filling transition occurs. Therefore, d does not affect Δp_c. This can be observed for d/h > 0.6 in Fig. 5. In this scenario, the liquid fills up the simulation when the filling transition occurs, as illustrated in Fig. 6(v). When d is very small (d/h < 0.6), however, Δp_c starts to be dependent on d, but it still does not obey Eq. (18). Instead, we find this variation is accompanied by nontrivial changes in the morphological pathway during the filling transition. With increasing groove width, five distinct pathways are identified, illustrated in Fig. 6, and indicated in Fig. 5: (i) liquid fills the gap forming a liquid–gas interface with a negative curvature while keeping the top of the wall dry, (ii) the same as scenario (i) except that top of the wall is covered by a liquid film, (iii) the same as scenario (ii) but the liquid film at the top of the wall is formed before the capillary filling stage, (iv) the same as scenario (iii) but the liquid–gas interface curvature is positive, and (v) liquid fills the system at the critical pressure.

In summary, we have presented systematic numerical studies of liquid filling on grooved surfaces using a phase field method. We consider both short-range and long-range liquid–solid interactions. The latter include purely repulsive and attractive interactions, and more complex interactions with short-range attraction and long-range repulsion. To the best of our knowledge, such versatility allows us to capture complex disjoining pressure profiles for the first time in a phase field approach, in agreement with previous studies using atomistic modeling^17,22 and analytical theory,^14,37 which in turn give rise to complete, partial, and pseudo-partial wetting states. In this work, we have also introduced a quartic polynomial to switch the interaction energy density between the liquid and gas phases [f(ϕ) in Eq. (8)]. This polynomial prevents enrichment of the liquid and/or gas phases on the solid surface, and it leads to more accurate contact angle calculations compared to the linear and cubic forms previously used in the literature.

We rationalize the liquid filling process on grooved surfaces into three stages: (i) pre-filling, corresponding to the growth of the thin film around the structure (for complete or pseudo-partial wetting), or liquid menisci in the bottom corners of the groove (for partial wetting); (ii) capillary filling, where there is a rapid increase of liquid volume in the groove marked by a critical pressure; and (iii) post-filling, typically signified by flattening of the liquid–gas interface before liquid completely fills the whole domain. Comparing the results for complete, partial, and pseudo-partial wetting, we find there is no contact line pinning for the complete wetting case and the liquid filling and emptying trajectories are reversible. In contrast, we observe clear hysteretic behavior for partial and pseudo-partial wetting, caused by the coexistence of two metastable states over a pressure range. In the partial wetting case, late in the filling transition, pinning of the interface on the top corner of the wall leads to a state that remains metastable over a range of positive pressures. Coexisting with this is the unpinned state, which at positive pressures sees liquid completely fill the system. In the pseudo-partial wetting case, the origin of metastability is different. Here, repulsive interactions in the center of the groove energetically penalize partial filling of the groove. Instead, either the groove remains almost empty or the groove is full. Considering the critical pressure, although the diffuse interface marginally rounds the first-order transitions, both the complete and partial wetting cases follow a Kelvin-like equation for its dependence on the groove width: large and negative at small widths and plateaus to zero for large widths. The pseudo-partial wetting case, however, is different. At large widths, the critical pressure is positive and constant. We find the critical pressure does depend on the groove width for smaller widths, and interestingly, this is accompanied by morphological changes in the trajectories of the liquid filling process.

There are a number of exciting avenues for future work. Here, we have considered grooved surfaces. It is straightforward to extend the study to other, more complex surface geometries, such as reentrant geometries, seesaws, hierarchical posts, or even nonsymmetric structures, which have extensively been harnessed in wetting applications. Another possible direction is to consider the liquid dynamics, beyond the quasi-static results presented here. There are some limitations, however, of the phase field method used in this study. The present method does not include the interfacial fluctuation effects, which have been shown to occur at nanoscale⁴⁶ and captured by atomistic simulations.⁴⁷ To capture these phenomena, one possible route is to couple the phase field model here with fluctuating hydrodynamics methods.^48,49 Finally, we hope the simulation results will inspire experimental studies to verify our theoretical predictions.

The supplementary material contains (i) our scheme for the calculation of Fs(r); (ii) the corresponding disjoining pressure profiles for the effective interfacial potential shown in Fig. 1; (iii) comparison of contact angle results for different forms of f(ϕs); (iv) a short discussion on the rounding of first-order phase transition due to the diffuse interface; and (v) liquid–gas interface configurations during the filling and emptying processes for partial and pseudo-partial wetting.

F.O. acknowledges a BPPLN scholarship from the Directorate General of Resources for Science Technology and Higher Education, Republic of Indonesia. H.K. and J.R.P. would like to acknowledge UKRI Engineering and Physical Sciences Research Council (Grant No. EP/V034154/1) for funding.

The authors have no conflicts to disclose.

Fandi Oktasendra: Conceptualization (equal); Investigation (lead); Methodology (equal); Writing – original draft (equal). Arben Jusufi: Conceptualization (equal); Funding acquisition (lead); Writing – review & editing (equal). Andrew R. Konicek: Conceptualization (equal); Funding acquisition (equal); Writing – review & editing (equal). Mohsen S. Yeganeh: Conceptualization (equal); Funding acquisition (equal); Writing – review & editing (equal). Jack R. Panter: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal). Halim Kusumaatmaja: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.