Modelling of water droplet dynamics on hydrophobic soils: a review

. The hydrophobicity of soils (or soil water repellency) can be naturally promoted by wildfires or synthetically induced by hydrophobic compounds (polydimethylsiloxane, tong oil, etc.). Soil phenomena can be related to hydrophobicity, such as soil erosion (splash erosion and rill erosion) and post-wildfire debris flows. The hydrophobicity of soils is characterized by the contact angle, and the interactions between water droplet and solid particles including spreading, oscillation, and infiltration. Early studies on soil water repellency mainly focus on the experimental aspects, while with the development of advanced numerical tools, numerical methods have been widely applied to study the hydraulic properties of hydrophobic granular materials in recent years. This paper comprehensively investigates the different numerical methods for modelling the interaction between water droplets and hydrophobic soils, i.e., smoothed particle hydrodynamics (SPH), lattice Boltzmann method (LBM), material point method (MPM), and volume of fluid (VOF). The features of different method are summarized, and the future work are discussed.


Introduction
Soil hydrophobicity is defined as a reduction of affinity between soil particles and water [1]. The hydrophobic properties of soils were first investigated in Californian soils [2]. At that time, hydrophobic soils were described as "could not be wetted, either by man, by rain, irrigation, or the movement of water from the subsoil". During the 1990s, soil hydrophobicity has been extensively reported [3], arousing people's attention to this problem. Persistent findings on this phenomenon indicate that hydrophobic soils do not occur in isolation but are widespread all over the world [4][5][6][7][8].
Early studies focused on the naturally generated hydrophobic soils [9][10]. The hydrophobicity can be formed due to the presence of organic coating related to the fungi, oily materials, or any waxy materials generally produced by plants. It is also reported that there is a strong link between wildfire and soil hydrophobicity [11,12], owing to the generation of oily substances after burning the plant. Soil hydrophobicity was generally considered as negative in slope processes, since it is regarded as an aggravating factor for some natural disaster such as soil erosion [13,14] and debris flows [15,16]. However, in recent years, researchers began to consider taking advantage of soil hydrophobicity to benefit mankind. Therefore, some work based on manually prepared hydrophobic soil has been carried out. For example, some researchers used synthetic compounds (polydimethylsiloxane, tong oil, etc.) to make the soil hydrophobic, and then the hydraulic and mechanical behaviours of hydrophobic soils were investigated [17,18]. It is expected that these hydrophobic soils can be used in ground infrastructure * Corresponding author: lourenco@hku.hk to delay or restrict water infiltration and thus prevent their failure [19]. Hence, soil hydrophobicity is crucial, whether from the perspective of avoiding its disadvantages or utilizing its merits.
Contact angle (CA) is commonly used to quantify the hydrophobicity. The CA at a granular surface can be describe at two different length scales, namely, intrinsic CA and apparent CA. The intrinsic CA is defined as the intersection angle between the water-air interface and the base plane, as shown in Fig. 1. The intrinsic CA ( 0 ) is determined by the Young's equation [20], which is defined as follows: where the , , and are three interfacial surface tension coefficients between solid-gas, solid-liquid, and liquid-gas. When 0 < 90°, the surface is considered to be hydrophilic or wettable; otherwise, the surface is hydrophobic or water-repellent.
The intrinsic CA at a granular surface is the CA on a single particle, while the CA on the surface consisting of many particles in called apparent CA (as shown in Fig. 2). The apparent CA can be measure in the laboratory based on the sessile drop method (SDM) [21,22]. If a droplet is placed on a hydrophobic granular surface, spreading and oscillation can be observed at the initial stage, and then the droplet may reach a stable state under the action of gravity, surface tension, and its interaction with solid particles [23]. These interaction behaviours between fluid and solid, i.e., spreading, oscillation and infiltration, are also significant for charactering the hydrophobicity of granular materials [23][24][25]. Early studies on the interaction between water droplet and hydrophobic surfaces remained at an experimental level [26][27][28]. In the past two decades, with the development of advanced numerical tools, numerical methods were applied to study the dynamic behaviour of droplets, such as the mesh-based method VOF [29], the particle-based method SPH [30,31], and other methods like the LBM [32] and MPM [33].
This paper aims to review the current state-of-the-art of droplet dynamics on hydrophobic surfaces based on numerical methods. First, the mathematical theories and algorithms that generally used in the simulation are summarized in section 2; then the research status and results from numerical modelling are described in section 3; discussions, conclusions, and future work are in section 4.

Theories and algorithms
Although there are several different numerical methods, they are basically based on the same theories and algorithms. In this section, the theories and algorithms that are commonly used for modelling droplet dynamics on hydrophobic surfaces are presented.

Governing equations
The movement of fluid phase (liquid and gas) are controlled by the Navier-Stokes (NS) equations. The NS equations in Lagrangian description are defined as follows: where is the density, is the velocity, is the dynamic viscosity. Generally, represents gravity in the single-phase flow, while also include surface tension forces in the multiphase flow. For different numerical methods, they use different means to discrete the governing equations (2)-(3).

Surface tension
The continuum surface force (CSF) model, which was developed by Brackbill et al. [34] in 1992, is widely used to calculate the surface tension in multiphase flows. The surface tension is modelled as a volume force in the CSF model, and a colour function c(x) is used to distinguish two fluids. As shown in Fig. 3, different fluids (fluid-1and fluid-2) are represented by different values (c1 and c2) in the colour function. The interface between fluids should be sharp so that the colour function is discontinuous. However, in CSF models, to calculate the surface tension of each element (mesh or particle) numerically, a transition region is used to bridge the discontinuity at the interface (as shown in Fig.  3). Therefore, the values of the colour function c(x) can be calculated continuously in the transition region. The computational domain can be discretized by grids or particles, depending on the different numerical methods used. The surface tension sv can be calculated based on the following equation: where is the coefficient of surface tension, ∇ ( ) is the gradient of colour function, and [ ] is the difference between the colour function values of two fluids ([ ] = | 1 − 2 |). is the curvature at the interface, which is the divergence of the normal direction n [34], as shown in the following equation: Based on the original CSF model [34], several new methods have been developed to calculate surface tension more accurately or in complex conditions. For example, Adami and Hu et. al. [35] reproduced a new divergence approximation and proposed a new formulation for the surface curvature based on the CSF model [34], which has advantages in the computation of problems with large density ratio and complex interfaces.
Another method generally used in particle-based method treats the surface tension as a particle-particle interaction force [36]. This method was original developed based on SPH method. When | − | ≤ ℎ, the inter-particle interaction force between particle-i and particle-j is given by the following equation [36]; where is the strength of the force acting between particle-i and particle-j, h is the support domain in SPH method and " − " is a vector pointing from particle-i to particle-j. Then, the item is introduced to the motion equation for including surface tensions. It should be noted that is a parameter related to surface tension, and the value of surface tension coefficient should be calibrated by numerical experiments [36].
Some methods solving the surface tension implicitly had also been developed. Hysing [37] used a semiimplicit strategy for calculating surface tension based on a variational CSF. Schroeder et. al. [38] treated the surface tension as a semi-implicit Lagrangian force, and Zheng et. al. [39] further developed their method and used an implicit technique for surface tension in the hybrid particle/grid-based method. Other numerical models for computing surface tensions can be found in the review paper [40].

Contact angle
Considering the surface roughness of granular materials, the apparent contact angle can be divided into two types: Wenzel's contact angle ( ) [41] and Cassie-Baxter contact angle ( ) [42]. The Wenzel state of droplet forms when the fluid fills surface indentations (as shown in Fig. 4a). The is defined as follows: where r represents the roughness ratio and can be calculated based on the true and apparent surface area of a granular surface (r = true surface area / apparent surface area). If a gas phase is present between the solid and liquid and fills the troughs on the granular surface, the Cassie-Baxter state is formed (as shown in Fig. 4b). The is defined as follows: where 1 and 2 are the fractional areas of the solid in contact with the liquid and air ( 1 + 2 = 1). When the troughs and depressions of a rough surface are partially filled by the liquid, the state CA is a combination of Wenzel and Cassie-Baxter models. Details can be found in the literature [19,43].
Several methods have been proposed to numerically model the contact angle with a solid boundary. The first one is by imposing a static CA at the wall [44]. In this way, the static CA and surface tension coefficient between liquid and gas are input parameters, and the current CA can be estimated at each calculation time step. If the current CA doesn't equal to the pre-set static CA, an extra force is applied to make it deform towards forming the input CA. In this method, the surface tension is only calculated at the liquid-gas interface.
This method is widely applied in CSF models and has been used by different numerical methods, such as VOF [45], SPH [44], LBM [23], MPM [46], and others.
Based on the study of Tartakovsk and Meakin [36], Shigorina et. al. [47] added to equation (6) to include the contact angle implicitly. The contact angle can be measured during or after the simulation. Chen et. al. [48] assumed that the surface tension coefficient between solid-gas is zero with a free-surface assumption, then the solid-liquid contact angle can be estimated by − ⁄ according to equation (1). Different contact angles can be produced by using different surface tension coefficients between solid-liquid and liquid-gas interfaces.

Numerical modelling results
In this section, recent developments on modelling droplet dynamics on hydrophobic soils using different numerical methods are summarized. It should be mentioned that very few studies are directly aimed at modelling soil surfaces, while a substantial body of work concentrates on flat/inclined smoothed surface, and rough surfaces consisting of microstructures or idealized particles such as spheres. Therefore, the modelling studies on liquid drops on smooth/rough surfaces are also present here, as an early simplification and to provide a basis to model soil surfaces.

Volume of fluid
VOF [29] is a widely used grid-based method to model multiphase flow. The VOF method considers the amount of gas and liquid in each grid based on a volume fraction function. In this method, incompatible fluid components share a set of momentum equations, where the volume fraction is used to track the interface between different phases in the calculation domain. CSF models are commonly used in the VOF method to compute the surface tension [49,50].
The VOF method has been applied to several fields concerned with the modelling of droplet dynamics. On smoothed hydrophobic surfaces, Šikalo et. al. [50] studied the droplet spreading and dynamic contact angle, where surfaces with different wettability were used to investigate the effect of static contact angle on the dynamic behaviour of droplets, such as dynamic contact angle and drop spreading diameter. Roisman et. at. [49] simulated the spreading of drops impacting a dry solid surface at low Weber numbers, with the predicted parameters (spread factor, height of the drop and apparent dynamic contact angle) validated against experiments. Gunjal et. al [51] modelled the dynamic behaviour including spreading, rebounding, splashing, and bouncing of the droplet on a flat solid surface. Bussman et. al. [52] investigated the droplet impacts on a 45° inclined surface based on the VOF method. Their numerical results were compared with experimental data, and based on their study, a model in which the contact angel is a function of contact line velocity was proposed.
Besides studies of droplets on smoothed surfaces, the VOF method has also been employed to simulate the droplet dynamics on a porous medium. For example, Das et. al. [45] placed a liquid droplet on a surface with a porous structure and investigated the droplet spreading behaviour. The effect of the porosity and the equilibrium contact angle on the spreading behaviour, liquid imbibition, and apparent contact angle were studied. Saha and Mitra [53] simulated the capillary flow in microchannels with integrated pillars, and the effect of dynamic contact angle on capillary filling was investigated.

Smoothed particle hydrodynamics
SPH is a particle-based method, where fluid flow is described as a group of interacting particles carrying various physical quantities, including mass and velocity, serving as the integration point for the governing equations. SPH has been proven to be a powerful tool in modelling fluid, solid, and fluid-solid interactions [54]. It is also able to simulate small scale processes and parameters such as surface tension and contact angles [55].
Modelling of droplet dynamics on hydrophobic smooth surfaces based on SPH has already been accomplished by various authors. Huber et. al. [56] simulated the wetting effect and analysed the advancing and receding contact angles based on CSF model. Breinlinger et. al. [44] studied the evolution of liquid drops with different intrinsic contact angles, and the pin effect of droplets on slopes with different angles was also investigated. Yeganehdoust et. al. [57] simulated the equilibrium contact angle and found that the droplet shape and equilibrium contact angle is influenced by the dimensionless Bond number. Furthermore, the droplet collision over an inclined surface with different inclined angles and Bond numbers have also been investigated. Tartakovsky and Meakin [36] modelled contact angles and droplet flow through a Y-shaped fracture junction. In their study, the surface tension was calculated based on pairwise fluid-fluid and fluid-solid particle-particle interactions rather that the CSF model.
Shigorina et. al. [47] further developed Tartakovsky and Meakin's [36] method and applied it to model droplets interactions with rough surfaces. They designed surfaces with different roughness (as shown in Fig. 5) and discussed the effect of surface roughness on the apparent contact and states of droplet (Wenzel state, Cassie-Baxter state). Meng et. al [58] studied the droplet spreading on porous substrates based on SPH, and two spreading stages (power-law inertial stage and viscous stage) were found, which were used to determine the evolution of the wetting radius with time. Sivanesapillai et. al. [59] carried out pore-scale simulations for multiphase flow in a porous medium based on the CSF model. They showed that their method is capable of modelling saturation-controlled primary drainage and main imbibition at heterogeneous pore spaces.

Lattice Boltzmann method
LBM is suitable for mesoscale and multiphase simulations. LBM is widely used to model the behaviour of a droplet on a flat or inclined smooth surface [60][61][62]. As for modelling water droplets on rough surfaces, Tang et. at. [63] studied the contact angle and wetting state with an array of microstructures as the substrate. A correlation between the microcolumn height and the contact angle was found in their work. Specifically, it was found that, with increasing microcolumn height, droplets transformed from the Wenzel to the mixed wetting state when = 105°− 125° ( : apparent contact under the Wenzel state), while droplets switched from the mixed wetting state to the Cassie state for = 125°− 140°. Kang et. al. [23] investigated the droplet dynamics on an idealized granular surface. The effects of solid particle size on the apparent contact and dynamic behaviour including oscillation, spreading and infiltration were investigated.
Furthermore, Suo et. al. [64] studied the behaviour of droplets on porous media, focusing on the spreading and imbibition process, and the curvature effect of the porous surfaces. Wang et. al. [65] proposed a theoretical model in which the apparent contact angle is a function of time. Then they introduced this model into LBM to simulate the multiphase flow in porous media, it was observed that the relative speed of the wetting transition was influenced by a dimensionless time ratio.

Material point method
MPM combines Lagrangian material particles (points) with Eulerian Cartesian grids [66]. The use of MPM in modelling droplet dynamics on hydrophobic surfaces is limited, especially in engineering applications. Chen et. al. [46] introduced the CSF model into the generalized interpolation MPM and modelled the capillary rise with an equilibrium contact angle. In computer graphics, several studies simulated contact angles based on the MPM. Hyde et. al. [67] developed an implicit updated Lagrangian formulation for simulating liquids with large surface tension. They showed that their method is capable of modelling the bouncing of a liquid drop on a hydrophobic surface. Based on the work of Hyde et. al. [67], Chen et. al. [48] modelled the surface tension with a momentum conserving implicit MPM. In their method, two different surface tension coefficients were used for a liquid drop to distinguish the liquid-gas and liquidsolid interfaces of a droplet on a surface. Thereafter, the contact angle can be implicitly defined by controlling the ratio of the two coefficients. Their method is able to model a droplet on a hydrophobic smoothed surface with different contact angles.

Concluding remarks
This review introduces a variety of numerical methods that can simulate droplet dynamics on hydrophobic surfaces. However, each method has different characteristics which may impair or facilitate their application on hydrophobic granular surfaces.  [60][61][62], rough surfaces [23,63], and porous media [64,65] MPM Combines the features of mesh methods and particle methods Validation is needed Droplet dynamics with smoothed surfaces [48,67] As a mesh-based method, one notable aspect about VOF is its high computational accuracy. However, the accuracy of VOF is dependent on the quality of the mesh, and its use on complex interfaces is inconvenient. In SPH, the mass conservation is automatic satisfied since the mass is carried by particles. Meanwhile, the calculation accuracy of SPH method is less dependent on particle arrangement, which is an advantage compared with VOF. Moreover, handling complex interface shapes including surface splitting and merging is straightforward in SPH. One disadvantage of SPH is that the boundary generally needs special treatment. LBM is suitable for mesoscale modelling and have an advantage in simulating complex materials such as porous media, while the computational efficiency in LBM is a challenging problem. MPM is considered to have the advantages of both a grid method and particle method to a certain extent [66]. For example, it can be used to simulate self-collision and fracture based on Cartesian grids and simulate splitting and merging behaviours based on particles [66]. However, the application of MPM in modelling droplet dynamics on hydrophobic surfaces is limited, and the results also need to be further validated. The pros and cons of different numerical methods, and the application of these methods in modelling droplet dynamics on hydrophobic surfaces, are summarized in Table 1.
Previous research mainly focused on droplet dynamics on a hydrophobic smoothed surfaces or on idealized granular surface. Since the shape of soil particles is irregular, further work is needed on realistic solid particle shapes. Moreover, the solid particles are generally fixed in the studies reported in this review. Therefore, incorporating the movement of soil particles is necessary, so that we can model erosional processes of soils including other phenomena in geotechnical engineering and earth surface processes.
In summary, if only considering the droplet dynamics on smoothed surfaces, all four methods can be used. If modelling irregular soil particles with complex surfaces, VOF is inferior to the other three methods. LBM is computationally intensive, and it is also not appropriate since the movement of solid particles should be incorporated in the future work. Although coupling LBM with discrete element method (DEM) can simulate the motion of solid particles, it is more straightforward to model the strong coupling of solids and fluids based on SPH. The application of MPM in this field is rare and haven't been fully verified. Therefore, in the follow-up study, we tend to use SPH method to study the droplet dynamics on hydrophobic soils.