Role of Structural Rigidity and Collective Behaviour in the Molecular Design of Gas Hydrates Anti-Agglomerants

Antiagglomerants (AAs) are surface active molecules widely used in the rubber and petroleum industry, among others. In the petroleum industry, it is believed that AAs strongly adsorb to the surface of hydrate particles to prevent the growth of clathrate hydrate within oil pipelines. Small changes in their molecular structures can strongly affect the thermodynamic and kinetic stability of the system as a whole. Here we employ molecular dynamics simulations to study the interplay between the modification of the molecular structure, rigidity and collective effects of AAs designed to prevent hydrate agglomeration in the conditions encountered in rocking cell experiments. The AAs are surface-active compounds with a complex hydrophilic head and three hydrophobic tails whose structural rigidity is enhanced with the attachement of a simple aromatic group. We observe that the aromatic group can positively or negatively affect the performance of the AAs, depending on its location along the hydrophobic tail. Our approach is based on first quantifying the molecular mechanisms responsible for the macroscopic performance. Although the mechanisms at play depend on the application, the methodology implemented could be applicable to other high-tech industries, where the agglomeration of small particles must be controlled.

Gas hydrates, also known as clathrate hydrates, are ice-like inclusion compounds consisting of polyhedral hydrogen-bonded water cages stabilized by guest gas molecules. [1][2][3] They are formed under high-pressure and low-temperature conditions such as those found in deep oceans and pipelines. [4] Clathrate hydrates are relevant in a variety of scientific and industrial contexts, including climate change modeling, [5] carbon dioxide sequestration, [6] hydrocarbon extraction, [7] hydrogen and natural gas storage, [7][8][9] separation and refrigeration technologies, [10] marine biology, [11] and planetary surface chemistry. [12] Of particular interest are the hydrocarbon hydrates that can form blockages in oil and gas pipelines. This phenomenon can severely affect the safety of pipeline flow assurance, potentially leading to large negative environmental consequences. [1,3,13,14] Among the range of approaches routinely used to manage gas hydrates is the stabilization of hydrate particles in hydrocarbon dispersions. Specifically designed surfactants, known as anti-agglomerants (AAs), are optimized to prevent hydrate plug formation in flow assurance. [2,[15][16][17] AAs are believed to adsorb on hydrate particles by their hydrophilic head groups, while the AAs tail groups are soluble in the hydrocarbon phase. The hydrate particles are expected to be covered by a film of AAs and oil, making them repel each other and remain dispersed. [18][19][20][21][22] Laboratory and field observations alike show that many phenomena determine the performance of AAs. For instance, small changes in the molecular structure of the surfactants can strongly affect their ability to prevent hydrate plugs formation. [22] Of particular interest in the present work is the fact that the stability of water-in-oil emulsions is * Corresponding author: francois.sicard@free.fr. dependent on the rigidity of the interfacial film, which in turn can be determined by the chemical structure and the collective effect of the AAs. [19,21,22] Here, we employ molecular dynamics (MD) simulations and enhanced sampling techniques to study the molecular-level properties along with thermodynamic and kinetic information of AAs specifically designed to prevent hydrate agglomerations in the conditions encountered in rocking cell experiments. The AAs considered are based on a compound, recently developed, which has shown good laboratory performance in preventing hydrate formation in light oils. [19,21,22] The molecular structure of the AA tail groups was computationally modified to improve the rigidity of the compound, while maintaining the ability of the polar head group to adsorb on the hydrate surface. We analyse the behaviour of the newly designed AA at the hydrate-oil interfaces (Figure 1), in terms of density profiles and preferential orientation, and quantify the thermodynamic and kinetic properties associated with the transport of the methane molecules across the AA film. Details regarding simulation models and algorithms are reported in the Supporting Information (SI). We observe that, at AA densities similar to those considered in rocking cell experiments, the attachement of the simple aromatic can positively or negatively affect the performance of the AAs depending on its location along the hydrophobic tail. In particular, collective and synergistic effects between the newly synthetic AAs and the liquid hydrocarbon can yield a free energy (FE) barrier for methane transport through the surfactant layer, whose height is significantly higher than the one measured for synthetic AAs routinely used in the oil and gas industry. Based on our interpretation of the simulation results, this suggests that the engineered AAs could have much higher performance in preventing the Green spheres represent methane molecules either in the sII methane hydrate or the hydrocarbon phase. Blue dots lines represent water molecules in the hydrate substrate. Silver lines represent n-dodecane molecules, either in the bulk or trapped within the AA layer. Yellow, red, blue, white, and cyan spheres represent chloride ions, oxygen, nitrogen, hydrogen, and carbon atoms in AA molecules, respectively. Only half of the simulation box is shown here for clarity.
formation of hydrate plugs than existing AAs. Experimental verification is required to test our expectations.
The two AA molecular structures considered throughout this work are shown in Fig. 2. The short tail R 2 is a linear hydrocarbon chain of four carbon atoms (n-butyl). The headgroup including boh amide and tertiary ammonium cation group is known for its ability to adsorb on the hydrate surface [19][20][21][22]. The two long R 1 tails are composed of a linear hydrocarbon chain of twelve carbon atoms with an aromatic ring (benzene) positioned either at the bottom (Fig. 2b) or in the middle (Fig. 2c) of the chain. The aromatic ring is intended to rigidify the structure of the AA tails thus providing significantly stronger mechanical strength and thermal stability. [23] The simulations were run at temperature and pressure maintained at 277 K and 20 MPa, similar to those encountered in laboratory experiments [19] (see details in the SI). Thus, the conditions chosen in this study are well within the gas hydrate stability zone and correspond to subcooled systems [20].
Visual Observation of Simulation Snapshots. We first considered in Fig. 3 the AA structure given in Fig. 2b, where the aromatic ring is positioned at the bottom of the long tail R 1 . Representative simulation snap-shots of the AA configurations adsorbed on the hydrate surface are shown in Fig. 3b at low (0.33 molecule.nm −2 ) and high (0.67 molecule.nm −2 ) surface densities. The snapshots are taken at the end of our simulations. We observe that the headgroups of the AA adsorb on the hydrate surface, as noticed previously [19][20][21][22]. The long tails of AAs are instead more likely to extend towards the alkane bulk phase. The snapshots shown in Fig. 3b suggest that different AA densities yield differences in the thin-film structure. The increase of the AA surface density yields a transition between disordered to ordered orientation of the AA long tails with a thin film within which the long tails of AAs and n-dodecane align parallel to each other and orient perpendicularly to the hydrate surface. This transition comes with the configurational change of the n-dodecane molecules from a Gauche conformation in the bulk phase to a nearly all-trans conformation within the AA film.
Although changing the position of the aromatic ring from the bottom to the middle in the long tails did not seem to affect the ability of the AA to adsorb at the hydrate surface, it significantly affected the ordering of the thin-film structure, even at the relatively high surface density (0.67 molecule.nm −2 ) studied here, as shown in Fig. 4a. The lack of ordering of the AA film also comes with n-dodecane molecules in a Gauche conformation either in the bulk phase or within the surfactant layer.
Density Profiles. To quantify the influence of AAs adsorbed on the hydrate surface on the distribution of methane in the system, we calculated the mass density profiles along the z-axis of the simulation box  for increasing AA surface density. In the following Z and Z box represent the projection of the Cartesian position of the methane molecules and the size of the simulation box along the z-axis, respectively. Figs. 3c and 4b show a high and constant density between Z/Z box = −0.1 and Z/Z box = 0.1, which corresponds to the methane molecules trapped in the hydrate cages. For Z/Z box > 0.3 and Z/Z box < −0.3 the results show a uniform density representative of the fluid hydrocarbon phase. The methane density profile in the thin region between the layer of AA headgroups and the bulk liquid hydrocarbon phase shows a pronounced dependence on the AA type and surface density. At low AA surface densities, the dendity of methane near the hydrate surface for systems containing the aromatic ring either at the bottom (Fig. 3c) or in the middle (Fig. 4b) of the AA long tails is similar to those found in the bulk. When the aromatic ring is positioned at the bottom of the AA long tails (Fig. 3c), the results show a pronounced depletion of methane at the interface (Z/Z box < 0.25) as the AAs surface density increases to 0.67 molecule.nm −2 , with the density profile being nearly 0. Combined with a visual observation of the simulation snapshots, these results suggests that the ordered layer successfully expels methane from the interfacial region. This phenomemom is similar to the one discussed in previous work for synthetic AA used in the gas and oil industry [19,21], which can be explained by the collective and the synergistic effects between the AAs and the hydrocarbon phase.
AA Orientation. To quantify the orientation of AAs at the interface, we considered the orientational angle, φ, formed between each tail and the direction perpendicular to the hydrate surface. We calculated the probability distribution of this angle as well as that of the conformational angle, θ, between the two long tails of one AA molecule. In Figs. 3d and 4c, we report the probability distribution of the orientational and configurational angles for the two AA structures represented in Figs. 2b and 2c at various surface densities. At low surface density, the orientational angle shows a wide probability distributions, from 0 to 90 • and above, irrespectively of the position of the aromatic ring in the hydrophobic long tail. Similarly, the conformational angle does not show preferential values at low surface coverage for either AA. These results suggest that the AAs are rather disordered at these conditions. When the AA surface density increases, the results show significant variations. While the results obtained for the compound with the aromatic ring positioned in the middle of the long tails do not show substantial changes compared to those obtained at low surfae density (Fig. 4c), the results obtained for the compound with the aromatic ring at the bottom show pronounced order. As shown in Fig. 3d, the orientational distribution shows a narrow peak at Φ ≈ 20 • when the surface densiy increases to 0.67 molecule.nm −2 , suggesting that the AA tails become almost perpendicular to the hydrate surface. At the same surface density, the conformational distribution shows pronounced peak at θ ≈ 10 • , suggesting that the AAs maintain their long tails almost parallel to each other at this conditions. As the only difference between the two AA structures simulated is the position of the aromatic ring in the hydrophobic long tail, the differences highlighted in Fig. 3 and Fig. 4 are likely due to collective effects and preferential interactions between the AA long tail and the hydrocarbon molecules in the fluid phase due to attractive chain-chain lateral van der Waals interactions.
Thermodynamic properties. We quantify the thermodynamic characteristics of the AA layer with the analysis of the free energy profile (FEP) associated with the transport of methane molecules. Following previous work [21], we focus our analysis on the diffusion of methane through an interfacial region made up of the largest cluster of hydrocarbons (characteristic size ≈ 20Å), after the free methane molecules were expelled from the interfacial region.
In this situation, the methane molecules interact mainly with hydrocarbon molecules trapped in the AA film. In Fig. 5a, it is shown the FEP obtained within the US/ABMD framework employed by Sicard and coworkers, [21,24] and using the dynamic histogram analysis method (DHAM) [25] plotted along the reduced units, Z/Z box (see details in the SI). Uncertainties were determined by dividing the data into four equal sections, determining the profiles independently, and calculating the standard error. We measured a difference in FE between the global (Z/Z box ≈ 0.26) and local (Z/Z box ≈ 0.36) minima ∆F 0 ≈ 8.5 kJ/mol. These two basins are well separated by activation energies associated with methane capture and escape, ∆F C ≈ 35 kJ/mol and ∆F E ≈ 26.5 kK/mol, respectively. As shown in Fig. 5c, the free methane molecule is initially in the bulk hydrocarbon phase, above the AA layer. When it comes closer to the interface, it is first trapped in a local FE minimum (Z/Z box ≈ 0.26). This minimum corresponds to a transition region between oil molecules isotropically oriented in the bulk and oil molecules parallel to the AA tails. The methane molecule then enters the interfacial film. As the methane travels farther across the interfacial layer, an energy barrier arises as the oil molecules are displaced from the methane pathway. The transport proceeds until one oil molecule cannot be pushed farther down (0.26 ≤ Z/Z box ≤ 0.31). Under this conditions, the oil molecule bends, eventually forming a cage surrounding the methane molecule (Z/Z box ≈ 0.31). This corresponds to the high-energy transition region in the FEP. Once the methane molecule overcomes this transition state, it is pushed down underneath the AA layer. The methane molecule then reaches the local minimum corresponding to the water layer between the AA layer and the hydrate (Z/Z box ≈ 0.36).
Kinetic properties. To complement the thermodynamic analysis, we estimated the position-dependent diffusion profile which provides molecular understanding of the transport of solute across three-dimensional heterogeneous media. [26][27][28] In this system, the variation of the solute diffusivity can be impacted by variation of the frictional environment as the solute moves from bulk hydrocarbon through interface, and into the water layer. We extended the standard scope of the US framework considering the method originated by Berne and co-workers [29] and elaborated by Hummer, [30] where the diffusion coefficient is calculated from the position autocorrelation function (PACF) obtained from harmonically restrained simulations In Eq. S2, z k is the average of the RC in the US window k, var(z) = z 2 − z 2 is its variance, and C zz (t) = δz(0)δz(t) the PACF calculated directly from the time series. In Fig. 5b, it is shown the position-dependent diffusion profile along the Z direction of methane across the AAs layer. Uncertainties were determined by dividing the data into four equal sections, determining the profiles independently, and calculating the standard error. The system shows a diffusion profile with two distinct plateaus located at positions D Z/Z box < 0.31 ≈ 0.6 × 10 −10 m 2 .s −1 and D Z/Z box < 0.31 ≈ 2.2 × 10 −10 m 2 .s −1 on both sides of the transition state Z/Z box ≈ 0.31. When methane enters the interfacial layer through two oil molecules, the effective diffusion coefficient is similar to the one measured experimentally in bulk hydrocarbons [31] (≈ 5.10 −11 m 2 .s −1 at 277 K and 20 MPa). As the methane goes farther across the interfacial layer, the hydrocarbon molecules surrounding it start to push it down underneath the AA layer, increasing the effective diffusion coefficient by almost an order of magnitude.
In conclusion, the extensive simulations discussed above highlight the role of the structural rigidity and collective behaviour in the molecular design of new gas hydrate AAs. Based on our simulation results, we studied the molecular-level properties of a newly synthetic AA specifically designed in silico to prevent hydrate plugs formation in rocking cell experiments. We studied the thermodynamic and kinetic characteristics of the system and quantified accurately the FEP experienced by one methane molecule travelling across the interfacial film along with the position-dependent diffusion coefficient, using a combination of MD simulations and enhanced sampling techniques. We showed that the FE barrier associated with methane transport across the AAs film is due to collective and synergistic effects between the AAs and the liquid hydrocarbon trapped in the interfacial layer.
Interestingly, the transport properties due to the AAs specifically designed in this work can be compared with those associated with synthetic AAs, which show good laboratory performance in preventing hydrate formation in light oils [19,21]. In the latter, the calculated FE of activation for methane capture was significantly lower (≈ 15 kJ/mol). It suggests that tuning the rigidity of AAs at the structural level, while preserving their ability to interact with the hydrocarbon molecules, can significantly improve their collective performance in preventing hydrate formation.
The results presented here follow from the hypothesis that transport trough the AAs film determines AAs performance in flow assurance. The computational methodology developed could be useful for a variety of hightech technologies in the petroleum, rubber latex, ink, and paint and coatings industries, where the agglomer-ation and flocculation of small particles must be controlled [32][33][34]. Although the mechanisms by which AAs are effective depends on the application of interest, several research and industrial groups focus on the interplay between the design and performance of AAs at the microscopic scale [35][36][37][38][39]. Accounting for the interplay between the structural rigidity and the collective effects of new synthetic AAs could allow us to infer their use in practical applications, which addresses current industrial needs.
Acknowledgments FS thanks Jhoan Toro-Mendoza, Denes Berta, and Tai Bui for useful discussions. Via our membership of the UKs HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). Financial support was graciously provided by the EPSRC to AS under grant number EP/T004282/1. The project was initiated under the support of EPSRC, via grant number EP/N007123/1.

Role of Structural Rigidity and Collective Behaviour in the Molecular Design of Gas
Hydrates Anti-Agglomerants

Molecular Dynamics (MD) simulations
Molecular dynamics (MD) simulations were performed with the GROMACS software package, version 5.1.4 [40] using the TIP4P/Ice water model [41]. Biased simulations were performed using version 2.3 of the plugin for FE calculation, PLUMED [42]. The TIP4P/Ice model has been successfully implemented to study hydrate nucleation and growth [43,44] and to investigate the performance of potential hydrate inhibitors [45]. This model yields an equilibrium temperature for the formation of gas hydrates at high pressure close to experimental values [46]. Methane and n-dodecane were represented within the united-atom version of the TraPPE-UA force field [47]. AAs were modeled using the general Amber force field (GAFF) [48], which is often implemented for modeling organic and pharmaceutical molecules containing H, C, N, O, S, P, and halogens. Atomic charges were calculated with the AM1-BCC method employed in Antechamber from the Amber 14 suite [49]. The chloride counterions (Cl − ) were modeled as charged Lennard-Jones (LJ) spheres with the potential parameters taken from Dang [50], without polarizability. The sII hydrates were considered to be the solid substrate, and they were not allowed to vibrate in this work. AAs, chloride counterions, n-dodecane, and methane composed the liquid phase. Dispersive and electrostatic interactions were modeled by the 12 − 6 LJ and Coulombic potentials, respectively. The Lorentz-Berthelot mixing rules [51,52] were applied to determine the LJ parameters for unlike interactions from the parameters of the pure components. The distance cutoff for all non-bonded interactions was set to 1.4 nm. Long-range corrections to the electrostatic interactions were described using the particle mesh Ewald (PME) method [53][54][55] with a Fourier grid spacing of 0.12 nm, a tolerance of 10 −5 , and fourth-order interpolation. Periodic boundary conditions were applied in three dimensions for all simulations.

Unbiased MD simulations
To construct the initial configurations, we followed the procedure described in previous work [19,21]. The sII methane hydrate was chosen to represent features of the experimental system considered, in which a small amount of gases other than methane is present. The underlying assumption is that the host gas does not affect the properties of the AAs film, which is the subject matter of this investigation. One unit cell of sII methane hydrates was adapted from the study of Takeuchi et al. [56]. The sII methane hydrate unit cell was replicated three times in the X and Y directions (5.193 nm) and two times in the Z direction (3.462 nm). It was then flanked by a thin liquid water film of approximately 0.5 nm on both sides along the Z direction, which represents the quasi-liquid interfacial layer identified in the experiments of Aman and coworkers [57]. The desired number of AA molecules was arranged near both sides of the hydrate substrate. The chloride counterions (Cl − ) were placed next to the AA headgroups. The n-dodecane and methane molecules were placed within the remainder of the simulation box. The time step used in all the simulations was 0.001 ps, and the list of neighbors was updated every 0.01 ps with the grid method and a cutoff radius of 1.4 nm.
The initial configuration was first relaxed using the "steepest descent minimization" algorithm to remove high-energy configurations, which might be related to steric hindrance between AAs and the hydrocarbon phase. During this step, the gas hydrate structure remained unaltered. Subsequently, to minimize the possibility that the initial configuration biased the simulation results, an N V T temperature-annealing procedure, as implemented in GROMACS [40], was conducted. The algorithm linearly decreased the system temperature from 1000 K to 277 K in 500 ps. In these simulations, the hydrate substrate and chloride ions were kept fixed in position. To relax the structure of n-dodecane and AAs, a N V T simulation was conducted at 277 K for 2 ns using the Berendsen thermostat [58], with the sII hydrate structure kept fixed in position. The equilibration phase was then conducted within the isobaric-isothermal (N P T ) ensemble under thermodynamic conditions favorable for hydrate formation (T = 277 K and P = 20 MPa) to equilibrate the fluid density. During the NPT simulation, all molecules in the system were allowed to move, including water and methane molecules in the hydrate substrate. The pressure coupling was applied only along the Z direction of the simulation box, which allowed the X and Y dimensions to be maintained constant. Temperature and pressure were maintained at 277 K and 20 MPa, respectively using the Berendsen thermostat and barostat [58] for 5 ns. This is considered the most efficient algorithm to scale simulation boxes at the beginning of a simulation [59]. We then switched to the Nose-Hoover thermostat [60] and the Parrinello-Rahman barostat [61] for 100 ns, which are considered more thermodynamically consistent algorithms [59]. This numerical protocol allowed the AAs to assemble and orient to form the interfacial layer depicted in Figure 1 in the main text. [19,21] The system was then equilibrated for 3 ns in N V T conditions coupling with the v-rescale thermostat [62] (T = 277 K, τ T = 0.1 ps). To define the position of the methane molecules in the simulation box with respect to the sII hydrate structure, the simulation was continued in N V T conditions holding in place the methane molecule enclathrated into the water cages defining the sII hydrate structure.

Biased MD simulations
The phenomenon of interest (i.e. the transport of methane across the interfacial layer) occurs on time scales that are orders of magnitude longer than the accessible time that can be currently simulated with classical MD simulations. A variety of methods, referred to as enhanced sampling techniques [63][64][65][66][67], can be implemented to overcome this limitation. These methods accelerate rare events and are based on constrained MD. In the present work, we used the numerical method employed by Sicard and coworkers, [21,24], which combined the adiabatic biased molecular dynamics (ABMD) [68][69][70][71] and umbrella sampling (US) [72] frameworks. To design the US windows, we used the projection of the Cartesian position of the methane molecules along the Z-direction as a reaction coordinate (RC). The starting configurations for the US simulations were obtained by pulling adiabatically the system along the RC generating 40 windows. Each US window was subsequently run for 2 ns to allow equilibration, followed by additional 16 ns of the production run. To control the accuracy of the sampling with respect to the RC orthogonal degrees of freedom (either X-or Y -direction) we implemented the flat-bottomed potential, [73][74][75] as already described in previous work [21].
Upon completion of the US simulations we obtained the free energy profile (FEP) associated with the transport of the methane molecules across the AA film using the dynamic histogram analysis method (DHAM). [25] This unbiasing method, originally derived by Rosta and Hummer, uses a maximum likelihood estimate of the Markov State Model (MSM) transition probabilities given the observed transition counts during each biased trajectory. To produce an MSM from the enhanced sampling simulations, we first discretized the RC into bins, where the number of bins is chosen sufficiently large to give a finely discretized coordinate but not so large as to give an under-sampling of transitions between bins. Once the bins have been determined, we count the number of observed transitions C k ij (t) between each pair of bins i and j in US simulation k at the chosen lagtime t, as well as the number of times each bin is occupied n k i = j C k ij (t) during each US simulation k. These values then provide the necessary conditional probabilities M ij (t) = P (j, t|i, 0). For biased simulations, where a biasing energy of u k i is applied to state i during simulation k, we computed the unbiased MSM from the biased data as given by (S1) Once the MSM has been constructed from simulation data, the equilibrium probabilities can be calculated as the eigenvector corresponding to eigenvalue 1 of the transition matrix M ij obtained in Eq. S1, that is, as its invariant distribution. To estimate the position-dependent diffusion coefficient associated with the methane transport across the interfacial layer, we extended the standard scope of the US framework considering the method originated by Berne and co-workers [29] and elaborated by Hummer, [30] where the diffusion coefficient is calculated from the position autocorrelation function (PACF) obtained from harmonically restrained simulations In Eq. S2, z k is the average of the RC in the US window k, var(z) = z 2 − z 2 is its variance, and C zz (t) = δz(0)δz(t) the PACF calculated directly from the time series.