Numerical Simulation of Droplet Impacting and Sliding on Hydrophobic Granular Surfaces

Droplet sliding naturally happens with practical significance in developing artificial self-cleaning surfaces or impermeable barriers. On water-repellent soil surfaces, such processes evolve at very small scales, typically at the particle level. To address this, this paper presents a two-dimensional Lattice Boltzmann (LB) study on the droplet sliding dynamics on a layer of regularly arranged particles with varying size and contact angle (CA) aimed at mimicking conditions comparable to those of real soils. ,e numerical droplet is initialized above the inclined granular surface with different lifting distances and deposited by gravity. ,e droplet hits the surface with different impacting velocities and subsequently slides down the slope. Four droplet-sliding behaviors were observed: a droplet sticks to the granular surface, a droplet moves by pinning and depinning of its interface (“stick-slip”), a droplet undergoes periodic elongation and shortening during sliding, and a droplet lifts off the granular surface and may be ruptured. For a droplet that displays the “stick-slip” behavior, the sliding velocity reaches a converged terminal velocity, which increases with a higher CA, a more inclined slope, and a smaller particle size. However, nonunique terminal velocities were identified to be affected by the impacting velocities, but their correlation is not continuous and may not be positive. Finally, we propose to quantify the rotational or translational movement by effective kinematic ratio (EKR), which is defined as the translational kinematic energy divided by the total kinematic energy. ,e unique relation between the EKR and the terminal velocity is suggested to be one practical indicator to intrinsically characterize the water repellency at the particle level.


Introduction
Droplet sliding on an inclined surface has practical significance in the development of artificial synthesized self-cleaning surfaces [1] such as antifouling clothes and kitchens [2]. In nature, wildfires [3] or soil microbiological activities [4] may also generate organic hydrophobic coatings on the soil particles, reducing the chemical affinity between the soil and the water, and thus making the soil waterproof. e waterproof property may be utilized for reducing the amount of water infiltration to improve the stability of slopes or landfills during rainfall. From the microscopic perspective, a droplet falling on an inclined soil surface will travel downslope, wherein the dynamics reveals the resistance of the soil to water infiltration. To better understand the fundamental mechanism behind, the paper utilizes an advanced numerical tool that can capture the evolution of droplet interfaces attached to one or several particles, and provides an insight into the droplet interaction and sliding dynamics with granular soil surfaces.
On an idealized, smooth, and homogeneous surface with a nonzero slope angle, the droplet slides down with a steady translational movement [5]. However, a realistic surface has microscale bulges or indentations, i.e., the rough substrate. Particularly, for granular materials, such as soils, the surface morphology is related to the size, shape, and the alignment of the particles [6]. In contrast to the idealized smooth surface, previous studies [7][8][9] indicate that a droplet is stuck on a rough surface with an inclined angle lower than a critical value, i.e., the sliding angle [8]. If the inclined angle exceeds the sliding angle, the droplet slides by regularly pinning and depinning of its interfaces on the substrate [10]. Similar stick-slip behaviors are also observed for a droplet sliding on a geometrically smooth surface with patterned varying wettability [11][12][13]. e stick-slip behavior triggers the oscillation of the droplet interfaces, resulting in a viscous dissipation. At the same time, gravitational potential energy contributes to the downslope movement. It is found that the stick-slip movement can reach a steady state with a constant sliding velocity [14], indicating that a balance is achieved between the viscous dissipation and the gravitational energy input. e constant sliding velocity, usually termed as the terminal velocity, is the focus of the studies on the droplet sliding dynamics.
Hyväluoma et al. [5] simulated the sliding dynamics for a droplet in the Wenzel state on a hydrophobic rough surface patterned by ridges. If a droplet maintains nearly a rounded shape during sliding, the linear relationship between the capillary number Ca � ρ]u drop /c lv , and the Bond number Bo � V 2/3 drop ρg/c lv proposed in [15] is verified, written as Ca ∝ Bo · sin θ sl − Bo c . Here, θ sl is the slope angle and g is the gravitational acceleration. Ca and Bo · sin θ sl indicate the dimensionless values of sliding velocity and the downslope gravitational driven force, respectively. e cutoff value Bo c reflects the minimum driven force in the tangential direction to the slope that can trigger the droplet sliding. Effects of wettability are also investigated. For a water droplet on a flat surface, the intersection angle between the interface and the base plane, i.e., the contact angle (CA) quantifies the wettability. A surface with a low CA smaller than 90°for the water droplet is considered as hydrophilic or wettable, while a high CA larger than 90°is considered as hydrophobic or water-repellent. In [5], by increasing the CA, Bo c decreases, and the ratio ΔCa/ΔBo increases, indicating that a higher CA can facilitate the droplet sliding. ey also found that it is easier for a droplet to move parallel than perpendicular to the patterned ridges. Wang et al. [16] numerically studied the effects of the solid fraction and the width of the micropillar for a droplet sliding in the Cassie-Baxter state on a patterned surface. ey found that increasing the spacing between the micropillars (i.e., decreasing the solid fraction) is beneficial for a droplet sliding down.
is is consistent with the experimental study in [8] that a larger slope angle is needed to initiate the droplet sliding for the surface with a higher solid fraction. However, further increasing the spacing between the micropillars may trigger the droplet configuration to transit from the Cassie-Baxter state to the Wenzel state [17,18]. In other words, the droplet penetrates the microscale indentations and the sliding is impeded. is indicates that there exists an optimal solid fraction for the patterned surface, on which the droplet achieves the maximum terminal velocity. Varagnolo et al. [13] used a high-speed camera to trace a droplet sliding on a geometrically smooth surface with alternating hydrophilic and hydrophobic strips. ey found that both the CA and the surface pattern of the domains with varying wettability affect the terminal velocity. e droplet undergoes a stick-slip motion and the terminal velocity is larger with a higher apparent CA. Zhao et al. [18] and Sbragaglia et al. [11] found that terminal velocity increases monotonically with the fraction of the hydrophobic regions. Concerning the detailed patterns, Varagnolo et al. [19] found that the terminal velocity depends mainly on the apparent CA and whether actual surface pattern is in the shape of strips, squares, or triangles has very limited effect. e above studies on the droplet sliding initiated the downslope movement by gradually increasing the slope angle of the tested surface, on which a droplet shall be deposited beforehand. In the real world, droplets may fall on an inclined surface with an impacting velocity. However, researches on the droplet impacting and rebounding dynamics only deal with the droplet distortion and the possible rupture [20,21] but ignore the subsequent sliding behavior. Up to now, the dependency of the terminal velocity on the impacting velocity remains unclear for a droplet, which can be clarified by conducting pore-scale numerical analysis. In the present work, the phase-field Lattice Boltzmann (LB) method is utilized as the numerical tool, due to its capability of capturing the evolution of numerous liquid-vapor interfaces interacting with irregular solid boundaries. On the contrary, traditional multiphase computational fluid dynamics (CFD) methods, e.g., the volume of fluid (VoF) method [22,23] and the level set method [24,25], can only trace the dynamics of a few predefined free surfaces in hydraulic problems, such as the open-channel flow. Another advantage is that surface tension and CA are accurately incorporated into the phase-field LB method by a thermodynamic approach. e feature is vital for accurately solving the dynamics of microfluids, where the capillary forces play an important role. By simulating the droplet impacting and sliding on inclined hydrophobic granular surfaces with varying CA and particle, various modes of droplet dynamics are identified with the sliding velocity monitored. In addition, we further propose an energy-based method to quantitatively characterize the rotational or translational type of droplet movement, by taking the advantages of the LB method that provides not only the global movement of the droplet but also the flow field inside the droplet.

Phase-Field Lattice Boltzmann Method
e numerical simulations of droplet sliding dynamics on a layer of particles need to deal with the movement of airwater interfaces on granular surfaces with controlled wettability. Such problems can be handled by adopting different multiphase LB methods, e.g., pseudopotential LB methods [26,27], free energy models [28][29][30], and phase-field LB method [31,32], noting that the LB method, originally developed from gas kinematics, is well-known for the efficient parallelization of the algorithm and flexibility in dealing with irregular boundaries. erefore, multiphase LB methods have been widely applied in simulating pore-scale flows, such as unsaturated permeability [33], wetting-drying in porous media [34], capillary rise [35], Pickering emulsion [36], and colloid hydrodynamics [37].
In the present paper, the simulation is performed with the phase-field LB method [32]. e system of binary fluids adopts the one-fluid formulation where a phase index field C is defined to describe the phase that occupies each computational cell. e phase index C is around one for the liquid phase and is around zero for the gas phase. e value gradually changes from zero to one in the interface. e evolution of the phase field is governed by the minimization of the free energy given by where Ψ is the total free energy, and the subscripts b and s represent the contributions from the bulk fluids and the solid boundaries, respectively. Ψ s is nonzero only for the fluid nodes adjacent to the solid nodes and C s is the composition at a solid surface and ϕ s (i.e., ϕ 1 , ϕ 2 , ϕ 3 , etc.) are constant coefficients. Ψ b consists of two parts, namely, the local bulk energy E 0 � βC 2 (1 − C) 2 and the interface-induced energy (κ/2)|∇C| 2 . e constants β and κ can be related to the surface tension σ and interface thickness ξ by Minimizing the bulk free-energy Ψ b for the fluid phase leads to the Cahn-Hilliard (CH) equation: where u is the velocity, µ is the chemical potential, and M is a numerical parameter controlling the rate of diffusion. e chemical potential is derived from the free energy as follows: where µ 0 is the local part of the chemical potential, and −κ∇ 2 C is the interface-induced chemical potential. Since E 0 takes the form of the double-well potential function, C will be close to either 0 or 1 in the pure single-phase region. e negative sign in the nonlocal term −κ∇ 2 C helps build a stable diffusive gas-liquid interface. e two competing terms mean that the phase index C is stable either completely separated as the pure single phase or well mixed as the interface. erefore, the evolution of the CH equation results in the phase separation while the total free energy in the system is decreasing.
For the fluid nodes adjacent to the solid surface, minimizing the Ψ b + Ψ s leads to a partially wetting boundary condition, where the CA can be defined at the fluid-solid interface by specifying the gradient of C [38]. e boundary condition is established as where θ eq is the equilibrium CA at the solid surface. Solving the CH equation (equation (4)) provides the spatial distributions of density ρ, the kinematic viscosity ], and the surface tension force F s . Here, density and viscosity are both necessary material properties for solving the incompressible Navier-Stokes (NS) equations. ey can be calculated from the index C simply through linear mapping, which are written as 1 where the subscripts l and g indicate the liquid phase and the gas phase, respectively. e surface tension force F s applied at the interfacial region is proportional to the gradient of the interface-induced chemical potential, given by e above equation is based on the chemical potential is called the "potential form" [31]. Together with the standard isotropic spatial discretion, the nonphysical velocities at the interfacial region in the multiphase LB method (also termed as "spurious currents") can be minimized [39].
With the information provided by the CH equation, the LB method computes the evolution of the mesoscale particle probability functions, from which the macroscale velocity and pressure are obtained. e standard discrete Boltzmann equation (DBE) with the body force term for macroscopic motion of fluid flows is where f α is the distribution function, e α is the discrete velocity for the direction α, c s is the lattice speed of sound, and λ is the relaxation time, which is related to the kinematic viscosity by ] � λc 2 s . F is the body force term F � F s + ρg, including the contributions from surface tension force and gravity. Γ α (u) � f eq α /ρ, and f eq α is the equilibrium distribution function defined by where t α is the weight for the direction α. e LB formulation has been proved to be equivalent to the artificial compression method in solving the incompressible NS equation by the Chapman-Enskog analysis [40]. However, the original formulation stores the base density with the hydrodynamic pressure as a small perturbation in the distribution Mathematical Problems in Engineering functions and applies only to the multiphase problem with a nonuniform density field. e problem is solved by defining a new distribution function g α � f α c 2 s + (p − ρc 2 s )Γ α (0) to split the density [41]. Similarly, the g eq α is transformed as g eq α � f eq α c 2 s + (p − ρc 2 s )Γ α (0). e DBE for the new distribution function is given as follows: from which the evolution of pressure and velocity in the macroscopic NS equations can be recovered as By integrating the DBE (equation (12)) with a trapezoidal approximation over a small time step δt [32], the Lattice Boltzmann equation for numerical implementation is where τ � λ/δt. With the temporal discretion, the truncation error is O(∆t 2 ). e velocity and the macroscopic pressure are obtained from the distribution functions by e same temporal discretion can also be applied to solve the CH equation [31,32]. However, the LB discretion can only recover the advection term ∇ · (uC) in equation (4), while the nonlinear diffusion term M∇ 2 μ cannot be recovered and related to the relaxation time λ. e advantages of the LB scheme to recover both the advection term and the linear diffusion as in solving the NS equations cannot be exploited. Instead, the M∇ 2 μ shall be included as an additional term, which still relies on the finite difference approximation with first-order accuracy. erefore, the CH equation can be more efficiently solved by a full finite difference (FD) scheme [42] and then coupled with the NS equations being solved by LB. We also adopt this hybrid LB-FD formulation. e advective term ∇ · (uC) is discretized by a second-order GVC (Group Velocity Control) scheme. In the Cartesian coordinate system with a uniform grid, the nonlinear diffusion term is approximated with the isotropic central difference scheme involving all the neighborhood nodes [42]. e third-order TVD (Total Variation Diminishing) Runge-Kutta scheme is used for the temporal discretion [43]. e computational details for solving the CH equation can be referred to computational fluid dynamics textbooks.

Numerical Test
e dynamic behavior of millimeter-sized droplet impacting and subsequently sliding on inclined granular surfaces is simulated with the LB method described in the previous section.
e numerical tool was validated by both the benchmark numerical tests and experimental tests, which can be referred to our previous work [44]. As it takes a significantly large number of time steps for a droplet to travel a long distance before reaching a steady state, a two-dimensional analysis is performed, so that the computation can be finished within a reasonable time.
e schematic diagram of the numerical setup is illustrated in Figure 1. e droplet diameter is set to 96 in LB units to mimic a 10-uL droplet, which is commonly used in the sessile drop method for quantifying the wettability of a granular surface [45]. e substrate is a layer of hydrophobic particles with varying CA, fixed on an inclined surface at different slope angles from 5°to 75°. e lattice nodes occupied by the particles are marked as solid nodes, implemented with the bounce-back boundary condition. e particle diameters are 24, 32, and 40 in LB units, corresponding to the physical values 0.67, 0.89, and 1.15 mm, respectively. e sizes are all within the range of coarsegrained sands, where the ratios D par /D drop are 1 : 4, 1 : 3, and 1 : 2.4. e relatively larger values of D par /D drop indicate that the details of the substrate geometry shall have a significant effect on the droplet dynamics, and hence a pore-scale numerical study is required. To facilitate pore fluid flows [46,47] in a two-dimensional analysis, artificial spacings are created to separate the particles. In all the simulations, the spacing ratio ϵ 0 , defined as the gap divided by the distance between two adjacent particles, is fixed as 0.25. e gap between two adjacent particles then approximately equals the equivalent circular diameter of the pores in a square packing. e CAs used in this work are 104.5°, 113.6°, 123.4°, and 134.4°, with the cosine values of −0.25, −0.4, −0.55, and −0.7, respectively. Only the hydrophobic particles with CA > 90°are studied, which ensures that the droplet will not infiltrate into the granular surface.
For each particle layer with different slope angles, the droplet hits the granular surface with multiple initial velocities. e impacting velocities are controlled by releasing the droplet from different distances above the granular surface so that the droplet accelerates under the gravitation force. e gravity is set as g � 4.27 × 10 −7 in the LB unit, corresponding to 9.81 m/s 2 in the physical unit. e falling distances H fall are chosen as 1/6, 2/ 3, 1, 2, 3, 4, 5 times the droplet diameter D drop , covering a wide range within the low-velocity impacts. Neither splitting nor splashing shall be induced. For the largest falling distance, the Weber number We, defined as the ratio of the inertial force over the surface tension ρu 2 D drop /c lv , is 9.84. e value is large to induce severe droplet deformation [24,25]. On the other hand, the smallest Weber number is only 0.328 for (H fall /D drop ) � (1/6), indicating a negligible inertial effect. e droplet almost remains circular during its deposition on the granular surface.
In each test, the liquid density is ρ l � 1, the liquid kinematic viscosity is ] l � 3.33 × 10 −3 , and the surface tension is σ � 4 × 10 −3 in LB units, corresponding to 10 3 kg/m 3 , 2.34 × 10 −6 m 2 s −1 and 72.8 mN/m, respectively, in physical units. To better illustrate the relationship between the LB units and the physical units, the input parameters are summarized in Table 1. Here, the kinematic viscosity is close to the value of 20% concentration saline water, which is 2.33 times that of pure water (1.006 × 10 −6 m 2 s −1 ). Further reducing the lattice viscosity ] l would trigger numerical instability, so a smaller physical viscosity can only be achieved by enlarging the computational domain, which would demand the exponential growth of the computational power. Despite the compromise, in terms of relative magnitude of inertial and capillary forces to viscous forces, the dimensionless Laplace number La � (σR/ρ l υ 2 l ) � 17280 ≫ 1 is sufficiently large so that the droplet is considered as having low viscosity [48]. Such a large value has been seldom numerically investigated in the literature. e simulation involves a preparation stage and a gravity deposition stage. First, the droplet is initiated in a circular shape above the particle layer with a certain distance to control the impacting velocity. e phase index C is predefined as 1 within the droplet and is predefined as 0 in the surrounding air. At the interface, a smooth transition is realized using a tanh function [32]. However, the distribution of the phase index C is not in the equilibrium state. e correct pressure field that fits the surface tension forces has not yet been built up. Hence, the simulation runs with no gravity until the system reaches equilibrium. en, gravity can be applied to the liquid phase so that the droplet moves downward and impacts the granular surface. At specific time intervals, the subsequent evolution process is recorded in the form of two-dimensional image files, showing the information of both the phase distribution and the velocity field. During the simulation, pressure boundary conditions are applied at the top and the bottom boundaries to maintain the ambient air pressure as constant. Periodic boundary conditions are imposed on the other two sides. With a limited computational domain, it is possible for a droplet to travel a sufficiently long distance to reach a steady state. As denoted by the dashed line in Figure 1, we notice that a smaller slope angle and a larger falling height require a larger computational domain. Hence, the width in the y dimension n y is determined by H fall cos α, plus the space to contain the particle layer and the gaps at the top and the bottom boundaries. e length n x in the sliding direction is defined as 480, which is 5 times the droplet diameter. e height n y is controlled by both the slope angle and the droplet falling height. e largest size is 480 × 720 for a 5°slope with 5D drop falling height, while smallest size is 480 × 250 for a 75°slope with (1/6)D drop falling height. For each configuration of CA and particle size, a complete sequence of 105 numerical tests with 15 slope angles and 7 impacting velocities require around 96 hours on an AMD R9 390 graphics processing unit (GPU) to finish the computation.

Droplet Sliding Dynamics on Granular Surfaces.
For clarity, a test is labelled by a combination of three segments. e first segment is the particle size and the intrinsic CA. e particle layer with the diameters 0.67 mm, 0.89 mm, and 1.15 mm is denoted by S, M, and L, respectively. e second segment, which begins with "sl," is the slope angle. e last segment, which begins with "h," is the falling height divided by the droplet diameter H fall /D drop . For example, "h1/6" means that the falling height is one-sixth of the droplet diameter and "h5" means 5 times the droplet diameter. erefore, "M113_sl60_h5" indicates that the droplet is falling from 5 times its diameter on a layer of mono-sized 0.89 mm spheres with an intrinsic CA of 113°, which is inclined at the slope angle 60°.
In the present paper, four modes for the droplet sliding on the inclined hydrophobic granular surfaces are identified: a droplet sticks to the granular surface (mode I), a droplet moves by pinning and depinning of its interface (mode II), a droplet undergoes periodic elongation and shortening during sliding (mode III), and a droplet lifts off the granular surface and may be ruptured (mode IV). Figure 2 illustrates the temporal evolution of the sliding velocity v drop of the droplet for L104_sl60_h5, M113_sl60_h5, S134_sl60_h5, and M134_sl65_h5, which are the representative cases for modes I to mode IV, respectively. In addition, for mode I and II, the evolution of the droplet shape in continuous  sequence showing the propagation of the droplet interfaces are plotted in Figure 3. For modes III and IV, since the droplet travels across the periodical boundaries for several times within a short period, we plot the snapshots in Figures 4 and 5 to show the droplet shape at selected moments. Mode I occurs if the droplet finally gets stuck on the slope after impacting and travelling a short distance. As shown in Figure 3, after impacting at 53.9 ms, the droplet spreads asymmetrically on the surface, developing a bulge at the downstream side and a thin tail at the upstream side. From 53.9 ms to 137.5 ms, inertia drives the droplet to move downslope. e interface undergoes pinning and depinning, from which the capillary wave (defined as a wave travelling along the phase boundary of a fluid) generates and propagates through the droplet. Severe distortion and oscillation are observed and hence, the kinematic energy is dissipated quickly and the sliding velocity slows down. Finally, the gravitational force, which drives the droplet down the slope, can be counterbalanced by pairs of surface tensions at the liquid-air, liquid-solid, and solid-air interfaces. From 137.5 ms, the droplet interfaces are fixed and the droplet is eventually stabilized.
Mode II occurs if the droplet sliding finally reaches a steady state featured by periodically pinning and depinning of its interfaces. Besides, the sliding velocity is constant, which is usually termed as the terminal velocity. e snapshots from 137.5 ms to 188.1 ms in Figure 3 demonstrate that regular pinning and depinning do not cause significant distortion of the droplet, and the interfacial areas for the liquid-air, liquid-solid, and solid-air interfaces remain dynamically constant at the steady state. In terms of its shape during the sliding, the style is called "stick-slip" movement [17]. At this moment, both the kinematic energy and the surface energy are dynamically stable, which indicates that the viscous dissipation is exactly compensated with the gravitational energy input.
Mode III in Figure 4 shows that the droplet undergoes periodic elongation and shortening, in terms of its shape during the sliding. e behavior is the transitional mode from mode II "stick-slip" to mode IV featured by the droplet completely lifting off from the granular surface. Compared to mode II, the gravitational energy input will not be counterbalanced by the viscous dissipation if the droplet is in a stick-slip movement. As a response to facilitate a higher rate of energy dissipation, the droplet undergoes elongation and shortening starting from 229.9 ms. For example, the droplet is shortened at 269.5 ms and elongated at 304.7 ms. Due to the difference in the surface energy, the sliding velocity is severely fluctuating.

Mathematical Problems in Engineering
Mode IV occurs in the extreme condition of a high CA with a high slope angle. e droplet lifts off the granular surface and jumps down over the slope. When the droplet is flying with no retardation from the granular surface, the droplet is accelerating at a high rate. Each time the droplet falls back on the granular surface (see ① and ③), Figure 5 shows that the higher impacting velocity will drive the droplet to be more and more distorted (see ② and ④). Finally, the droplet is ruptured to halves at 502.7 ms with no terminal velocity being achieved.

Effect of Impacting Velocity.
For the first time, the dependence of the terminal velocity on the impacting velocity is confirmed for a droplet that displays the stick-slip sliding dynamics. e terminal velocity is not unique, as indicated by the temporal evolution of the sliding velocity in Figures 6-9. Here, we select the representative cases S113_sl55, S123_sl30, M113_sl40, and S113_sl35, which all belong to mode II and converge to terminal velocities. It can be seen that most cases display multiple terminal velocities, although some converge to one single value. Figures 6 and 7 indicate that two terminal velocities are observed for the cases S113_sl55 and S123_sl30. Figure 8 shows that for M113_sl40, three terminal velocities are possible. In contrast, there can be one single value irrespective of the impacting velocity, for example, in S113_sl35. However, it should be noted that the terminal velocity does not have a continuous positive-correlated relationship with the impacting velocity. With as many as seven different impacting velocities, only a few discrete values of the terminal velocities can be obtained. e phenomenon of a droplet sliding in multiple discrete terminal velocities can be understood through an energy approach. e internal energy related to the droplet sliding dynamics consists of the kinematic energy and the surface energy. In the steady state, the constant terminal velocity indicates that the rate of change in the kinematic energy ΔE kine is zero. e rate of change in the surface energy depends on the magnitude of the surface tensions c lv , c sl , and c sv , multiplying the rate of change in the surface areas ΔΩ lv , ΔΩ sl , and ΔΩ sv , for the liquid-vapor, solid-liquid, and solid-vapor interfaces, respectively. e ΔE surf is written as On average, the segments of the main liquid-vapor droplet interface and the particle surface attached to liquid ) ) S113_sl55_h1/6 (1) S113_sl55_h2/3 (2) S113_sl55_h1 (3) S113_sl55_h2 (4) S113_sl55_h3 (5 S113_sl55_h4 (6 S113_sl55_h5 (7) (1) (2) (3) Time (ms) Figure 6: Temporal evolution of the sliding velocity for the series of case S113_sl55 with two terminal velocities. S123_sl30_h1/6 (1) S123_sl30_h2/3 (2) S123_sl30_h1 (3) S123_sl30_h2 (4) S123_sl30_h3 (5) S123_sl30_h4 (6) S123_sl30_h5 (7) (1) or vapor phases remain dynamically constant, indicating that ΔE surf is also zero. Next, the energy input and loss are analyzed. e gravitational energy W grav is the only input energy source, written as where u dr op and A dr op are the sliding velocity and the area of the droplet, respectively. e viscous dissipation rate is Equations (19) and (20) indicate that the higher sliding velocity at a higher energy level results in a higher gravitational energy input rate, and also higher viscous dissipation rate, which increases with the shear rate of the liquid inside the droplet. e opposite is true for a lower energy level. Despite various energy levels, since ΔE kine and ΔE surf are both zero, the gravitational energy input and the viscous dissipation shall always cancel each other out. Otherwise, the sliding state is unstable, triggering the droplet acceleration or deceleration through its interaction with the granular surface. When the sliding is in such an unstable state, even a small discrepancy of the initial condition may have great influence on the converged energy level. Figure 6 shows that the sliding velocities are very close to each other during the period from 148.5 ms to 336.6 ms, for the cases S113_sl55_h1/6 and S113_sl55_h2/3. However, the small difference causes their divergence at 336.6 ms and S113_sl55_h2/3 turns into a higher energy level by a sudden acceleration. It is worth noting that the droplet may be able to run in an unstable state without significant changes in the velocity, for a relatively longer duration. For instance, Figure 7 shows that the velocity is almost constant from 337.7 ms to 608.3 ms in the case S123_sl30. Finally, the unstable state breaks up and the sliding velocity transits into a higher energy level.
Another important finding is that the terminal velocity may not be positively correlated with the impacting velocity at relatively low slope angles. To better clarify this finding, Figure 10 plots the terminal velocities against various impacting velocities for the two illustrative series of M113 and M134. For example, in the series of M113_sl40 with three energy levels, "h2/3" and "h1" correspond to the low energy level; "h1/6" and "h3" correspond to the middle energy level; and "h2", "h4," and "h5" correspond to the high energy level. We see that a higher impacting velocity may not necessarily lead to a higher energy level for a relatively lower slope angle. is is because the impact is more developed in the normal direction to the slope as the spreading and recoiling processes on the complex granular surface. e relatively small downslope component of the impacting velocity does not efficiently drive the droplet to slide down at a higher velocity. Moreover, the recoiling process results in a larger viscous dissipation rate, which retards the downslope movement of the droplet. However, in the case of a higher slope angle, for example, M113_sl60, "h3", "h4," and "h5" slide at the higher energy level, while all the other cases with lower impacting velocities correspond to the low energy level. e relatively larger downslope component of the impacting velocity can overcome the energy barrier, so that the droplet tends to slide at a higher energy level. A similar phenomenon can also be observed in the series M134 with a higher CA in the figure. However, the threshold slope angle for the terminal velocity to be positively correlated with the impacting velocity is only 15°, compared with the value of 60°i n M113.  (7) (1) (3) S113_sl35_h1/6 (1) S113_sl35_h2/3 (2) S113_sl35_h1 (3) S113_sl35_h2 (4) S113_sl35_h3 (5) S113_sl35_h4 (6) S113_sl35_h5 (7) (1) (3) 200 400 600 800 1000 0 Time (ms) Figure 9: Temporal evolution of the sliding velocity for the series of case S113_sl35 with a single terminal velocity. e terminal sliding velocities are summarized in Figure 11 for all the combinations of CA and particle size. At each slope angle, since multiple terminal velocities may be achieved for the simulations conducted with seven different impacting velocities, the values are averaged together with the error bars indicating the upper and lower bounds. e figure shows that the terminal sliding velocity increases with the slope angle since the energy input is higher. For example, in S134, the average terminal velocity increases from 0.0178 m/s on a 5°slope to as large as 0.534 m/ s on a 45°slope. e CA is found to have a significant positive effect on the sliding velocities. A larger CA causes less distortion and oscillation every time the depinning of the droplet interface happens at the upstream, leading to less viscous damping and a higher velocity. Take the 0.67 mm particle layer inclined at 40°as an example; the droplet sticks to the granular surface for the CA 104°, irrespective of the impacting velocities. With larger CAs 113°, 123°, 134°, the average terminal velocities are 0.0619 m/s, 0.277 m/s, and 0.492 m/s, respectively. e CA also affects the static behavior of critical slope angle to trigger the downslope sliding. As the depinning of the droplet interface is more difficult, an inclined surface with a smaller CA can retain the droplet statically at a larger slope angle. For example, the minimum slope angles to induce the droplet sliding for the particle size 0.67 mm are 5°, 20°, 35°, and 60°for the CAs 134°, 123°, 113°, and 104°, respectively. e effect of particle size in terminal velocity is small but not negligible, within the relatively narrow range of coarse sandy sizes. Figure 11 shows that at the same CA, the terminal velocity is generally smaller for a larger particle size, which may be interpreted as the larger geometric roughness and hence result in a higher energy dissipation rate. is is consistent with Wang et al. [20] where increasing the width of micropillars slows down the droplet sliding. Next, the larger geometric roughness for the larger particle creates a larger energy barrier between the higher and the lower energy levels. In other words, it is possible for the droplet to run in significantly different terminal velocities on the particle layer formed by larger particles, as indicated by the larger error bars. Figure 11 also shows the deviation between the upper and lower bounds of the multiple terminal velocities for each configuration of CA and particle size. e deviation is generally larger for the larger particle, except for L104 with only a few samples. Take CA � 113°for instance, the maximum deviations are 0.0816 m/s at the slope angle of 50°, 0.0934 m/s at the slope angle of 40°, and 0.121 m/s at the slope angle of 65°, for the particle sizes 0.67 mm, 0.89 mm, and 1.15 mm, respectively.

Quantification of the Rotational or Translational Behavior.
In this subsection, we propose a new method to quantify the motion style to describe whether a droplet is undergoing rotational or translational movement. For a droplet sliding at a lower velocity, the droplet appears to rotate over the particles in a more rounded shape. With a higher sliding velocity, the droplet is more elongated and slides in a translational motion. Previous studies describing the droplet rotation or translation need to assume a droplet radius so that the velocity can be decomposed into the rotational component and the translational component [18,49]. However, this approach only applies to a droplet that maintains its rounded shape during sliding.

Mathematical Problems in Engineering
From the perspective of energy, we propose to use the ratio of the translational kinematic energy in the total kinematic energy to quantitatively describe whether a droplet is rotating or translating. We define it as the effective kinematic ratio EKR, which is written as u drop and e drop are the sliding velocity and the kinematic energy unit area of the droplet as a whole. ey are averaged over the region that is occupied by the liquid phase (C > 0.5), based on the Euler description of the fields of phase index C and the velocity u. Such a kinematic ratio from the energy perspective can be applied to a droplet with an arbitrary shape. e value ranges between 0 and 1, from rotation to translation. Figure 12 plots the effective kinematic ratio EKR against the sliding velocity for CA � 104°, 113°, 123°, and 134°. EKR increases with the sliding velocity in a high rate at low velocities but finally converges to a constant value. Hence, a hyperbolic equation EKR � (u drop /a + bu drop ) is adopted to empirically fit the relationship. e CA � 104°is not fitted since it does not have enough data points. e particle size has a negligible effect on the EKR. is may due to the fact that the selected particle sizes are all within a narrow range of coarse sandy sizes. Hence, for each CA, the results from all the particle sizes are analyzed together.
e table in Figure 12 shows that for the CAs of 113°, 123°, and 134°, the values of a are 0.05978 m/s, 0.04316 m/s, and 0.02227 m/s, respectively, while the values of b are 1.231, 1.218, and 1.187, respectively. Considering the geometric interpretation, the inverse of parameter a means the rate of effective kinematic ratio increasing with the terminal velocity at the point of origin, while the inverse of parameter b means the maximum effective kinematic ratio that the droplet is converging to. e results indicate that the droplet sliding is getting closer to the rotational motion with a lower CA.
is is because the depinning of the droplet interface is more difficult for a lower CA, and therefore additional kinematic energy is required to drive the droplet to rotate over the particles. Compared to the multiple terminal velocities that depend on the impacting velocity, we demonstrate that the unique relationship between the EKR and the terminal velocity intrinsically characterizes the particle-level water repellency.

Conclusions
e present work investigated the interaction between an inclined hydrophobic granular surface and a millimeter-sized droplet. e droplet dynamics, including the impacting and the subsequent downslope sliding, was simulated with the LB S104 M104 L104 S113 M113 L113 S123 M123 L123 S134 M134 L134 method. Granular surfaces with different geometric configurations and water repellency were constructed. To elucidate the mechanisms controlling the interaction of droplets and hydrophobic granular slopes, a series of two-dimensional numerical droplet sliding tests were carried out on a granular surface formed by a single layer of spherical particles with different sizes and intrinsic CAs. Particularly, the droplet was impacting the granular surface with different velocities by controlling the falling height. e major findings can be drawn as follows: (1) Under different conditions, the droplet dynamics can be described within four modes in the downslope sliding: (I) a droplet sticks to the granular surface, (II) a droplet moves by pinning and depinning of its interface, (III) a droplet undergoes periodic elongation and shortening during sliding, and (IV) a droplet lifts off the granular surface and may be ruptured. A higher CA and a larger slope angle favor the downslope sliding from mode I to mode IV. e higher slope angle means a higher gravitational energy input, and a larger CA causes less energy dissipation, since the easier depinning process triggers less distortion and oscillation. (2) Mode II, a droplet moving steadily by pinning and depinning of its interface (stick-slip) can reach a converged terminal velocity through droplet-particle interaction, which means the viscous dissipation is exactly balanced with gravitational energy input. (3) e dependence of the terminal velocity on the impacting velocity is confirmed, and the relationship is not continuous. We first demonstrate that the droplet has multiple discrete choices of its terminal velocity under different impacting velocities. Moreover, the terminal velocity is positively correlated with the impacting velocity only at relatively high slope angles, since the large downslope component of the impacting velocity can overcome the energy barrier, and the droplet tends to slide at a higher energy level. At low slope angles, the convergence to a higher or lower terminal velocity is unpredictable, since the impact develops in the normal direction to the slope, featured by the spreading and recoiling. (4) e terminal velocity increases with both the slope angle and the CA, which, respectively, indicates a higher gravitational energy input and a lower viscous energy dissipation. (5) e larger geometric roughness induced by a larger particle size results in a higher energy dissipation rate during the droplet-particle interaction, and therefore a smaller terminal velocity. e larger roughness also creates a larger energy barrier, so the droplet may converge to significantly different terminal velocities on granular surfaces formed with larger particles. (6) e effective kinematic ratio, EKR, defined as the translational kinematic energy divided by the total kinematic energy, is proposed as the quantitative indicator to describe a rotational or translational type of movement. Within the "stick-slip" sliding mode, the EKR increases with the terminal velocity at a high rate at small velocities but finally converges to a maximum value. An empirical equation (v drop /a + bv drop ) can be adopted to empirically fit the relationship. e fitted results capture the phenomenon that for a larger CA, the EKR increases faster with the terminal velocity and the maximum EKR that can be achieved is also larger. Besides, compared to multiple terminal velocities, which may be affected by the impacting velocity, the unique relationship between the EKR and the terminal velocity is a better indicator to intrinsically characterize the particle-level water repellency.
e current study provides a basis to investigate the factors that affect the efficient removal of a droplet under its gravity on an inclined granular surface, which is fundamental to applications of synthetic water-repellent soils in the geotechnical engineering discipline. e study can also extend to other applications involving self-cleaning rough surfaces, for e.g., water-repellent clothes and antifouling kitchenware.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.