Surface morphology effects on clathrate hydrate wettability.

HYPOTHESIS
Clathrate hydrates preferentially form at interfaces; hence, wetting properties play an important role in their formation, growth, and agglomeration. Experimental evidence suggests that the hydrate preparation process can strongly affect contact angle measurements, leading to the different results reported in the literature. These differences hamper technological progress. We hypothesize that changes in hydrate surface morphologies are responsible for the wide variation of contact angles reported in the literature.


EXPERIMENTS
Experimental testing of our hypothesis is problematic due to the preparation history of hydrates on their surface properties, and the difficulties in advanced surface characterization. Thus, we employ molecular dynamics simulations, which allow us to systematically change the interfacial features and the system composition. Implementing advanced algorithms, we quantify fundamental thermodynamic properties to validate our observations.


FINDINGS
We achieve excellent agreement with experimental observations for both atomically smooth and rough hydrate surfaces. Our results suggest that contact line pinning forces, enhanced by surface heterogeneity, are accountable for altering water contact angles, thus explaining the differences among reported experimental data. Our analysis and molecular level insights help interpret adhesion force measurements and yield a better understanding of the agglomeration between hydrate particles, providing a microscopic tool for advancing flow assurance applications.


Introduction
Gas hydrates naturally form under certain conditions of pressure and temperature because of a balance between water-water hydrogen bonds and host-guest (i.e., water-methane) dispersive interactions [1]. Over the last decade, gas hydrates have attracted significant focus because of their possible applications in various sustainable and environmental technologies [2][3][4][5][6][7][8][9]. Natural gas hydrates could also provide an alternative energy source; notwithstanding their environmental effects should be understood, prevented and mitigated when necessary [8,10]. The uncontrolled formation of hydrate plugs in oil/gas pipelines [1] and offshore energy operations [11] also causes disruptions in oil and gas production, leading to accidents, as well as large negative environmental consequences [12][13][14].
As the formation of hydrates occurs at interfaces, it has been proposed that hydrate formation and agglomeration could be inhibited by altering the water wettability of hydrates [15,16] using low dosage hydrate inhibitors, thereby preventing the interaction/coalescence among hydrate particles and water droplets [12,[17][18][19][20]. Understanding hydrate wettability is therefore crucial in optimizing flow assurance activities, where a remarkable change is taking place from ''hydrate avoidance" to ''hydrate management" [21,22]. At the macroscopic scale, wetting properties can be assessed by measuring the Young contact angle of a liquid droplet on a rigid and smooth surface [23]. However, this is challenging for hydrates, which form at high-pressure & low-temperature conditions following processes that are stochastic in nature [24]. Notwithstanding many practical obstacles, recent research advances have been achieved in the direct measurement of contact angles for hydrate systems [16,25]. For example, Brown et al. [16] reported an average water contact angle on a cyclopentane (C5) hydrate in liquid C5 of 94.2°± 8.5°at atmospheric pressure and temperatures of 273 -280 K. However, the reported contact angles vary significantly in the literature, perhaps because the preparation of the hydrate surfaces impacts the measured contact angles. For example, in contrast to Brown et al. [16], who used an annealing time of 30 min, Thomas et al. [25] implemented temperature cycling protocols to create C5 hydrates over a period of $24 h, inspired by the experimental procedures employed by Zylyftari et al. [26] Using the latter substrates, Thomas et al. [25] studied the spreading dynamics of a water droplet and showed a fully wetted C5 hydrate 7 s after droplet deposition. Supporting the possible effect of preparation on wetting properties, Li et al. [27] showed significant differences in the hydrate surface morphologies when changing gas compositions and degree of subcooling in the process of creating ethane and methane-ethane (C1-C2) hydrates. It is possible that both hydrate porosity and surface roughness affect the spreading dynamics of a water droplet on a hydrate surface, explaining the different results reported in the literature [16,25,28]. Reconciling the differences in wetting properties is essential, as these fundamental interfacial properties are related to the cohesive/adhesive forces responsible for hydrate agglomeration and pipelines plugging. It is expected that hydrates formed in practical field applications would differ in morphology and surface features from hydrates formed in the laboratory under controlled physical conditions.
Because of the practical hurdles faced by experimental approaches, we implement here atomistic molecular dynamics (MD) simulations to probe the wetting properties of atomically smooth and rough C5 at atmospheric pressure. As realistic flow assurance activities occur at higher pressure and lower temperature, here we also examine the wettability of C1-C2 hydrates in liquid hydrocarbons under high-pressure & low-temperature conditions (3.45 MPa and 274 K). We quantify contact angles of a water droplet on hydrate surfaces with different morphologies, surface tensions of hydrate-liquid and liquid-liquid interfaces, contact line pinning forces as induced by surface heterogeneity, and work of adhesion (these quantities are defined later) between water droplets and hydrate particles. The simulation results, in general agreement with several experimental observations, help us reconcile at the molecular-level the differences among contact angles reported in the literature, and further highlight the applicability of simulated interfacial properties to examine adhesion force measurements and predict properties of experimental relevance.
In what follows, we first describe the simulation methodologies and force fields implemented, then discuss our simulated results, compare them against a variety of experimental datasets, and summarize our fundamental observations with main practical implications.

Methods
Model Setup. Our simulation setup (see Fig. 1A) mimics the sessile contact angle measurement for a water droplet on a hydrate surface in the presence of a liquid hydrocarbon. In this study, two systems are considered: (1) a water droplet on a C5 hydrate immersed in bulk liquid C5 at atmospheric pressure and 275 K, and (2) a water droplet on a C1-C2 hydrate surface immersed in bulk liquid n-dodecane (C12) at 274 K and 3.45 MPa. The sII hydrate structure was employed as the solid substrate for both C5 and C1-C2 hydrates. To construct the hydrate configurations, we adopted the sII unit cell from Takeuchi et al. [29]. Tiling the unit cell in X and Y directions, the X , Y and Z dimensions of the hydrate substrate were 10.386, 5.193, and 1.731 nm, respectively. We carried out simulations with a droplet of 940 water molecules starting from various initial configurations (such as rectangular or cylindrical) placed on the hydrate substrates and surrounded by hydrocarbon (C5/C12) molecules. The resultant simulation box length was $ 10.8 nm in the Z direction, with the hydrate substrate aligned parallel to the X-Y plane. The hydrate substrate model becomes infinite along the X and Y directions when applying periodic boundary conditions in all directions.
Upon equilibration, we used a cylindrical water droplet periodic along the Y direction -the cylinder axis, as this configuration guarantees that the line tension of the three-phase boundary does not influence contact angle estimates from simulations [30,31].
To quantify the hydrate-liquid interfacial free energy, we simulated systems in which a thin film of liquid molecules, i.e., C5, C12, or water, is placed on the hydrate substrate (see Fig. 2A). The number of water, C5, and C12 molecules on the hydrates were 7000, 2600, and 1200, respectively, yielding a water film of $ 4 nm and a liquid hydrocarbon film of $ 9 nm in thickness. In these simulations, we used the same sII hydrate substrates employed in the contact angle measurements, with their surface parallel to the X-Y plane. Because the length of the simulation box in the Z direction was 17.731 nm, and the thickness of solid substrate was 1.731 nm, there is an empty gap between the periodic image of the hydrate substrate and the thin film [32,33]. To evaluate liquid-liquid interfacial tensions, we conducted simulations for systems composed of the water-hydrocarbon (C5, C12 or C5toluene) interface at selected temperature (T) and pressure (P) conditions. The water-hydrocarbon interface model was generated by combining aqueous and hydrocarbon phases into 4.0 Â 4.0 Â $13.0 (nm 3 ) simulation boxes, as shown in Figure S1 in Supporting Information (SI). Additional simulations for systems with waterhydrocarbon interfacial areas of 5.0 Â 5.0 nm 2 were conducted to verify that using smaller surface areas (4.0 Â 4.0 nm 2 ) was sufficient to extract reliable liquid-liquid surface tension results. The results were comparable. We have employed similar approaches to quantify the surface tensions of liquid-liquid interfaces for a Fig. 1. A) Representative simulation snapshot for the final configuration for a water droplet on a C5 hydrate surface immersed in C5 at 275 K and 0.1 MPa. Red and white spheres represent water oxygen and hydrogen atoms, respectively. Cyan, grey, and blue wireframes represent C5 in the liquid, C5 in the hydrate, and water in the hydrate, respectively. B) 2D density profiles (top) and circular best fit (bottom) of the water droplet on the C5 hydrate surface in the X-Z plane. The color bar shows water density in the units of 1/nm 3 . It should be noted that, in our representation, the purple regions are those where no water molecules of the droplet are present, i.e., these regions could be occupied by molecules from the oil phase and/or the hydrate surface. C) Spreading dynamics of the droplet on the C5 hydrate surface in C5 at ambient pressure and 274.7 K as observed experimentally. Adapted with permission from Ref. [25]. Copyright 2021 Elsevier B.V. Fig. 2. A) Representative simulation snapshot for a system composed of a thin C5 film on a C5 hydrate surface at 275 K and 0.1 MPa. Cyan and white spheres denote C5 carbon and hydrogen atoms, respectively. Grey and blue wireframes represent C5 and water molecules in the hydrate phase, respectively. B) Profiles of three components of the stress tensor r(z): r xx and r yy are parallel to the surface while r zz perpendicular to it. C) Surface tensions of hydrate-C5, hydrate-water and C5-water interfaces, and water contact angle on the C5 hydrate surface in C5 calculated using Young's equation using data from the present simulations (left) and from experiments (right) [60].
Force Fields. Water was characterized by the TIP4P/Ice model [37], which has been used successfully in studying hydrate systems [18,38,39]. Methane, ethane, and n-dodecane were represented by employing the united-atom version of the TraPPE-UA force field [40], which correctly describes the critical properties and the vapor-liquid coexistence of linear alkanes far from the critical point. Cyclopentane (C5) and toluene (C7) were modelled implementing the General Amber Force Field (GAFF), often employed to study cyclic, organic, and pharmaceutical compounds containing H, C, N, O, S, P, and halogens [41,42]. Previous computational studies suggest that the combination of GAFF force field for modelling cyclic and organic molecules, the TIP4P/Ice model for describing water, and TraPPE-UA force field for simulating linear alkanes yields excellent agreement against experimental studies [18,19]. All non-bonded interactions were described via both electrostatic and dispersion forces using the Coulombic and 12-6 Lennard-Jones (LJ) potentials, respectively, with the cutoff distance of 14 Å. We utilized the particle-particle particle-mesh (PPPM) approach for the treatment of long-range corrections [43] and the Lorentz À Berthelot combining rules to quantify unlike LJ interactions [44].
Algorithms. Equilibrium MD simulations were performed employing the GROMACS package [45], version 2016.3. To quantify contact angles on hydrates, we first performed simulations in the NVT canonical ensemble (constant volume, temperature and number of particles) for 1 ns to relax the initial configurations, while the hydrate layer was kept rigid, and then the simulations were carried out within the NPT ensemble (T = 275 K and P = 0.1 MPa for water droplet on C5 hydrate in C5 as well as C5-C7 and T = 274 K and P = 3.45 MPa for water droplet on C1-C2 hydrate in C12) implementing the Nose-Hoover thermostat and the Berendsen/Parrinello-Rahman barostat [45]. We applied the pressure coupling algorithm only along the Z direction, which allows X and Y dimensions of the simulation box to be kept constant.
For solid-liquid interfacial free energy calculations, all simulations were performed only within the NVT canonical ensemble while the simulations for the estimation of liquid-liquid interfacial tensions were conducted in the NVT and subsequently NPT ensembles.
Using the leapfrog algorithm the equations of motion were solved with the time step of 1.0 fs [45]. We implemented a harmonic restraint force constant of 2000 kJ/mol/nm on water, C5 as well as C1-C2 molecules in the hydrate phase to tether them to their initial positions [45] while other molecules in the system move freely. We conducted each NPT simulation for ! 150 ns until both fluids appeared to be stable, and the droplet shape remained unchanged within a simulation time interval of 20 ns.

Results and discussion
Wettability of Cyclopentane Hydrates immersed in Cyclopentane. We extracted contact angles from 2D density profiles obtained for the simulated water droplets (see Fig. 1B, top).
The isodensity contours at water density q o , obtained as halfway between the water density in the hydrocarbon phase and the water bulk density [31] were used to determine the contact angle for all systems considered. Once the droplet contours were identified following the procedure described in previous studies [31,46], we fit them with a circular function (see Fig. 1B, bottom). The droplet base was identified at the second hydration layer away from the hydrate surface, and the slopes of the tangent lines on both sides of the droplet were determined to calculate the contact angle.
Our simulation results yield contact angles of $ 21°on the C5 hydrate at 275 K and 0.1 MPa, which is in good agreement with macroscopic experimental data reported by Thomas et al. [25] (18.3 -23.1°) under similar T and P conditions (see Fig. 1C). These results suggest the C5 hydrate is remarkably hydrophilic, as expected, as it is largely made of water and it is structurally similar to ice. We also conducted additional simulations for a system in which the length of hydrate slab in the X direction was enlarged to L x = 15.579 nm to verify that periodic boundary conditions did not significantly affect the results presented. The water contact angle found on the enlarged hydrate surface was similar to the one obtained on the smaller one, suggesting that the systems are large enough to minimize the effect of periodic boundary conditions on the water contact angle. It should be noted that the contact angle of a simulated spherical nanodroplet [47] or of a cylindrical droplet [48] depends on the droplet size. Based on our experience, the system size used here is sufficiently large to minimize system-size effects on the estimated contact angles [30,31]. The simulation results also depend on the determination of the position and geometry of the water-liquid and water-solid interfaces, which is somewhat arbitrary. Therefore, these calculations are subject to some uncertainty, especially for highly hydrophilic surfaces [49].
As an alternative, which in some cases could yield more reliable results, one could calculate the contact angle directly from Young's equation: which relates the contact angle to the surface free energies of solid-liquid (c sf2 and c sf1 ) and liquid-liquid (c f1,2 ) interfaces [23,50]. Simulated Interfacial Tensions. Molecular modelling has been shown promising in investigations of interfacial forces at solid-liquid interfaces and in fact, it has been widely employed to estimate surface free energies [51][52][53][54]. In recent years, multiple approaches have been proposed to quantify solid-liquid interfacial free energy. These methods can be categorized into either thermodynamic or mechanical routes. The thermodynamic approaches, which provide accurate results, require large sets of simulations at a huge computational cost [51]. Alternatively, the surface free energy can be determined from the mechanical approach, where it is extracted from the interfacial stress anisotropy [51]. Kirkwood and Buff [55] were the first to describe the stress tensor components against the intermolecular forces, a method which is broadly used to calculate liquid-liquid interfacial tensions [56,57]. This approach, however, is challenging when applied to solid-liquid systems [53,54], because of the anisotropy of the solid surface, which yields an alternative interpretation of the stress tensor and the surface free energy in accordance with the Shuttleworth relationship [58], applicable only when the transverse stresses have no effects on the interface. Because the hydrate surface is kept rigid in our work, such limitations can be overcome. Hence, because of its effectiveness in comparison to the thermodynamic one, in this study we employ the mechanical approach following the Irving-Kirkwood formalism [59]. This method presents a local definition of the stress tensor when velocities and interaction forces of all atoms are known.
For a system built with isotropic and planar symmetrical surfaces, such as a flat hydrate substrate exposed to water/liquid hydrocarbon, the stress tensor r(z) can be described with its three components: r zz is perpendicular to the surface, and r xx and r yy are parallel to it. To quantify the local r(z), we employed the stresses at each frame, as it has been suggested that the stress profiles show little difference beyond this cutoff distance [61]. To compute normal and tangential stresses using GROMACS-LS, we divide the simulation box into multiple rectangular grids with cell size L grid = 0.1 nm in three directions and calculate the averaged r (z) over each grid cell. In Fig. 2B, we report the three components of stress tensor profiles r xx (z) (top), r yy (z) (middle), and r zz (z) (bottom), which arise from the local interaction forces exerting on the liquid C5 molecules in the normal direction of the C5 hydrate surface at 275 K and 0.1 MPa.
Considering the stress tensor profile along the Z direction perpendicular to the hydrate surface, we describe P N z ð Þ ¼ Àr zz ðzÞ, corresponding to the normal pressure along the Z direction, while P T z ð Þ ¼ À 1 2 r xx z ð Þ þ r yy ðzÞ À Á , which is the tangential pressure. Using these quantities, we obtain the surface tension c sl through the computation of the stress tensor profile [59]: In Fig. 2C, left panel, the results of surface tensions of hydrateliquid and liquid-liquid interfaces are shown. We estimated the surface free energies of hydrate-C5 and hydrate-water interfaces as 42.8 ± 2.6 and 0.24 ± 0.16 mN/m, respectively. The simulated hydrate-C5 interfacial free energy is remarkably consistent with the estimates of Aman et al. [60] (47 ± 5 mN/m), who employed a revised cohesive force model combined with experimental data obtained from micromechanical force measurements (MMF) between C5 hydrate particles. Our simulation results also compare well with an approximation based on wetting angle data proposed by Kwok and Neumann [64,65], who reported a value of 45 mN/m. It is worth pointing out that these hydrate-C5 interfacial free energy values are comparable to the C5-water interfacial tension data (46.0 ± 2.4 mN/m obtained from the simulation results presented here and $51 mN/m from literature experiments [60]) under similar T and P conditions. This agreement suggests that the hydrate surface in contact with a liquid hydrocarbon may exhibit similar energetic properties as a water surface. Aman et al. [60] also reported a hydrate-water interfacial free energy of 0.32 ± 0.05 mN/m (see Fig. 2C, right), which is close to our findings (0.24 ± 0.16 mN/m). For completeness, we point out that this estimate for the surface tension of hydrate-water interface is two orders of magnitude less than some literature data, in the range of 14 -45 mN/m, reported for methane, ethane, and propane hydrates [66][67][68]. However, these estimates were indirectly determined from porous media measurements, where huge discrepancies in wetting angle and pore size may lead to uncertainties [69].
(1)] with the values of computational and experimental surface free energies shown in Fig. 2 yields the equilibrium contact angles of 22.3°and 23.8°, respectively, for the droplet on the C5 hydrate surface immersed in C5. These values are in excellent agreement with the data obtained directly from our simulations (see Fig. 1B, top), as well as from recent experimental observations [25] (see Fig. 1C).
Wettability of Methane-Ethane Hydrates Immersed in n-Dodecane. Most of the wetting studies on clathrate hydrates were conducted on hydrates stable at atmospheric pressure, i.e., using C5 hydrates [16,25], even though realistic flow assurance applications typically encounter high P and low T (i.e., > 2 MPa and < 277 K) in liquid-hydrocarbon-dominated phases [70]. To examine somewhat realistic conditions, we investigate the wetting properties of C1-C2 hydrates immersed in bulk C12 at 274 K and 3.45 MPa (see Fig. 3A). After analysing 2D density profiles obtained for the simulated water droplets (Fig. 3B, top), we determined the contact angles (Fig. 3B, bottom). Our simulations yield a contact angle of 34°on C1-C2 hydrates immersed in C12, which agrees well with the simulated data reported by Naullage et al. [71] (34 ± 2°), who considered sI methane hydrates at slightly higher T (277 K) and P (10 MPa) conditions. We also simulated the water contact angle on C1-C2 hydrates immersed in C12 under atmospheric pressure. The result ($30°) is comparable to data obtained at higher pressures, suggesting that moderate changes in pressures up to $ 20 MPa only marginally impact the wetting properties of hydrate surfaces, a result which is consistent with data from other studies [72,73]. Our results suggest that, similarly to the C5 hydrate, the C1-C2 hydrate is also highly hydrophilic. For completeness, it should be mentioned that Molinero et al. [74] showed that wetting a clathrate hydrate surface with a thin water film Fig. 3. A) Representative simulation snapshot for the final configuration for a water droplet on the C1-C2 hydrate surface immersed in C12 at 274 K and 3.45 MPa. Red and white spheres denote water oxygen and hydrogen atoms, respectively. Grey spheres and blue wireframes represent C1-C2 and water in the hydrate, respectively. Grey wireframes denote C12 molecules in the bulk phase. B) 2D density profiles (top) and circular best fit (bottom) of the droplet on the C1-C2 hydrate surface along the X-Z plane. The color bar shows water density in the units of 1/nm 3 , similar to those shown in Fig. 1B. C) Simulated surface tensions of hydrate-C12, hydrate-water, and C12-water interfaces, as well as the corresponding water contact angle on the C1-C2 hydrate in C12 as calculated using Young's equation from the surface tension data. could lead to strong adsorption of alkanes, suggesting that the resultant interface is oleophilic.
(1)] with the simulated water contact angle of 34 ± 2°, the hydrate-water interfacial free energy of 33 ± 4 mN/m [67,75], and the experimental C12-water interfacial tension of 53 mN/m [76], Naullage et al. [71] reported the hydrate-C12 interfacial free energy of 76.9 ± 2 mN/m. Contrarily, we estimated the surface tensions of hydrate-C12 and hydratewater interfaces at 42.1 ± 2.2 and 1.65 ± 0.01 mN/m, respectively, directly using GROMACS-LS [61][62][63] (see Fig. 3C). Our calculations yield similar values for the hydrate-C12 (42.1 ± 2.2 mN/m) and hydrate-C5 (42.8 ± 2.6 mN/m) interfacial free energies at similar T (274 -275 K). Our simulation results for the surface free energy of the C1-C2 hydrate-water interface (1.65 ± 0.01 mN/m) obtained at 274 K are comparable to those for the C5 hydrate-water interface (0.32 ± 0.05 mN/m) at 275 K (see Fig. 2C, left). These values are, however, much lower than those used by Naullage et al. (33 ± 4 mN/m) [71], which were obtained by determining equilibrium temperatures of hydrate in hydrate-filled porous media and then employing the Gibbs-Thompson equation to estimate the surface free energy for the hydrate-water interface (using 4 -100 nm pore diameters with an assumed hydrate wetting angle of 0°) [60]. Employing GROMACS-LS [61][62][63], we obtain values for the simulated C12-water interfacial tension (53.7 ± 1.8 mN/ m) that are remarkably consistent with experimental data (53 mN/m) [76]. Applying Young's equation using these values yields equilibrium contact angles of $ 41°for a water droplet on C1-C2 hydrates immersed in C12, which is in excellent agreement with the values obtained directly from our simulations (see Fig. 3B, top). We conclude that the value for the C1 hydrate-C12 surface free energy found by Naullage et al. [71] (76.9 ± 2.0 mN/m) is unrealistically high, compared to our estimates, likely as a consequence of the approximations invoked.
Although our simulation results agree quantitatively with the contact angle measurements reported by Thomas et al. [25], they do not agree with Brown et al. [16], who reported water contact angles of 94.2°± 8.5°on a C5 hydrate immersed in C5. Very recently, Stoner et al. [28] suggested that surface abnormalities as well as hydrate film growth may slow the spreading of a water droplet deposited on the hydrate surface, thereby affecting the observed contact angle (see Fig. 4A). We hypothesize that surface heterogeneities, such as roughness (at various length scales), could affect the measured contact angle, particularly noting that the hydrate surfaces used by Thomas et al. were very smooth, rigid, and flat [25], and noting that Young's equation (see Figs. 2 and 3), is only applicable to smooth and flat surfaces [23,50]. Our hypothesis is supported by the experimental observations by Stoner et al. [28], and also by Spencer et al. [77], who quantified the effect of heterogeneity on contact line pinning. Factors considered included chemical heterogeneity, surface roughness, and interfacial wetting states, which possibly create contact line pinning [78,79]. These authors conducted experiments to measure water contact angles on various rough surfaces with precisely controlled surface chemistry, showing that pinning can shift the contact angles on hydrophilic rough surfaces to more hydrophobic values (reported in detail in ref. [77]). Note that pinning effects induced by roughness become less significant under conditions of strong capillary forces that promote complete wetting [77].
Pinning Effects on Cyclopentane Hydrates. Employing MD, we examine the impact of contact line pinning on the wetting properties for C5 hydrate surfaces immersed in C5. First, we generated rough hydrate surfaces (see Fig. 4B, left) through hydrate growth from water droplets immersed in liquid hydrocarbons spread on flat hydrates (see Fig. 1A/3A for the various configurations). We then performed simulations for water droplets of different sizes on the produced atomically rough C5 hydrate surfaces at 275 K and 0.1 MPa (see Fig. 4B, right). In Fig. 4C, we present 2D density profiles obtained for water droplets of various sizes on the rough C5 hydrate. Extracting the droplet contours and further fitting them with a circular function, we obtained contact angle values (see Fig. 4D). Analysis of the simulation results allows us to distinguish two contact angles for each droplet, i.e., the 'local' and the 'apparent' contact angles. Note that the local contact angle, as defined here, has been referred to as 'intrinsic' contact angle in the literature [81]. On real surfaces, chemically heterogeneous and geometrically rough [82][83][84], these different contact angles could be characterized [85], which can be difficult to assess experimentally. In our nomenclature, we also distinguish local and apparent contact angle from the 'Young' contact angle, which would correspond to the contact angle obtained in an atomically flat and chemically homogeneous substrate. On a rough surface, the local contact angle is formed between the tangent to the local solid surface and the tangent to liquid-liquid interface. The apparent contact angle is instead formed by the macroscopicgeometrical solid surface and the tangent to liquid-liquid interface (see Figure S2 of the SI) [85]. On an atomically smooth and chemically homogeneous surface, in the absence of pinning phenomena, the apparent contact angle is expected to correspond to the local and to the Young contact angles [86]. Following Shuttleworth and Bailey [87] (see Figure S2), we propose that the local angle is the sum of the apparent angle and the angle of the surface slope at the liquid-solid contact point [81].
Within the size limits accessible to our simulations, we observe that the effect of the contact line pinning created by surface heterogeneity becomes insignificant when the droplet volume increases from 28 to 68 nm 3 . This is probably because the base radius of the droplet increases above the length scale controlled by the local roughness of the surfaces considered in our simulations (see Figure S3 of the SI). After a certain droplet size (68 nm 3 ), we observe a sharp increase in the observed contact angle, which is likely due to pronounced contact line pinning effects, which maintain the constant base of the droplet despite the increased droplet volume, yielding a local contact angle as high as $ 75°, which is comparable to the experimental data reported by Brown et al. [16] (94.2°± 8.5°) and by Stoner et al. [28] (94°± 17°). Although the length scales of the two approaches are different, we find this agreement encouraging. Note that the estimated characteristic length of the rough feature (defined as the ratio of the height of roughness to the base radius of the droplet) obtained from experiments ($0.083, estimated from Fig. 4A) and this simulation study (0.085 -0.096) are comparable, justifying the credibility of our comparison.
Pinning Effects on Methane-Ethane Hydrates. We also investigate the influence of contact line pinning on wetting properties on rough C1-C2 hydrates immersed in C12 at 274 K and 3.45 MPa (see Fig. 5A). We implemented simulation and analysis procedures similar to those employed to study water droplets on C5 hydrates. In Fig. 5B, the results of 2D density profiles for the simulated water droplets of various sizes on the rough C1-C2 hydrate are shown. We obtained local and apparent contact angles (see Fig. 5C) by fitting the extracted droplet contours using a circular function. In contrast to the results shown for water droplets on C5 hydrates (see Fig. 4D), the local/apparent contact angles on rough C1-C2 hydrates immersed in C12 increase significantly (from 30°to 45°) when the droplet volume increases from 28 to 42 nm 3 . The pinning effects appear to be more pronounced for small droplets on rough C1-C2 hydrates in C12 compared to the ones on rough C5 hydrates (see Figure S3). This difference could be attributed to the lower water wettability of the C1-C2 hydrates in C12 (h C1-C2inC12 (34 -41°) > h C5inC5 (21 -23°), see Figs. 1, 2 and 3). Recently, when studying the variation of contact angle due to pinning forces, Ozcelik et al. [80] observed that pinning effects at a microscopic level become predominant with the increase of surface heterogeneity and lower surface wettability. When the droplet volume continues to increase, the contact angle increases slightly on rough C1-C2 hydrates, leading to local contact angles of $ 78°, consistent with those measured experimentally in the high-pressure MMF apparatus (85°± 5°) [88] under similar conditions (Fig. 5D).
Contact-line Pinning Force. The results in Figs. 4 and 5 are consistent with the observation that a pinned contact line maintains the base of the droplet constant as the droplet volume increases, which ultimately yields an increase in the measured contact angle. The so-called pinning force is responsible for an energy barrier that the contact line must overcome to move from one metastable state to another [80]. To quantify the pinning effects, we calculate the pinning force F pin on the three-phase contact line using the following expression [80,89,90]: In Eq. (3), c f1,2 is the water-hydrocarbon interfacial tension, h 1 is the contact angle on an atomically flat and chemically homogeneous surface without contact line pinning ('Young' contact angle) and h is the local contact angle measured from simulations (or experiments). In Fig. 6, left, we show the results of pinning force as a function of the droplet volume on the rough C5 hydrates immersed in C5 (blue) at 275 K and 0.1 MPa and on the rough C1-C2 hydrates immersed in C12 (red) at 274 K and 3.45 MPa. We observe that the pinning force on the contact line of water-C5-C5 hydrate barely changes for droplet volumes up to $ 70 nm 3 but quickly grows when further increasing the droplet volume, while the pinning force on the water-C12-C1-C2 hydrate contact line initially increases rapidly, and then grows only slightly when the droplet volume continues to increase. This suggests that the interactions between droplet fluid and the surrounding liquid strongly influence the relation between pinning force and droplet size when (as in our simulations) the rough surface morphologies are similar.  [28]. Copyright 2021 American Chemical Society. B) Representative simulation snapshots for the configurations for the rough C5 hydrate surface and systems composed of a water droplet on the rough C5 hydrate surface in C5 at 275 K and 0.1 MPa. The color code is analogous to that used in Fig. 1A. C) 2D density profiles of water droplets of different size on the rough C5 hydrate surface. The color bar shows water density in the units of 1/nm 3 . D) Variation of the apparent and local contact angle as a function of droplet volume on the rough C5 hydrate surface. The droplet volume is determined from the number of water molecules, following the algorithm described elsewhere [80]. Lines between data points are guides to the eye. In our nomenclature, the local contact angle is formed between the tangent to liquid-liquid interface and the tangent to the local solid surface. The apparent contact angle is instead formed by the tangent to the liquid-liquid interface and the macroscopic-geometrical solid surface.  Fig. 3A. B) 2D density profiles of water droplets of increasing size on the rough C1-C2 hydrate. The color bar shows water density in the units of 1/nm 3 . C) Variation of the apparent and local contact angle as a function of droplet volume on the rough C1-C2 hydrate surface (see Figure S2 of the SI for schematic representation of apparent vs. local contact angles). D) Contact angle measurement in high-pressure MMF, where a water droplet was placed on the C1-C2 hydrate surface immersed in C12 at 274 K and 3.45 MPa, reported by Hu. Adapted with permission from Ref. [88]. Fig. 6. Pinning force as a function of droplet volume (left) and as function of the ratio between the free-standing radius of the given cylindrical droplet from its volume and the size of roughness (right). The results were obtained for the water droplets on a rough C5 hydrate immersed in C5 at 275 K and 0.1 MPa (blue) and on a rough C1-C2 hydrate immersed in C12 at 274 K and 3.45 MPa (red). Lines are guides to the eye. We notice that the pinning forces obtained for both simulated systems are comparable when the droplet volume reaches $ 82 nm 3 . When the droplet volume increases further, eventually the droplets completely wet the rough hydrate surfaces due to the significant increase of capillary force to values much greater than the pinning force [77,91]. It is expected that droplets that are much smaller than the roughness length scale would fully wet the hydrate surfaces between two pinning spots, yielding water contact angles similar to those found on flat hydrates. It is perhaps interesting to point out that our results show a linear correlation between the pinning force and the ratio between the freestanding radius of the cylindrical droplet (independent of the droplet's wetting shape) and the extent of roughness (estimated as the distance between two pinning points) for both rough C5 and C1-C2 hydrate surfaces (see Fig. 6, right). This linear correlation is consistent with previous studies [80,90]. It should be noted that the ratio between the free-standing radius of the cylindrical droplet and the extent of roughness eventually yields a parameter that could be used to tune the pinning force through the design of appropriate surface features [90]. Work of Adhesion. One of the most useful quantities to characterize wetting properties for a solid surface is the work of adhesion, often defined for a liquid droplet on a solid substrate, W 12S . W 12S is described as the energy per unit area required to separate the droplet from the solid. It is determined from the Young contact angle h 1 and water-hydrocarbon interfacial tension c f1,2 using the Young-Dupré equation [92]: With the known values of the water Young contact angle on C5 hydrates immersed in C5 (22.3°) and the C5-water interfacial tension (46.0 ± 2.4 mN/m), the work of adhesion was found to be $ 88.6 mN/m. Similarly, the work of adhesion for the water droplet on the C1-C2 hydrate surfaces immersed in C12 was found to be $ 98.2 mN/m using the water contact angle and interfacial tension for C12 systems (see Fig. 3B and 3C). In Fig. 7, we compare the work of adhesion calculated from our simulations (yellow columns) to experimental results. For example, Chen et al. [93] recently reported that the equilibrium attractive forces measured experimentally between a water droplet (D d = 2.00 ± 0.05 mm) and a C5 hydrate particle (D h = 4.00 ± 0.22 mm) immersed in C5 are of 0.23 ± 0.01 mN at 275 ± 0.5 K, yielding a work of adhesion of $ 86.3 mN/m (see Fig. 7, blue), which is in excellent agreement with our simulation results (88.6 mN/m). To test applicability and reliability of the method, we constructed additional simulation systems composed of the droplet on C5 hydrates immersed in (50 %C5 + 50 %C7) mixtures at 275 K. We obtained a Young contact angle of $ 30°and an interfacial tension of $ 42.7 ± 2.1 mN/m (close to the experimental value of $ 42.2 mN/m obtained for this interface) [94], leading to a work of adhesion of $ 79.7 mN/m, which agrees well with experimental data reported by Chen et al. [93] (78.8 mN/m). This suggests that the direct simulation of water contact angle on an atomically smooth, chemically homogeneous hydrate surface immersed in a liquid, combined with the surface free energy estimated for water-surrounding liquid interfaces could be used to predict, rather accurately, experimental adhesion force measurements.

Conclusions
The literature shows discrepancies on contact angles measured for water on hydrate particles [16,25,28,88]. Such discrepancies hamper technological progress, specifically in flow assurance. We hypothesize that surface heterogeneity is responsible for the wide variation of contact angles reported. To test this hypothesis at the molecular-level, we implement an arsenal of computational materials chemistry methods. We compute directly contact angles and we estimate surface free energies of solid-liquid and liquid-liquid interfaces. The simulation results show excellent agreement with macroscopic experimental observations for both atomically smooth and rough hydrate surfaces under similar temperature and pressure conditions. We conclude that contact line pinning, prompted by surface heterogeneities such as roughness, is responsible for hindering the spreading of the water droplet on the hydrate surface, leading to contact angles larger than those expected based on thermodynamic arguments alone (i.e., the 'Young' contact angle in our nomenclature). This interpretation allows us to reconcile the discrepancy between experimental contact angle values reported as low as 20°, and as high as 80-90°in the literature [16,25,28,88]. We found that our estimates for the work of adhesion are consistent with experimental measurements [93], suggesting that simulations such as those presented here could be useful for the estimation of adhesion force interactions between hydrate particles and water droplets. To our knowledge, this is the first determination of the hydrate-water adhesive forces achieved by employing molecular simulations under realistic conditions, as most previous experimental studies were performed using hydrates stable at atmospheric pressure, such as cyclopentane [93,[95][96][97]. Our work contributes to closing the gap between experimental and realistic models of clathrate hydrates. Because understanding the wetting characteristics of hydrates will provide insights into their formation, growth, and agglomeration, our analysis will be beneficial to quantify the impact of additives on hydrate-forming systems in future studies. For example, it could be of interest to conduct a molecular level study on the effects of various additives adsorbed on hydrate surfaces on wetting properties of relevance to flow assurance.   7. Work of adhesion between a water droplet and a hydrate particle immersed in liquid hydrocarbons. The results were obtained for the systems composed of the water droplet and C5 hydrate particles in either pure C5 or C5-C7 mixture at 275 K and 0.1 MPa, and of the water droplet and the C1-C2 hydrate particle in C12 at 274 K and 3.45 MPa.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.