Molecular Dynamics Analysis of Anti-Agglomerant ... - ACS Publications

Subscriber access provided by READING UNIV

Article

Molecular Dynamics Analysis of Anti-Agglomerant Surface Adsorption in Natural Gas Hydrates Michael A. Bellucci, Matthew R. Walsh, and Bernhardt L Trout J. Phys. Chem. C, Just Accepted Manuscript • Publication Date (Web): 08 Jan 2018 Downloaded from http://pubs.acs.org on January 8, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics Analysis of Anti-Agglomerant Surface Adsorption in Natural Gas Hydrates Michael A. Bellucci,† Matthew R. Walsh,‡ and Bernhardt L. Trout∗,† Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA, and Chevron Energy Technology Company, Flow Assurance, Houston, Texas 77002, USA E-mail: [email protected]

Abstract We have used molecular dynamics simulations to examine the surface adsorption of a model anti-agglomerant inhibitor (quaternary ammonium salt) binding to a hydrate surface in both aqueous and liquid hydrocarbon phases. From our molecular dynamics simulation data, we were able to identify the preferred binding sites on the (111) crystal face of a methane-propane sII hydrate as well as characterize the equilibrium binding configurations of the inhibitor and their associated binding energies. In the aqueous phase, we observed that the inhibitor proceeds through a two-step surface adsorption mechanism whereas in the liquid hydrocarbon phase, surface adsorption occurs through a single-step mechanism. To characterize the extent of surface adsorption in each liquid phase, we calculated the standard binding free energy using the free energy perturbation method following a double decoupling thermodynamic cycle. We found that the surface adsorption in the liquid hydrocarbon phase is an exergonic process whereas the surface ∗

To whom correspondence should be addressed Massachusetts Institute of Technology ‡ Chevron Corporation †

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

adsorption in the aqueous phase is an endergonic process. Our results demonstrate that the extent of surface adsorption is much larger in the liquid hydrocarbon phase relative to the aqueous phase, and suggest that the inhibitor is less effective in the aqueous phase because the surface adsorption is less favorable. Finally, we examine the effect of the inhibitor on the water structure in the liquid phase and in the hydrate phase, with the results highlighting the difference between the nature of anti-agglomerant/hydrate interactions as compared to kinetic inhibitor/hydrate interactions.

Introduction Natural gas hydrates are non-stoichiometric crystalline compounds in which small hydrocarbons or other gases (e.g. CO2 , H2 S) are trapped within a lattice structure of polyhedral cages by hydrogen bonded water molecules. 1–3 Increasingly, they are of interest to many research groups, industries, and government agencies due to their impact on global climate trends 4–8 and because they have a broad range of technological applications such as carbon dioxide sequestration, gas storage and transportation, and water desalination. 9–13 Moreover, since hydrates naturally form at low temperatures and moderate pressures, vast reservoirs of hydrates exist under the sediments of the ocean floors, resulting in an enormous potential resource of fossil-fuel energy. 13–15 However, in the oil and gas industry, hydrates can cause major flow assurance challenges by blocking transmission pipelines, leading to serious safety and environmental problems as well as significant economic loss. 1,16–18 The failure of the containment dome to capture leaking oil in the 2010 Macondo oil spill in the Gulf of Mexico is a recent example of the serious environmental consequences of hydrate blockage. 19 With industry trends moving toward deepwater oil and gas production, over the past few decades preventing hydrate plug formation in pipelines has become an active area of research. To prevent the formation of hydrate plugs in pipelines via chemical injection, the traditional approach is to use thermodynamic inhibitors, such as methanol or ethylene glycol. 1–3,21 These inhibitors work by shifting the hydrate stability region away from the temperature 2


Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and pressure conditions of the pipeline, thereby reducing the thermodynamic driving force and rendering hydrate formation impossible at fully-dosed conditions. However, large quantities of thermodynamic inhibitors are often required for full dosing, which can exceed 50% by volume in the aqueous phase, and therefore their application can be prohibitively costly for both full and partial dosing. 20 As a representative example, one leading industry expert estimates that the cost for full continuous thermodynamic inhibition of a single oil well (tied into a long subsea pipeline) can exceed one million dollars per day. 21 Another chemical approach to hydrate plug prevention is to use low-dosage hydrate inhibitors (LDHIs), which can prevent hydrate plugs at smaller dosages when compared to thermodynamic inhibitors, i.e. between 0.5-2.5 vol % in the aqueous phase. There are two classes of LDHIs: kinetic hydrate inhibitors (KHIs) and anti-agglomerants (AAs). The KHI class of molecules are thought to delay hydrate nucleation, but many KHIs also bind to hydrate surfaces and interfere with hydrate growth. 1–3 A significant limitation of KHIs is that they are often ineffective in deepwater operations since there are large thermodynamic driving forces for hydrate formation in these conditions. In contrast, the AA class of molecules can bind to hydrate surfaces and change the surface chemistry, morphology, and size of the solids to prevent hydrate particles from agglomerating and forming into large masses. In effect, AAs keep hydrate particles dispersed in produced fluids, allowing transport of hydrates as a slurry. As a result, AAs are less sensitive to thermodynamic driving forces for hydrate formation and are more effective in deepwater operations. The AA class of molecules are generally surfactants that contain a hydrophilic headgroup and a hydrophobic tail group. For example, some of the most effective and prototypical AAs are quaternary amines (QA), which were first developed by Shell in the early nineties. As operational oilfields can range in watercut from 0% to 99%, it is important for AAs to be effective over the entire range of watercut, i.e. the water volume fraction in the condensed phases. Across this wide range of watercuts, hydrates can form at the interfaces between gas and water or between oil and water, and can reside at any interface (including

3



the oil/gas interface). Moreover, they can reside within bulk free oil, bulk free water, oilcontinuous emulsions, water-continuous emulsions, and in many cases hydrates can partition between two or more coexisting phases simultaneously. For a variety of reasons, such as a potentially higher hydrate volume fraction, free water breakout, and enhanced ability for capillary attraction between hydrate particles, many AAs become ineffective when the watercut is above 70%. However, recent work has suggested that certain AAs can be effective when the watercut is above 70%. 22,23 Of note in the literature is the work of Sun and Firoozabadi who report an AA (cocamidopropyl dimethylamine) that is possibly effective at low doses and without an oil-phase. 23 In the absence of an oil phase, they postulate that their AA molecules bind to hydrate surfaces and the long hydrocarbon chain of their molecule covers the hydrate particle surfaces, which in turn prevents the hydrates from agglomerating. Similarly, it has been suggested that adsorption of surfactants at the liquid-solid interface is a crucial component of the mechanism of anti-agglomeration, 24,25 but with limited mechanistic studies it remains unclear what features, geometries, properties and moieties could constitute an effective AA across a range of water cuts. Molecular simulations have contributed much to our understanding of gas hydrates, including their phase behavior, 26–28 nucleation, 29–42 growth, 43,44 and transport properties. 45–47 English and MacElroy recently reviewed the many different molecular-simulation-based hydrate investigations over the last two decades. 48 In addition, a number of studies using molecular simulation have been performed to elucidate the mechanism of kinetic inhibition, 49–57 and the investigation of AA’s using molecular simulation has recently begun in earnest. 58,59 To design more effective AA inhibitors and to elucidate AA-hydrate interaction mechanisms, it is necessary to study these systems on a molecular level to understand the binding mechanisms and energetics, and eventually establish microscopic quantities that can be measured from simulation which are predictive of experimental performance, perhaps allowing for a reduction in trial and error in the innovation process. Furthermore, as mentioned above, hydrates can exist in both oil-dominated and water-dominated environments,

4


Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and so a molecular understanding of AA-hydrate binding is needed in both water and liquid hydrocarbon phases. In order to gain a molecular-level understanding of AA-hydrate interactions, in this work we have utilized molecular dynamics simulations to investigate the surface adsorption of a typical AA molecule (quaternary ammonium salt) at the liquid-hydrate interface in both an aqueous phase (100 % watercut) and a liquid hydrocarbon phase (0 % watercut) at state conditions similar to deepwater operating conditions. Using molecular dynamics, we can directly probe how different liquid phase environments influence the surface adsorption of AA molecules, which is a necessary and perhaps key step in AA effectiveness. To characterize the surface adsorption, we performed molecular dynamics simulations in each of the different liquid phase environments, i.e. for the 100% watercut and 0% watercut systems. From this simulation data, we used order parameters (described below) to analyze our trajectories and identify the polyhedral cage binding sites on the hydrate surface. Using this approach, we were able to characterize the binding configurations of the AA molecule as well as the binding energies to the various sites of the hydrate surface. We then utilized representative equilibrium binding configurations and calculated the binding free energies with free energy perturbation (FEP) following a double decoupling thermodynamic cycle. 60–62 Finally, we analyzed the structural effects of AA binding to the hydrate surface using relevant radial distribution functions for bound and unbound states of the AA molecule.

Methodology Systems. Our choice of AA is n-dodecyl-tri(n-butyl)-ammonium chloride (see Fig. 1), which is a prototypical and effective AA for hydrate systems. 1,3,18 We have performed simulations in each liquid phase environment using the QA salt (QAS) and the QA cation (QAC), i.e. with and without an explicit chloride anion, respectively. Throughout this manuscript we refer to the simulations performed without an explicit chloride anion as the QAC system

5



and to the simulations performed with an explicit chloride anion as the QAS system. When using lattice-sum electrostatics, in the absence of an explicit counterion the system is neutralized with a uniform background charge. 63,64 Utilizing a neutralizing background charge can be advantageous because performing molecular dynamics simulations with explicit counterions in a finite-size simulation can cause a number of sampling problems. This is due to their slow equilibration and a low demixing rate between the ions that can lead to hysteresis and results that can have a strong dependence on the initial placement of the ions. 65–68 However, for systems with a nonhomogenous dielectric constant, such as a hydrate in a liquid hydrocarbon phase, the use of a uniform background charge is unphysical, because in reality the counterion would mostly be found in the portion of the system with a higher dielectric. Furthermore, Hub et. al. demonstrated that the use of a uniform background charge in systems with a nonhomogenous dielectric constant can lead to varying artifacts in the energetics. 69 For these reasons, we have elected to perform our simulations with and without an explicit counterion for comparison. The sII hydrate utilized in our simulations is a particular polymorph of natural gas hydrates that is the predominant polymorph found in oil and gas pipeline blockages. It is comprised of two different polyhedral cages: a smaller pentagonal dodecahedral cage, commonly denoted as the 512 cage, and the larger hexakaidecahedral cage, commonly denoted as the 512 64 cage. For both systems, initially a crystalline sII hydrate made of 3×3×3 unit cells was generated from sII unit cell data. 1 We prepared the crystal structure such that the (111) face of the hydrate was exposed in the z direction by removing a cubic cell from the 3×3×3 supercell using an in-house program. The cubic cell had orthogonal lattice vectors oriented along the [-110], [-1-12], and [111] directions with respect to the original unit cell. The cubic cell had periodic dimensions of 12.23×7.06×10.00 Å and the hydrate slab was built as a 4×6×4 supercell from this cubic cell. The (111) face was chosen because it has previously been identified as crucial in the growth process of hydrates and has been used in number of KHI surface adsorption, growth, and inhibition studies, both computational and

6


Page 6 of 34

Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


experimental. 1,49,51,52 In addition, since the (111) face of the hydrate is the slowest-growing face, 1 it should have the largest surface area, and therefore be the most relevant crystal face for adsorption of an anti-agglomerant molecule on the crystal surface. Once the initial hydrate crystals were created, each of the liquid phases were constructed such that they had the same dimensions of the hydrate crystal in the x-y plane and had a width of 40 Å in the z direction. The liquid hydrocarbon phase utilized in our simulations consisted of a mixture of decane, methane, and propane. A single QAC or QAS was solvated in the center of each of the liquid phases and the composition of methane and propane in the initial liquid phases were prepared according to solubility data calculated from the SRK Peneloux equation of state. 70,71 Subsequently, each liquid phase was placed on top of a hydrate crystal in the z direction, leading to a distance of about 20 Å between the AA and the hydrate surface in the z direction. During the equilibration of the systems, a portion of the cages near the interface dissociated resulting in excess methane and propane in the liquid phases. Consequently, we removed these excess methane and propane molecules and equilibrated the system once again to avoid supersaturated bulk solutions. The final configuration of the aqueous system had dimensions of 49 Å × 43 Å × 80 Å and consisted of 1 AA, 205 methane molecules, 90 propane molecules, and 4939 water molecules. The final configuration of the liquid hydrocarbon system had dimensions of 49 Å × 43 Å × 83 Å and consisted of 1 AA, 246 methane molecules, 135 propane molecules, 227 decane molecules, and 2472 water molecules. After equilibration of the systems, the composition of the hydrate slab (including the interfacial region) in both phases was 8.2 mol % methane, 4.8 mol % propane, and 87.0 mol % water, whereas the composition of the purely crystalline portion of the hydrate slab (excluding the interfacial region) was 10 mol % methane, 5 mol % propane, and 85 mol %water. Force Fields. The TIP4P/2005 water model was used as the force field for water molecules in our simulations. This force field has been shown to accurately reproduce a number of liquid water properties as well as solid state properties. 27,28 In addition, it has

7



Figure 1: Structure of n-dodecyl-tri(n-butyl)-ammonium chloride anti-agglomerant. been shown to accurately reproduce the density and stability of hydrates for temperatures above 150 K as well as the density of water in the range of 270-370 K. 28 The general Amber force field (GAFF) was used for the AA, decane molecules, and chloride anion. GAFF has been used previously for the simulation of quaternary ammonium cations. 72,73 The atomic partial charges for the AA and decane were computed using the Merz-Singh-Kollman scheme at the B3LYP/6-31G(d,p) level of theory. In addition, constraints were implemented in the electrostatic potential fitting such that the atomic partial charges reproduced the dipole moment of the molecules. For methane and propane, we used all-atom potentials where the intramolecular interactions of these molecules were modeled with GAFF. However, we modeled the methane-water intermolecular interactions using the partial charges and LennardJones parameters from Ref. [ 49 ] and the propane-water intermolecular interactions using the Lennard-Jones parameters from Ref. [ 74 ]. The Lorentz-Bertholet combining rules were used for the Lennard-Jones parameters between different types of atoms. These intermolecular potentials have been parameterized using high level ab initio methods taking into account hydrate guest-host interactions and they have been shown to reliably reproduce experimental data for hydrate systems. A table of these guest-host interaction parameters is given in the Supporting Information. Molecular Dynamics Simulations. Molecular dynamics simulations were carried out with the Gromacs package and all restraints were enforced using the PLUMED plugin. 75

8


Page 8 of 34

Page 9 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For each of the liquid phase environments, we performed 35 different molecular dynamics simulations using the QAS and the QAC. Each individual simulation was performed for 150 ns for a total of 5.25 µs of simulation data for both the QAC and QAS systems. For a complete description of the simulation settings and equilibration procedure, please refer to the Supporting Information. Analysis. In order to determine which kind of cage type the AA binds to throughout the molecular dynamics trajectories, either 512 or 512 64 for a sII hydrate, we first needed to be able to distinguish between these cage types. To distinguish the cage types on the hydrate surface, we have utilized the order parameters of Jacobson et. al., 76 which are order parameters specifically developed for this purpose. These order parameters utilize a geometric based search algorithm that first finds the five- and six-membered rings of water molecules in a 3.5 Å radius, and then searches for five- and six-membered rings that are connected on each edge. Since 512 and 512 64 cage types each have a unique geometry, the algorithm is able to discern the two cage types as well as the water molecules that comprise the cages. Furthermore, it can accurately classify partially formed cages of each type. A complete description of the order parameters is given in Ref. [ 76 ]. For each molecular dynamics trajectory, we have used these order parameters to classify the cage types of the binding sites on the hydrate surface. In addition, once we have identified all the binding sites on the hydrate surface, we monitored the distance between the centroid of each cage type binding site on the hydrate surface and each of the four terminal carbon groups on the AA, i.e. the three terminal carbon groups on the short chain butyl groups as well as the terminal carbon group on the longer hydrocarbon chain. We classified the binding configurations based on what cage type each terminal carbon group was bound to at any given moment in time. The binding sites within the cages on the surface of the hydrate were defined as the region of space within 2.0 Å of each cage centroid. In addition, we defined each of the terminal carbon groups to be bound to the hydrate surface if they were within 2.0 Å of the cage centroid for longer than a 100 ps timeframe. We chose to define that a hydrocarbon chain was bound when it was within

9



2.0 Å of the centroid because we found that this distance was approximately the distance of fluctuations with respect to the centroid when the hydrocarbon chains were bound to various cage structures. This distance also led to the best classification of bound configurations when we tested known bound structures. Moreover, by requiring that a structure be bound for longer than 100 ps, our definition eliminates the possibility of classifying random fluctuations of the hydrocarbon chains as bound configurations. This procedure allowed us to obtain a set of binding configurations and their corresponding energies in each liquid phase environment. After analyzing all the molecular dynamics simulation data, we computed the binding energy, or surface adsorption energy, of each binding configuration. Typically, the binding energy is computed as the difference between the ensemble-averaged total system energy of the AA away from the hydrate surface (unbound state) and the ensemble-averaged total system energy of the AA bound to the hydrate surface (bound state). However, with this approach the unbound state of the AA, generally defined to be an arbitrary distance from the hydrate surface, can be difficult to define in a finite-sized system as the AA inevitably interacts with the hydrate surface through long-range electrostatics. This approach can yield binding energies which are dependent on the definition of the unbound state. Furthermore, it can be difficult to obtain adequate sampling of the AA at this arbitrarily specified distance if the molecule diffuses quickly through the liquid phase environment, which proved to be the situation in the liquid hydrocarbon environment. Consequently, we have computed the binding energy of the AA using a slightly different approach that allows us to define the unbound state of the AA to be within the bulk of the liquid phase, i.e. significantly far from the hydrate surface. This approach, which is described in detail in the Supporting Information, has the advantage that the binding energy of the AA is not arbitrarily dependent on the definition of the unbound state and it allows us to adequately sample the unbound state. Binding Free Energy Simulations. The binding free energy of the AA binding to a hydrate surface can be obtained by calculating the free energy difference associated with

10


Page 10 of 34

Page 11 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


decoupling the electrostatic and van der Waals (VDW) potentials of the AA from both the liquid phase (∆GLP ) and the hydrate surface (∆GHS ). These calculations can be performed with FEP and the binding free energy can be computed from the difference between these quantities, i.e. ∆Gb = ∆GHS − ∆GLP . However, care must be taken when performing these calculations. When decoupling a molecule from its environment, the molecule is free to drift away from its binding site, and consequently, for the reverse coupling transformation, the molecule is unlikely to form in the relevant binding site, thereby violating microscopic reversibility. This problem is commonly known as the wandering-ligand problem in the literature. 60–62 To circumvent this problem, a suitable set of geometric restraints (conformational, translational, and rotational) are introduced. These sets of restraints prevent the molecule from wandering during alchemical transformations and also help to aid in convergence of the calculations by restricting the conformational search space. 60 However, these restraints bias the free energy calculation, and so to properly account for this bias, a thermodynamic cycle is constructed incorporating these restraints so that their corresponding free energy contributions are taken into account. This approach is known as the double decoupling method and has been used in numerous studies to calculate the binding free energy of ligands to protein receptors 60,61 as well as to membrane surfaces. 62 To calculate the standard binding free energy, ∆Gob , of the AA to the hydrate surface in each liquid phase environments, we have utilized the double decoupling procedure outlined by Deng and coworkers. 60,61 The thermodynamic cycle for this procedure is shown in Fig. 2. From this thermodynamic cycle, the standard binding free energy can be calculated as

∆Gob = ∆Gint + ∆Gc + ∆Got + ∆Gr ,

(1)

HS LP LP LP o where ∆Gint = [(∆GHS − ∆GHS int + ∆Gcorr ) − (∆Gint + ∆Gcorr )], ∆Gc = [∆Gc c ], ∆Gt = o HS HS LP [−∆GHS t −kB T ln(Ft C )], and ∆Gr = [−∆Gr −kB T ln(Fr )]. The ∆Gcorr and ∆Gcorr terms

are analytic correction terms that are applied to correct for finite-size effects of the raw ionic 11



Figure 2: Thermodynamic cycle used to compute the binding free energy ∆Gb according to the alchemical double decoupling method with restraining potentials. This thermodynamic cycle is a linked sequence of thermodynamic processes that connects two end-point states: the unbound state S+AA with the hydrate surface (S) and the unrestricted AA in the liquid phase (AA), and the bound state SAA. In the figure, c, t, and r in parenthesis denote when conformational, translational, and rotational restraints are applied and vac denotes when the AA molecule is in a noninteracting decoupled (vacuum) state. See Ref. [ 60 ] and the Supporting Information for a complete description of each term. charging free energies of the thermodynamic cycle, removing their dependence on the size of the simulation cell. 63,64 These terms correct for finite-size effect errors involving interactions and undersolvation of the AA and hydrate surface in the reference computational box, its periodic replicas, and the homogeneous neutralizing background charge density filling the infinite periodic system. In addition, Ft and Fr are the factors associated with applying translational and rotational restraints on the AA when the intermolecular interactions are decoupled, and the constant C o ensures conversion to the standard state concentration (C o =1 M or 1/1661 Å−3 ). Given the extensive number of studies that have utilized this method, we refer the reader to the literature for a detailed mathematical description 60,61 and we present a more detailed description of our free energy calculations following this scheme in the Supporting Information.

12


Page 12 of 34

Page 13 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results and Discussion The binding configurations, binding energies, and residences times for the QAC system and QAS system in each of the different liquid phase environments are shown in Table 1. The configurations are ordered from strongest binding energy to weakest binding energy and they are labeled with a configuration key. The configuration key designates what cage type each of the hydrocarbon chains of the AA are bound to in the configuration. The first position of the configuration key (prior to the hyphen) designates the long hydrocarbon chain of the AA while the remaining three positions of the key (after the hyphen) designate the three short hydrocarbon chains of the AA. The key has a marker U if the hydrocarbon chain is unbound, it has a marker S if the hydrocarbon chain is bound to a smaller 512 cage, and it has a marker L if the hydrocarbon chain is bound to a larger 512 64 cage. For example, the key L−USL indicates that the long chain is bound to a 512 64 cage, one short chain is unbound, one short chain is bound to a 512 cage, and one short chain is bound to a 512 64 cage. Note that since the short chains are all of the same length there can be degenerate configurations, such as L−USL and L−SLU, and therefore the configuration key is not unique. Consequently, when calculating the binding energy of each configuration, we have combined all energies corresponding to degenerate configurations into the ensemble averages. The negative binding energies shown suggest that the adsorption of the AA to the hydrate surface is an exothermic process. In addition, in comparing the results for the QAC system and QAS system in Table 1, the binding energies and residence times are very similar for each binding configuration in either liquid phase, indicating that the presence of the chloride anion has no appreciable impact on the surface adsorption of the QAC, or at the very least, on the overall energetic order of the binding configurations. However, there is a noticeable trend that each binding configuration in Table 1 (c)-(d) has a lower energy when compared to its corresponding counterpart in Table 1 (a)-(b). This suggests that the use of an explicit chloride anion leads to a non-negligible energy stabilization of the binding configurations, or equivalently, that the use of a uniform background charge understabilizes the binding con13



(a)

S−USL L−USL S−UUS L−SLL U−USL S−UUL U−SLL

∆E (kcal/mol) -16.4 ±0.5 -13.3 ±0.5 -10.9 ±2.8 -9.5 ±1.9 -5.2 ±0.4 -3.5 ±2.5 -0.3 ±0.7

(b)

Time (ns) 765.8 440.0 26.4 23.6 862.5 23.7 179.4

U−UUL U−ULL U−UUS U−USL

∆E (kcal/mol) -27.4 ± 0.6 -26.9 ± 1.1 -21.9 ± 0.9 -4.7 ± 0.4

(c)

S−USL L−USL U−USL S−UUS S−UUL U−SLL

∆E (kcal/mol) -20.8 ±0.7 -13.7 ±0.5 -8.6 ±0.4 -7.3 ±0.8 -3.8 ±2.2 -1.3 ±1.2

Page 14 of 34

Time (ns) 1611.7 172.6 647.9 756.3

(d)

Time (ns) 438.3 608.6 899.2 183.8 14.8 47.3

U−UUL U−UUS U−USL

∆E (kcal/mol) -35.7 ±0.9 -27.6 ±0.8 -3.3 ±0.8

Time (ns) 1762.3 1472.9 1111.4

Table 1: Binding configurations, binding energies, and residence times for the QAC system in (a) aqueous solution and (b) liquid hydrocarbon solution, and for the QAS system in (c) aqueous solution and (d) liquid hydrocarbon solution. The binding configurations are denoted by a configuration key which designates the polyhedral cage type each of the hydrocarbon chains of the AA are bound to in the configuration (U for unbound, S for 512 , and L for 512 64 ). The first position of the key represents the long hydrocarbon chain, while the remaining three positions of the key represent the short hydrocarbon chains. figurations. This makes intuitive sense since the neutralizing background charge introduces no electrostatic gradient and has a uniform distribution whereas the chloride anion exerts an electrostatic gradient and can directly interact with the QAC. Although the use of an explicit chloride anion makes little difference on the surface adsorption, it is evident that there are clear differences between the energies and binding configurations in either of the liquid phases. As is shown in Table 1 (a) and Table 1 (c), the most energetically favorable binding configurations in the aqueous phase are those with the long chain and two short chains of the AA bound to the hydrate surface (see Fig. 3 (a)-(b)). These results are consistent with the binding configurations proposed in the mechanism of Sun and Firoozabadi for AA inhibition at 100% watercut, whereby the long chain of the AA

14


Page 15 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecule binds to the hydrate surface (thereby covering the surface) and prevents agglomeration. 23 By taking the difference in the binding energies of the configurations L−USL and U−USL in Table 1 (a), we see that configurations with the long chain of the AA bound to the hydrate surface are roughly 5 kcal/mol lower in energy than configurations with the long chain unbound. This energy stabilization occurs because the binding of the long chain to the hydrate surface minimizes the disruption of the hydrogen bond network and allows for more hydrogen bonds to form in the aqueous phase, thereby lowering the energy of the system. Despite the fact that binding of the long chain is more energetically favorable, we observed that the U−USL binding configuration (see Fig. 3 (c)) has the longest residence time. In all simulations, we observed that the short chains of the AA bind first due to their proximity to the ammonium cation, which drives the surface adsorption of the AA through electrostatic interactions. This is supported by our results in Table 1 where it is apparent that configurations where only the long chain is bound (such as S-UUU or L-UUU) are noticeably absent. In addition, at no point during any of our simulations did we observe a binding configuration where the hydrocarbon chains of the AA bound to a cage that was occupied by a methane or propane molecule. Since the hydrocarbon chains generally penetrate deep into the cage structures, the steric clashes with guest molecules prevent it from binding to cages that are occupied. Consequently, after the initial binding event of the short hydrocarbon chains and given the finite length of the long hydrocarbon chain, it can only bind to the surface under two circumstances: if there is already an unoccupied cage within the vicinity of the initial binding site or if at some point after the primary binding event of the short chains, a gas molecule diffuses out of a nearby cage, thus allowing the long chain to bind. Therefore, the surface adsorption of the AA to the lowest energy binding configuration in the aqueous phase occurs through a two-step mechanism. Our results in Table 1 (a) and Table 1 (c) demonstrate that the U−USL binding configuration is the most likely initial binding configuration in the two-step mechanism, followed by secondary binding of the long chain to a 512 cage, which yields the lowest energy S−USL binding configuration. Given the mechanism

15



Page 16 of 34

(a)

(b)

(c)

(d)

Figure 3: The lowest energy binding configuration in the aqueous phase system, S-USL, is shown in (a) (side-view) and (b) (top-view), whereas the binding configuration, U-USL, with the longest residence time in the aqueous phase is shown in (c). The lowest energy binding configuration in the liquid hydrocarbon system, U-UUL, is shown in (d). The thin liquid water film that formed in the liquid hydrocarbon simulations is clearly present in (d). In addition, the binding configurations in (c) and (d) were used as the initial configurations in the binding free energy calculations in the aqueous and liquid hydrocarbon phases, respectively. Note that the liquid water and liquid hydrocarbon phases are omitted for clarity. proposed by Sun and Firoozabadi for AA inhibition in a 100% watercut environment, the two-step nature of the surface adsorption process implies that the affinity of the long chain for the hydrate surface may significantly influence the AA effectiveness. In Table 1 (b) and Table 1 (d) the binding configuration results from our molecular dynamics simulations in the liquid hydrocarbon phase are presented. In contrast to the binding configurations in the aqueous phase, in the liquid hydrocarbon phase the long chain of the AA remained unbound in all binding configurations. This suggests that surface adsorption in the liquid hydrocarbon system proceeds through a single-step mechanism, since there are no secondary binding events of the long hydrocarbon chain. Furthermore, in all binding configu-

16


Page 17 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


rations, the long chain remains solvated in the liquid hydrocarbon phase, which is consistent with the commonly proposed AA mechanism in the literature whereby the QA molecules adsorb onto the hydrate surface forming a reverse micelle. 1,2 We found that a thin liquid water film (∼ 5-7 Å thick) formed on the surface of the hydrate during equilibration in the liquid hydrocarbon system in all simulations and is visible in Fig. 3 (d). Thin liquid films have been reported in previous studies 77 of hydrates involving two-phase systems, particularly in hydrate/gas or hydrate/liquid-hydrocarbon systems. In addition, there is evidence throughout the literature that this quasi-liquid layer at the hydrate interface can drive hydrate particle aggregation through capillary liquid bridges that form between adjacent hydrate particles. 25,78–80 Interestingly, we observed that this liquid water layer influenced the binding configurations in the liquid hydrocarbon system. The preferred binding configuration, U-UUL, in this system only had one short chain bound to a larger 512 64 cage on the hydrate surface. This was the lowest energy configuration and had the longest residence time of all binding configurations. The U-UUL configuration was more energetically favorable because it allowed for the ammonium cation to be more exposed to negatively charged oxygen atoms in the liquid water layer. The water molecules in the liquid layer can solvate the ammonium cation more effectively than the water molecules in the hydrate surface (thereby shielding the charge from the liquid hydrocarbon phase) since they are mobile and not restrained to the hydrate cage geometries. These additional electrostatic interactions with the ammonium cation outweigh any energetic benefit of binding more short chains to the hydrate surface. In addition, in the QAS system, we found that this interfacial water layer was important for facilitating the binding of the QAS to the hydrate surface. In our simulations, the QAC and chloride anion remained associated as a salt when dissolved in the liquid hydrocarbon phase and diffused toward the hydrate surface. Upon interaction with the interfacial water layer, the QAS dissociated into a QAC and chloride anion, and subsequently, the QAC bound to the hydrate surface while the chloride anion remained solvated in the interfacial water layer. Although there are differences in the structures of the binding configurations in either

17



Page 18 of 34

(a)

∆GHS int AQ -30.9 HC -44.0

∆GLP int -28.3 -28.0

∆GHS int

∆GLP int

-103.2 -92.2

-101.7 -15.1

∆Gel -1.0 -17.1

∆GV DW -1.6 1.2

∆Gc 1.2 0.3

∆Got 0.1 0.9

∆Gr 2.7 2.1

∆GHS corr 0.1 4.3

∆GLP corr 0.2 5.3

∆Gob 1.4 ±1.4 -13.7 ±1.1

∆Got 0.5 0.6

∆Gr 2.9 2.6

∆GHS corr 0.1 4.3

∆GLP corr 0.2 5.3

∆Gob 3.2 ±1.3 -74.5 ±2.0

(b)

AQ HC

∆Gel 0.8 -78.5

∆GV DW -2.3 1.4

∆Gc 1.3 0.4

Table 2: Binding Free Energies for (a) QAC system and (b) QAS system in the aqueous and HS LP liquid hydrocarbon phases. Note that ∆Gel = ∆GHS el − ∆Gel and ∆GV DW = ∆GV DW − ∆GLP V DW . These quantities represent the free energy change of the electrostatic and VDW intermolecular interactions upon surface adsorption of the AA to the hydrate from the liquid phase. All energies are reported in kcal/mol. liquid phase, the main difference between these systems is in the energetics of the binding configurations. In particular, the lowest energy binding configuration in the liquid hydrocarbon system is roughly 11-15 kcal/mol lower in energy than the lowest energy binding configuration in the aqueous phase. The stronger binding affinity implies increased surface adsorption of the AA in the liquid hydrocarbon phase relative to the aqueous phase. However, to properly characterize the extent of the surface adsorption in each liquid phase environment, we computed the binding free energies for representative equilibrium binding configurations. Since the equilibrium constant for the adsorption process is a function of the binding free energy, the binding free energy should give us an understanding of the surface coverage of the AA on the hydrate. For our binding free energy calculations, we chose the U−USL and U−UUL binding configurations as the initial configurations in the aqueous and liquid hydrocarbon phases, respectively (see Fig. 3 (c) and (d)). Although the U−USL binding configuration is not the lowest energy binding configuration in the aqueous phase, we chose it as the initial configuration for the binding free energy calculation because it has the largest residence time and is therefore the most likely initial binding configuration in the two-step binding process in the aqueous phase. Consequently, the binding free energy associated with this configuration is more indicative of the surface adsorption of the AA. 18


Page 19 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The results of the binding free energy studies of the QAC system and QAS system in both liquid phases are shown in Table 2. In comparing the results for the QAC system and QAS system in Table 2 for the aqueous phase (∆Gob = 1.4 kcal/mol and ∆Gob = 3.2 kcal/mol) we see that our results are very similar. For the aqueous phase, we found that the QAS remained dissociated as a QAC and a chloride anion when the AA was bound to the hydrate surface and when dissolved in the bulk aqueous phase. Moreover, since the chloride anion interacts with liquid water in both the bound and unbound states, its contribution to the free energy largely cancels out between both alchemical legs of the thermodynamic cycle. As a result, the binding free energy for a QAC and a QAS absorbing to the hydrate surface are very similar, and taking into consideration the error in the calculation, these results do not differ significantly. However, in comparing our results in Table 2 for the QAC and QAS systems in the liquid hydrocarbon phase (∆Gob = -13.7 kcal/mol and ∆Gob = -74.5 kcal/mol) the results do differ significantly. In the bulk liquid hydrocarbon phase, we found that the QAS remained associated as a salt whereas when the QAS absorbed to the interface, it first dissociated into a QAC and chloride anion in the interfacial water layer, and subsequently the QAC bound to the hydrate surface while the chloride anion remained solvated within the water layer. Since the chloride anion primarily interacts with a liquid hydrocarbon phase in the bulk and an aqueous phase upon absorption, there is a large contribution to the free energy that is approximately equal to the difference of the solvation free energies of the chloride anion in these two liquid environments. This contribution is only approximately equal to the difference of solvation free energies because in the bulk liquid hydrocarbon phase the AA remains associated in salt form, and so the chloride anion strongly interacts with the QAC in addition to the liquid hydrocarbon phase. Given that the chloride anion is much more soluble in the aqueous phase, there is a large energy stabilization upon adsorption of the QAS to the hydrate interface, which accounts for the significantly different binding free energies of the QAC and QAS systems. Regardless if we consider the AA to be a QAC or a QAS, we find that surface adsorption

19



of the AA in the liquid hydrocarbon phase is an exergonic process and highly favorable, whereas surface adsorption in the aqueous phase is an endergonic process and unfavorable. Furthermore, since the binding free energy is significantly lower in energy in the liquid hydrocarbon phase, the extent of surface adsorption is much larger, leading to a greater degree of surface coverage of the hydrate and likely to improved AA performance. To gain insight into these interaction energies and the underlying driving forces for surface adsorption, we have calculated the free energy change of the electrostatic and VDW interaction terms LP upon surface adsorption of the AA from the liquid phases, i.e. ∆Gel =∆GHS el − ∆Gel and LP ∆GV DW = ∆GHS V DW − ∆GV DW . As shown in Table 2, the free energy change in the elec-

trostatic interactions upon surface adsorption is significantly lower in energy in the liquid hydrocarbon phase, ∆Gel = -17.1 kcal/mol and ∆Gel = -78.5 kcal/mol, than in the aqueous phase, ∆Gel = -1.0 kcal/mol and ∆Gel = 0.8 kcal/mol. This stronger electrostatic driving force is not surprising given the non-polar nature of the liquid hydrocarbon phase and the polar nature of the hydrate surface. Interestingly, our results show that there is a larger VDW driving force for binding in the aqueous phase. However, the unfavorable interactions of the long chain of the AA with the aqueous phase imposes an energetic penalty, which is apparent from the large free energy values in Table 2 associated with imposing the restraints. These energetic penalties outweigh the favorable intermolecular interactions of the AA with the hydrate surface, leading to an overall positive binding free energy. In addition to analyzing the energetics of the surface adsorption, we have examined the structural effects of AA binding by calculating the radial distribution function, g(r), from our molecular dynamics simulation data in both liquid phases. Figure 4 (a) shows the radial distribution function between the nitrogen atom on the inhibitor molecule and the oxygen atoms of water in the aqueous phase and in the hydrate phase. It is clear from Fig. 4 (a) that the ammonium cation of the AA interacts very differently with water when it is bound to the hydrate surface as opposed to being solvated in the aqueous phase. The larger and broader peak around 4.35 Å suggests that there are more water interactions with the ammonium

20


Page 20 of 34

Page 21 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

Figure 4: Radial distribution functions calculated from molecular dynamics simulation data. (a) Comparison of the radial distribution functions of the nitrogen on the AA and the oxygen atoms of water in the aqueous and hydrate phases. (b) Comparison of the radial distribution functions of the nitrogen on the AA and oxygen atoms of water when the AA is bound to the hydrate surface as well as the oxygen-oxygen radial distribution function of the water molecules in the hydrate. (c) Comparison of the radial distribution functions of the nitrogen on the AA and the oxygen atoms of water when the AA is bound to the hydrate surface in both the aqueous and liquid hydrocarbon phases. cation when it is solvated in the aqueous phase. However, the smaller but more sharply defined peak, which is shifted to the left at 4 Å, indicates slightly stronger interactions when the inhibitor is bound to the hydrate surface. In Fig. 4 (b) we compare the nitrogen-oxygen (N -O) radial distribution function when the inhibitor is bound to the hydrate surface and the oxygen-oxygen radial distribution function of the hydrate. These two radial distribution functions are quite different, which demonstrates that the ammonium cation does not fit into the water structure of the hydrate lattice when it is bound. These results are in contrast

21



to previous studies with various kinetic inhibitors, such as PVCap, which do insert into the water structure of the hydrate upon binding to the hydrate surface. 49 Since each of the N -O radial distribution functions have peaks that begin at around 4 Å (See Fig. 4(c)), it is likely that the ammonium cation of the AA does not insert into the hydrate lattice due to steric hindrance stemming from the short hydrocarbon chains. Finally, in Fig. 4 (c), we compare the N -O radial distribution functions when the inhibitor is bound to the hydrate surface in each liquid phase. There is a correlation between these radial distribution functions but it is apparent that the N -O radial distribution function associated with the liquid hydrocarbon phase has slightly larger and broader peaks by comparison. Moreover, the N -O radial distribution function of the inhibitor bound to the hydrate surface in the liquid hydrocarbon phase also correlates well with the N -O radial distribution function of the unbound inhibitor in the aqueous phase (see Fig. 4 (a)). This suggests that the broadening of the radial distribution function in the liquid hydrocarbon phase is a result of the increased interactions of the ammonium cation with the thin liquid water film that is present on the hydrate surface. Furthermore, the larger peak suggests there are more interactions with the water molecules in the thin film on the hydrate surface, and provides evidence that the thin liquid film influences the binding configurations in the liquid hydrocarbon phase.

Conclusions We have examined the surface adsorption mechanism of a model anti-agglomerant inhibitor to a hydrate surface in both aqueous and liquid hydrocarbon phases using molecular dynamics simulations. The model anti-agglomerant used in our simulations was n-dodecyl-tri(nbutyl)-ammonium chloride, a quaternary ammonium (QA) molecule that has been shown to be an effective anti-agglomerant inhibitor. 1,3,18 Moreover, we performed simulations in each liquid phase environment utilizing the QA salt (QAS) and the QA cation (QAC), i.e. with and without an explicit chloride anion, respectively, and showed that the inclusion of

22


Page 22 of 34

Page 23 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


an explicit chloride anion did not make a significant impact on our results. Using molecular dynamics simulations, we were able to gain a molecular-level understanding of how the different liquid phase environments influenced the surface adsorption of the inhibitor, which is likely a crucial component of AA effectiveness. From our simulation data, we classified all the binding configurations of the AA on the hydrate surface and calculated their corresponding binding energies. In the aqueous phase, we demonstrated that the lowest energy surface adsorption pathway of the inhibitor proceeds through a two-step mechanism: the AA initially binds to the hydrate surface with two of the short hydrocarbon chains bound to the hydrate surface, followed by a secondary binding of the long hydrocarbon chain, which is contingent upon there being an empty cage structure within the vicinity of the original binding site. Furthermore, we found that binding configurations with the long chain bound to the surface are roughly 5 kcal/mol lower in energy than binding configurations with the long chain unbound. The two-step nature of the surface adsorption process in the aqueous phase implies that the affinity of the long chain for the hydrate surface may significantly influence the AA effectiveness. In the liquid hydrocarbon phase, we demonstrated that surface adsorption occurs through a single-step mechanism. In all binding configurations, we found that the long chain remains unbound and solvated in the liquid hydrocarbon phase while the ammonium cation head of the molecule bound to the surface and was solvated by a quasi-liquid layer that formed naturally as a result of equilibration of the system. These results support the commonly proposed AA mechanism in the literature whereby the QA molecules adsorb onto the hydrate surface to stabilize slurries of the hydrates in the produced fluids. 1,2 In addition to stabilizing the water-in-oil emulsions, the AA may destabilize capillary liquid bridges that form between hydrate particles by disrupting hydrogen bond networks that form within these liquid bridges, thus preventing agglomeration. However, further studies will be required to evaluate this mechanism. To characterize the extent of surface adsorption in each liquid phase, we calculated the standard binding free energy using the free energy perturbation (FEP) method

23



following a double decoupling thermodynamic cycle. We found that the surface adsorption in the liquid hydrocarbon phase is an exergonic process whereas the surface adsorption in the aqueous phase is an endergonic process. Our results demonstrate that the extent of surface adsorption is much larger in the liquid hydrocarbon phase relative to the aqueous phase, and suggest that the AA is less effective in the aqueous phase because the surface adsorption is very low. The high affinity of the AA for absorption to the hydrate interface in the liquid hydrocarbon phase suggests that a large concentration of AAs would necessarily bind to and occupy the hydrate interface. In the aqueous phase, the low affinity of the AA for adsorption to the hydrate surface is likely due to the similar electrostatic nature of the aqueous phase and the hydrate surface, and suggests there is a minimal concentration of AAs that bind to and occupy the hydrate interface. Consequently, there is little to prevent hydrate particles from agglomerating in the aqueous phase. Finally, we examined the structural effects of the inhibitor on the water structure in the liquid phase and in the hydrate phase using radial distribution functions. In contrast to many kinetic inhibitors, 49 we have shown that the ammonium cation does not incorporate into the local water structure of the hydrate lattice upon binding due to steric hindrance stemming from the short hydrocarbon chains of the AA. To the extent that surface adsorption corresponds to AA performance, our results are consistent with the trend in experiments in which AA’s perform better in oil-dominated rather than water-dominated systems. Future work should include benchmarking specific AA’s or a single AA in specific environments against experimental conditions of known performance to search for microscopic quantities that are predictive of AA performance, perhaps streamlining the innovation process. Additionally, the effect of surface coverage, as investigated in Ref. [ 58,59 ] on binding mechanism and anti-agglomeration should be studied across a range of conditions (fresh water, brine, oil, AA concentration, AA type, hydrate crystal type and crystal face), ideally including energetic analysis and mechanistic insight into the AA-hydrate interaction and the nature of anti-agglomeration itself.

24


Page 24 of 34

Page 25 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Associated Content Supporting Information. The Supporting Information is available free of charge on the ACS Publications website (http://pubs.acs.org). Additional information related to force fields; Molecular dynamics simulation settings; Binding energy and statistical error calculations; Binding free energy calculations.

Acknowledgements The authors thank Andrew Huy Nguyen and Dr. Valeria Molinero for the order parameter code that classifies polyhedral cage sites of the hydrate. The authors would also like to thank Drs. Gregg Beckham, Tony Spratt, Peter Conrad, Eric Grzelak, Hari Subramani, Siva Subramanian, Douglas Estanga, Nicole Bernstein, and Sadegh Aryanpanah for useful discussions about this work, and Drs. Amadeu Sum, David Wu, James Ely, and Scott Wierzchowski for early suggestions on model selection. The authors would kindly acknowledge financial support from Chevron-MIT Energy Initiative program Alliance.

References (1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases; CRC Press: Boca Raton, FL, 2008. (2) Koh, C. A. Towards a Fundamental Understanding of Natural Gas Hydrates. Chem. Soc. Rev. 2002, 31 , 157-167. (3) Kelland, M. A. History of the Development of Low Dosage Hydrate Inhibitors. Energy Fuels 2006, 20, 825. (4) Buffett, B. A. Clathrate Hydrates. Annu. Rev. Earth Planet. Sci. 2000, 28, 477-507.

25



(5) Kvenvolden, K. A. Gas Hydrates-Geological Perspective and Global Change. ReV. Geophys. 1993, 31, 173-187. (6) Park, Y.; Kim, D. Y.; Lee, J. W.; Huh, D.G.; Park, K. P.; Lee, J.; Lee, H. Sequestering Carbon Dioxide into Complex Structures of Naturally Occurring Gas Hydrates. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 12690-12694. (7) Glasby, G. P. Potential Impact on Climate of the Exploitation of Methane Hydrate Deposits Offshore. Mar. Pet. Geol. 2003, 20, 163-175. (8) Reagan, M. T.; Moridis, G. J. Global Climate and the Response of Oceanic Hydrate Accumulations. Fire in the Ice: Methane Hydrate Newsletter 2010, 10, 9-12. (9) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Stable Low-Pressure Hydrogen Clusters Stored in a Binary Clathrate Hydrate. Science 2004, 306, 469-471. (10) Mao, W. L.; Mao, H. K.; Goncharov, A. F.; Struzhkin, V. V.; Guo, Q.; Hu, J.; Shu, J.; Hemley, R. J.; Somayazulu, M.; Zhao, Y. Hydrogen Clusters in Clathrate Hydrate. Science 2002, 297, 2247-2249. (11) Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Fundamentals and Applications of Gas Hydrates. Annu. Rev. Chem. Bio. Eng. 2011, 27, 237-257. (12) Chatti, I.; Delahaye, A.; Fournaison, L.; Petitet, J. P. Benefits and Drawbacks of Clathrate Hydrates: A Review of Their Areas of Interest. Energy Convers. Manage. 2005, 46, 1333-1343. (13) Sloan, E. D. Fundamental Principles and Applications of Natural Gas Hydrates. Nature (London, U.K.) 2003, 426, 353-363. (14) Collett, T. S.; Kuuskraa, V. A. Hydrates Contain Vast Store of World Gas Resources. Oil Gas J. 1998, 96, 90-95. 26


Page 26 of 34

Page 27 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Lee, S. Y.; Holder, G. D. Methane Hydrates Potential as a Future Energy Source. Fuel Process. Technol. 2001, 71, 181-186. (16) Bagirov, E.; Lerche, I. Hydrates Represent Gas Source, Drilling Hazard. Oil Gas J. 1997, 95, 99-104. (17) Englezos, P. Clathrate Hydrates. Ind. Eng. Chem. Res. 1993, 32, 1251-1274. (18) Chua, P. C.; Kelland, M. A. Study of the Gas Hydrate Anti-Agglomerant Performance of a Series of n-Alkyl-tri(n-butyl)ammonium Bromides. Energy Fuels 2013, 27, 12851292. (19) Graham, B.; Reilli, W. K.; Beinecke, F.; Boesch, D. F.; Garcia, T. D.; Murray, C. A.; Ulmer, F. Deep Water: The Gulf Oil Disaster and the Future of Offshore Drilling. Report to the President; 2011. (20) Koh, C. A.; Westacott, R. E.; Zhang, W.; Hirachand, K.; Creek, J. L.; Soper, A. K. Mechanisms of Gas Hydrate Formation and Inhibition. Fluid Phase Equilib. 2002, 194, 143-151. (21) Creek, J. L. Efficient Hydrate Plug Prevention. Energy Fuels 2012, 26, 4112-4116 (22) Gao, S. Q. Hydrate Risk Management at High Watercuts with Anti-Agglomerant Hydrate Inhibitors. Energy Fuels 2011, 11, 3343-3351. (23) Sun, M; Firoozabadi, A. New Surfactant for Hydrate Anti-Agglomeration in Hydrocarbon Flowlines and Seabed Oil Capture. J. Colloid Interface Sci. 2013, 402, 312-319. (24) Somasundaran, P; Krishnakumar, S. Adsorption of Surfactants and Polymers at the Solid-Liquid Interface. Colloids Surf., A. 1997, 123, 491-513. (25) Aman, Z. M.; Koh, C. A. Interfacial Phenomena in Gas Hydrate Systems. Chem. Soc. Rev. 2016, 45, 1678-1690. 27



(26) Jensen, L.; Thomsen, K.; von Solms, N.; Wierzchowski, S.; Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. Calculation of Liquid Water-Hydrate-Methane Vapor Phase Equilibria from Molecular Simulations. J. Phys. Chem. B 2010, 114, 5775-5782. (27) Conde, M. M.; Vega, C. Determining the Three-Phase Coexistence Line in Methane Hydrates using Computer Simulations. J. Chem. Phys. 2010, 133, 064507. (28) Conde, M. M.; Vega, C.; McBride, C.; Noya, E. G.; Ramirez, R.; Sese, L. M. Can Gas Hydrate Structures be described using Classical Simulations? J. Chem. Phys. 2010, 132, 114503. (29) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Gas Hydrate Nucleation and Cage Formation at a Water/Methane Interface. Phys. Chem. Chem. Phys. 2008, 10, 4853. (30) Vatamanu, J.; Kusalik, P. G. Observation of Two-Step Nucleation in Methane Hydrates. Phys. Chem. Chem. Phys. 2010, 12,15065. (31) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science 2009, 326, 1095. (32) Walsh, M. R.; Beckham, G. T.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. Methane Hydrate Nucleation Rates from Molecular Dynamics Simulations: Effects of Aqueous Methane Concentration, Interfacial Curvature, and System Size. J. Phys. Chem. C 2011, 115, 21241. (33) Jimenez-Angeles, F.; Firoozabadi, A.; Nucleation of Methane Hydrates at Moderate Subcooling by Molecular Dynamics Simulations. J. Phys. Chem. C 2014, 118, 1131011318 (34) Barnes, B. C.; Knott, B. C.; Beckham, G. T.; Wu, D. T.; Sum, A. K. Reaction Coor-

28


Page 28 of 34

Page 29 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dinate of Incipient Methane Clathrate Hydrate Nucleation. J. Phys. Chem. B 2014, 118, 13236-13243. (35) Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. J. Am. Chem. Soc. 2010, 132, 11806. (36) Jacobson, L. C.; Molinero, V. Can Amorphous Nuclei Grow Crystalline Clathrates? The Size and Crystallinity of Critical Clathrate Nuclei. J. Am. Chem. Soc. 2011, 133, 6458. (37) Debenedetti, P. G.; Sarupria, S. Chemistry. Hydrate Molecular Ballet. Science 2009, 326, 1070. (38) Radhakrishnan, R.; Trout, B. L. A New Approach for Studying Nucleation Phenomena using Molecular Simulations: Application to CO2 Hydrate Clathrates. J. Chem. Phys. 2002, 117, 1786. (39) Moon, C.; Taylor, P. C.; Rodger, P. M. Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 2003, 125, 4706. (40) Guo, G. J.; Li, M.; Zhang, Y. G.; Wu, C. H. Why Can Water Cages Adsorb Aqueous Methane? A Potential of Mean Force Calculation on Hydrate Nucleation Mechanisms. Phys. Chem. Chem. Phys. 2009, 11, 10427. (41) Christiansen, R. L.; Sloan, E. D. Mechanisms and Kinetics of Hydrate Formation. Ann. N.Y. Acad. Sci. 1994, 715, 283-305 (42) Knott, B. C.; Molinero, V.; Doherty, M. F.; Peters, B. Homogeneous Nucleation of Methane Hydrates: Unrealistic Under Realistic Conditions. J. Am. Chem. Soc. 2012, 134, 19544-19547. (43) Vatamanu, J.; Kusalik, P. G. Unusual Crystalline and Polycrystalline Structures in Methane Hydrates. J. Am. Chem. Soc. 2006, 128, 15588. 29



(44) Liang, S.; Kusalik, P. G. Crystal Growth Simulations of H2 S Hydrate. J. Phys. Chem. B 2010, 114, 9563. (45) Chihaia, V.; Adams, S.; Kuhs, W. F. Molecular Dynamics Simulations of Properties of a (001) Methane Clathrate Hydrate Surface. Chem. Phys. 2005, 317, 208. (46) Peters, B.; Zimmermann, N. E. R.; Beckham, G. T.; Tester, J. W.; Trout, B. L. Path Sampling Calculation of Methane Diffusivity in Natural Gas Hydrates from a WaterVacancy Assisted Mechanism. J. Am. Chem. Soc. 2008, 130, 17342. (47) English, N. J.; Tse, J. S. Mechanisms for Thermal Conduction in Methane Hydrate. Phys. Rev. Lett. 2009, 103, 015901. (48) English, N. J.; MacElroy, J. M. D. Perspectives on Molecular Simulation of Clathrate Hydrates: Progress, Prospects and Challenges. Chem. Eng. Sci. 2015, 121, 133-156. (49) Anderson, B. J.; Bazant, M. Z.; Tester, J. W.; Trout, B. L. Properties of Inhibitors of Methane Hydrate Formation via Molecular Dynamics Simulations. J. Phys. Chem. B 2005, 109, 8153. (50) Gomez Gualdron, D. A.; Balbuena, P. B. Classical Molecular Dynamics of ClathrateMethane-Water-Kinetic Inhibitor Composite Systems. J. Phys. Chem. C. 2007, 111, 15554-15564. (51) Carver, T. J.; Drew, M. G. B.; Rodger, P. M. Characterisation of the (111) Growth Planes of a Type II Gas Hydrate and Study of the Mechanism of Kinetic Inhibition by Poly(vinylpyrrolidone). J. Chem. Soc., Faraday Trans. 1996, 92, 5029-5033. (52) Storr, M. T.; Taylor, P. C.; Monfort, J. P.; Rodger, P. M. Kinetic Inhibitor of Hydrate Crystallization. J. Am. Chem. Soc. 2004, 126, 1569-1576. (53) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Adsorption Mechanism of Inhibitor and

30


Page 30 of 34

Page 31 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Guest Molecules on the Surface of Gas Hydrates. J. Am. Chem. Soc. 2015, 137, 1207912085. (54) Oluwunmi, P. A.; Finney, A. R.; Rodget, P. M. Molecular Dynamics Screening for New Kinetic Inhibitors of Methane Hydrate. Can. J. Chem. 2015, 93, 1043-1049. (55) Bagherzadeh, S. A.; Alavi, S.; Ripmeester, J. A.; Englezos, P. Why Ice-Binding Type I Antifreeze Protein Acts as a Gas Hydrate Crystal Inhibitor. Phys. Chem. Chem. Phys. 2015, 17, 9984-9990. (56) Bagherzadeh, S. A. Molecular Mechanisms of Methane Hydrate Dissociation and Inhibition, PhD Thesis, The University of British Columbia, March 2015. (57) Makogon, T. An Experimental and Computer Study of Gas Hydrate Formation and Inhibition, PhD Thesis, The Colorado School of Mines, January 1994. (58) Phan, A; Bui, T.; Acosta, E.; Krishnamurthy, P.; Striolo, A. Molecular Mechanisms Responsible for Hydrate Anti-Agglomerant Performance. Phys. Chem. Chem. Phys. 2016, 18, 24859-24871. (59) Bui, T.; Phan, A.; Monteiro, D.; Lan, Q.; Ceglio, M.; Acosta, E.; Krishnamurthy, P.; Striolo, A. Evidence of Structure-Performance Relation for Surfactants Used as Antiagglomerants for Hydrate Management. Langmuir 2017, 33, 2263-2274. (60) Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 2234-2246. (61) Wang, J.; Deng, Y.; Roux, B. Absolute Binding Free Energy Calculations using Molecular Dynamics Simulations with Restraining Potentials. Biophys. J. 2006, 91, 2798-2814. (62) Zhang, L.; Yethiraj, A.; Cui, Q. Free Energy Calculations for the Peripheral Binding of Proteins/Peptides to an Anionic Membrane. 1. Implicit Membrane Models. J. Chem. Theory Comput. 2014, 10, 2845-2859. 31



(63) Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hunenberger, P. A. Calculating the Binding Free Energies of Charged Species based on Explicit-Solvent Simulations Employing Lattice-Sum Methods: An Accurate Correction Scheme for Electrostatic Finite-Size Effects. J. Chem. Phys. 2013, 139, 184103. (64) Hummer, G.; Pratt, L. R.; Garcia, A. E. Free Energy of Ionic Hydration. J. Phys. Chem. 1996, 100, 1206-1215. (65) Lin, Y.; Baumketner, A.; Song, W.; Deng, S.; Jacobs, D.; Cai, W.; Ionic Solvation Studied by Image-Charge Reaction Field Method. J. Phys. Chem. 2011, 134, 044105. (66) Pfeiffer, S.; Fushman, D.; Cowburn, D. Impact of Cl- and Na+ Ions on Simulated Structure and Dynamics of βARK1 PH Domain. Proteins: Struct., Funct., Bioinf. 1999, 35, 206-217. (67) Ibragimova, G.; Wade, R.; Importance of Explicit Salt Ions for Protein Stability in Molecular Dynamics Simulation. Biophys. J. 1998, 74, 2906-2911. (68) Ibragimova, G.; Wade, R.; Stability of the β-Sheet of the WW Domain: A Molecular Dynamics Simulation Study. Biophys. J. 1999, 77, 2191-2198. (69) Hub J. S.; de Groot, B. L.; Grubmuller, H.; Groenhof, G.; Quantifying Artifacts in Ewald Simulations of Inhomogeneous Systems with a Net Charge. J. Chem. Theory Compute. 2014, 10, 381-390. (70) Soave, G. Equilibrium Constants From a Modified Redlich-Kwong Equation of State. Chem. Eng. Sci. 1972, 27 1197-1203. (71) Peneloux A.; Rauzy E.; Freze, R. A Consistent Correction for Redlich-Kwong-Soave Volumes. Fluid Phase Equilib. 1982, 8, 7-23. (72) Sarode, H. N.; Lindberg, G. E.; Yang, Y.; Felberg, L. E., Voth, G. A.; Herring, A. M.

32


Page 32 of 34

Page 33 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Insights into the Transport of Aqueous Quaternary Ammonium Cations: a Combined Experimental and Computational Study. J. Phys. Chem. B 2014, 118,1363-1372. (73) Babiaczyk, W. I.; Bonella, S.; Guidoni, L.; Ciccotti, G. Hydration Structure of the Quaternary Ammonium Cations. J. Phys. Chem. B 2010, 114, 15018-15028. (74) Klauda, J. B.; Sandler, S. I. Ab Initio Intermolecular Potentials for Gas Hydrates and their Predictions. J. Phys. Chem. B 2002, 106, 5722-5732 (75) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961-1972. (76) Jacobson, L. C.; Matsumoto, M.; Molinero, V. Order Parameters for the Multistep Crystallization of Clathrate Hydrates. J. Chem. Phys. 2011, 135, 074501. (77) Jimenez-Angeles, F.; Firoozabadi, A. Induced Charge Density and Thin Liquid Film at Hydrate/Methane Gas Interfaces. J. Phys. Chem. C 2014, 118, 26041-26048. (78) Dieker, L. E.; Aman, Z. M.; George, N. C.; Sum, A. K.; Sloan, E. D.; Koh, C. A. Micromechanical Adhesion Force Measurements between Hydrate Particles in Hydrocarbon Oils and Their Modifications. Energy Fuels 2009, 23, 5966-5971. (79) Aman, Z. M.; Dieker, L. E.; Aspenes, G.; Sum, A. K.; Sloan, E. D.; Koh, C. A. Influence of Model Oil with Surfactants and Amphiphilic Polymers on Cyclopentane Hydrate Adhesion Forces. Energy Fuels 2010, 24, 5441-5445. (80) Aman, Z. M.; Joshi, S. E..; Sloan, E. D.; Sum, A. K.; Koh, C. A. Micromechanical Cohesion Force Measurements to Determine Cyclopentane Hydrate Interfacial Properties. J. Colloid Interface Sci. 2012, 376, 283-288.

33



For Table of Contents Use Only Molecular Dynamics Analysis of Anti-Agglomerant Surface Adsorption in Natural Gas Hydrates Authors: Michael A. Bellucci, Matthew R. Walsh, and Bernhardt L. Trout

Figure 5: Table of Contents Figure.

34


Page 34 of 34

Molecular Dynamics Analysis of Anti-Agglomerant ... - ACS Publications

Recommend Documents