Towards the design of new and improved drilling fluid additives using molecular dynamics simulations

During exploration for oil and gas, a technical drilling fluid is used to lubricate the drill bit, maintain hydrostatic pressure, transmit sensor readings, remove rock cuttings and inhibit swelling of unstable clay based reactive shale formations. Increasing environmental awareness and resulting legislation has led to the search for new, improved biodegradable drilling fluid components. In the case of additives for clay swelling inhibition, an understanding of how existing effective additives interact with clays must be gained to allow the design of improved molecules. Owing to the disordered nature and nanoscopic dimension of the interlayer pores of clay minerals, computer simulations have become an increasingly useful tool for studying clay-swelling inhibitor interactions. In this work we briefly review the history of the development of technical drilling fluids, the environmental impact of drilling fluids and the use of computer simulations to study the interactions between clay minerals and swelling inhibitors. We report on results from some recent large-scale molecular dynamics simulation studies on low molecular weight water-soluble macromolecular inhibitor molecules. The structure and interactions of poly(propylene oxide)-diamine, poly(ethylene glycol) and poly(ethylene oxide)-diacrylate inhibitor molecules with montmorillonite clay are studied.

The development of effective clay swelling inhibitors is an important goal of the oil and gas exploration industry.
Owing to the disordered nature of clay minerals, and the variability of natural clay composition, laboratory based analysis and characterization of the action of swelling inhibitors in these minerals is highly challenging.It is difficult, if not impossible, to experimentally replicate the interaction of swelling inhibitors with clay minerals under borehole conditions.With recent advances in computational hardware, and the development of increasingly efficient algorithms, computer simulation has become an extremely useful, if not essential tool for understanding the underlying principals behind clay swelling (Bougeard and Smirnov 2007) and for determining how clay swelling inhibitor molecules interact with clay minerals (Bains et al. 2001).
This present paper highlights the effectiveness of molecular dynamics (MD) simulation techniques in the design of improved swelling inhibitors for use in waterbased drilling fluids (WBDFs).We report on results from some recent large-scale MD simulation studies (Greenwell et al. 2005(Greenwell et al. , 2006b)), in which the structure and interactions of low molecular weight, water-soluble poly(propylene oxide)-diamine (PPO-DiAm), poly-(ethylene glycol) (PEG) and poly(ethylene oxide)-diacrylate (PEO-DiAc) inhibitor molecules with montmorillonite (Mmt) clays are considered.
In the next sections of this article we introduce the structure of clay minerals, the mechanisms of clay swelling and briefly consider the evolution of technical drilling fluids.We then survey the literature assessing the environmental impact of drilling operations in the marine environment and some of the legislation affecting drilling fluid design.The concepts of MD simulations are briefly reviewed followed by the application of these techniques to study the three different swelling inhibitor polymers described above.

CLAY MINERAL STRUCTURE
Clay minerals are layered aluminosilicates consisting of stacks of negatively charged two-dimensional layers.Each layer comprises fused sheets of octahedra of Al 3+ oxides and tetrahedra of Si 4+ oxides (Wyckoff 1968).These silica tetrahedra and aluminium octahedra can be arranged in a variety of different ways in clay sheets.Iso-morphic substitutions of the atoms in the tetrahedral and octahedral sites, e.g., Si 4+ by Al 3+ and Al 3+ by Mg 2+ , can lead to an overall negative charge on the clay layers.This is balanced by sorption of charge-balancing cations into the interlayer, and on surfaces and edges of the clay sheet.The interlayer cations are exchangeable and may be replaced with other cations under appropriate conditions.The cation exchange capacity (CEC) of clay minerals depends upon crystal size, pH and the type of exchangeable cation (Pinnavaia 1983).In natural clays there is variability in the charge balancing ions, which tend to be small inorganic species such as Na + and Ca 2+ cations, the charge distribution and structure.

CLAY SWELLING AND ITS INHIBITION IN OIL AND GAS PRODUCTION
In the presence of water, interlayer cations have a tendency to hydrate forcing the clay layers apart.This can occur via two different mechanisms.The first of these is crystalline swelling and is known, both experimentally (Mooney et al. 1952) and through simulations (for example Boek et al. 1995 andHensen andSmit 2002), to occur in a discrete fashion through the stepwise formation of integer-layer or mixtures of integer-layer hydrates.Following this, osmotic swelling may occur in some clay minerals where the concentration of cations present in the interlayer is higher than that in the surrounding water (Norrish 1954).Swelling of this variety leads to substantial increases in the interlayer spacing.The tendency of the sodium smectites, frequently encountered during oil exploration, to swell macroscopically is the principle cause of shale instability.
The type, size and charge of cations present in the interlayer have long been known to greatly impact upon the magnitude of clay swelling (Falconer and Mattson 1933, Lutz 1935, Norrish 1954).Cations with high valences have been shown to be more strongly associated with the surfaces of the clay sheets than low valence cations (Norrish 1954, Young andSmith 2000).They have also been shown to reduce swelling when compared to clays containing lower valence cations (Chávez-Páez et al. 2001).In addition larger, weakly hydrating cations such as K + and Cs + , which can replace interlayer Na + cations, have been shown to significantly reduce swelling (Norrish 1954, Boek et al. 1995, Young An Acad Bras Cienc (2010) 82 (1) DESIGNING DRILLING FLUID ADDITIVES 45 and Smith 2000).Weakly hydrating cations have been shown, through computer simulation, to be reluctant to hydrate fully and act to screen the mutually repelling clay surface more effectively than strongly hydrating cations (Boek et al. 1995).Reduction of clay swelling in oilfield applications has been achieved by this method, however high concentrations of toxic KCl is required, containing about 1 wt.%K + cation, which fail the mysid-shrimp bioassay (O'Brien and Chenevert 1973).Consequently, K + fluids currently find low acceptance for offshore drilling in many waters.
Nevertheless, KCl was commonly used in association with partially-hydrolyzed polyacrylamide (PHPA) as a swelling inhibitor additive in WBDFs throughout the 1960s (Bloys et al. 1994) and its use continues to this day.PHPA helps stabilize reactive shale deposits by coating them with a protective layer of polymer.In the 1970s, the drilling industry switched their attention from WBDFs to the development of oil based drilling fluids (OBDFs) as a means of controlling reactive shales.These fluids proved to be very effective as they had excellent lubricating properties, temperature stability and resulted in minimal swelling.However, the use of OBDFs carries a high environmental price.Even with carefully selected low-toxicity mineral oils, the disposal of treated cuttings still produce a lasting negative environmental impact (Davies et al. 1984).As a consequence, the design and development of environmentally acceptable, biodegradable, WBDFs with OBDF swelling inhibitor properties is currently an area of great interest.
To this end, PEG type polymers have been used in WBDFs to prevent swelling enhancing clay shale stability (Bloys et al. 1994).A relatively low concentration of PEG is required, thus making them relatively cheap to use; and at sufficiently low concentrations they do not significantly affect other important drilling fluid properties such as viscosity and fluid-loss control.High molecular weight (M w ) PEGs also reduce the amount of solids dispersion (non-swelling clays particles) that can occur.Moreover PEGs have a relatively low toxicity, though have a long persistence in the environment.It is thought that PEGs work by three main methods to restrict swelling: i) they diffuse into and displace water molecules from the interlayer.This is thought to be driven through the increased entropy of the resul-tant systems, with the collective non-constrained water molecules having higher degrees of freedom than the polymer replacing them; ii) enthalpic terms, with the energy of the clay-inhibitor interaction driving the reaction.Recent work (Chen and Evans 2005) has suggested this may in-fact be more important than the entropic contributions; iii) At a macroscopic level, clay swelling inhibition by high M w PEG polymers has been attributed to osmotic dewatering of the clay system.However, Smalley in his recent book (Smalley 2006) suggests a bridging flocculation model is far more likely on the basis of neutron scattering data.This is also confirmed by data by Chen and Evans (2004) showing that the d-spacing of clay-PEG systems is irrespective of the M w of the PEG used.
To improve the efficiency of clay swelling inhibition, Coveney and co-workers devised an in situ polymerisation, within the reactive shale, of small soluble molecules with reactive end groups (Stackhouse et al. 2001, Coveney et al. 2004).An example of such a system is the PEO-DiAc system reported here (Greenwell et al. 2006a).In these systems it is important to know how the reactive acrylate groups align between the clay sheets.Computer simulation gives molecular insight that is experimentally unavailable.
There have been several simulation studies of the interlayer behaviour and arrangement of quaternary alkyl ammonium species within clay minerals.Zeng and co-workers (Zeng et al. 2003(Zeng et al. , 2004) ) have shown that, depending on the alkyl ammonium species used, a range of interlayer arrangements including monolayers, bilayers, pseudo-trilayers, and pseudo-quadrilayers of inhibitor molecules may be observed.Pospísil and co-workers (Pospísil et al. 2001) examined the different intercalation behaviour of ammonium surfactants in Mmt paying particular attention to the host-guest interaction energies.Heinz and co-workers (Heinz et al. 2003, Heinz andSuter 2004) have also examined the structural transitions in organo-ammonium clays where the ammonium has long alkyl chains, comparing simulation results with experimental measurements.
A further approach to swelling control, which has received much attention from the experimental chemistry group of Lin and co-workers (Lin et al. 2001, Chou et al. 2003, 2004, Lin and Chen 2004), is to use a water-soluble PPO oligomer terminated with an amine group at each end.It is thought that these amine groups become protonated in the clay interlayer, and as organic hydrophobic cations replace the hydratable metal cations, thus prevent swelling.We have previously studied di-amine systems in a combined experimental and simulation study to understand the effect of organic cation substitution on the structure of the clay interlayer (Greenwell et al. 2005).
Recently, interest has increased in the development of improved and increasingly environmentally friendly swelling inhibition additives for WBDFs based upon a variety of polymers incorporating amine and quaternary alkyl ammonium species (for example, Patel et al. 2007, Kippie and Gatlin 2007, Patel and Stamatakis 2007).
IMPACT OF DRILLING FLUIDS ON THE MARINE ENVIRONMENT Increasing environmental concerns and changes in legislation has limited the use of some drilling fluids with known environmental impacts (Geehan et al. 1990), however knowledge of the impact of many types of drilling fluids is still limited.A recent study (Holdway 2002) assessed the acute and chronic toxic effects of drilling fluids on temperate and tropical marine ecological processes, through evaluation of prior laboratory and field studies.Detrimental effects were indicated on a wide variety of ecosystems including pelagic communities, benthic communities, plankton communities, coral reefs, mangroves and seagrass beds.The environmental response close to a contamination source was identified as a toxic or smothering effect, causing species populations and biodiversity decline, or an enrichment effect in which species abundance increased.
Earlier literature recounts the sublethal effects of chromium or ferrochromium lignosulphate based drilling fluids on 35 species of marine organisms (Hinwood et al. 1994).Close to the contamination source, sublethal effects of ester-based drilling fluids are of moderate concern, indicated by the relatively low toxicity measured in the burrowing marine sand snail, Polinices conicus (Gulec 1994).Marine microbial processes are also adversely affected, with changes in substrate specificity of heterotrophic bacteria, as indicated by studies on subgroups including lipolytic and cellulolytic bac-teria (Okpokwasili and Nnubia 1995).Recent studies have found that whilst dissolved contaminants are expected to dilute in seawater to non-harmful levels, potentially toxic contaminants may accumulate through flocculation, potentially resulting in rapid transport of waste material to the seafloor at levels that can impact benthic organisms (Cranford et al. 1999).The benthic boundary layer has characteristics that differ from the overlying water column that may allow particles from drilling fluids to remain in concentrated suspension close to the seabed (Muschenheim and Milligan 1996).Concentrations of barite as low as 0.5 mg/l (ppm) caused significant effects limiting adult scallop growth of Placopecten magellanicus in the North Atlantic Ocean (Gordan et al. 1992).Studies also indicated that oil droplets sequester contaminating particles, thus affecting levels of contaminants at both the seabed and the surface microlayer (Cranford et al. 1999).Gray et al. (1990) indicated increased abundance patterns of some species and altered presence and absence of rare species due to barium, hydrocarbons and drilling fluids up to 3 km from drilling sites, indicating impacts at a distance from platforms not previously considered.Daan et al. (1990Daan et al. ( , 1994Daan et al. ( , 1995Daan et al. ( , 1996) ) indicates significant effects of OBDFs and ester based drilling fluids (EBDFs) on macrobenthos species abundance up to 1 km from the drilling site, with the greatest effects shown in the echinoderm Echinocardium cordatum and its symbiont, Montacuta ferruginosa.WBDFs commonly indicated a lack of significant environmental effect in comparison to OBDFs and EBDFs for most species.However, whilst the abundance of most of the indicator species was reduced near OBDFs and EBDFs sources, the abundance of the opportunistic polycheate Capitella capitata was shown to increase indicating organic enrichment owing to changes in sediment properties.Larval settlement of red abalone, Haliotis rufescens, in addition to adult mortality, tissue loss and viability of brown cup coral, Paracyathus stearnsii, were shown to be adversely affected by even low toxicity WBDFs of 0.002-200 mg/l (Raimondi et al. 1997).In addition, whilst components of WBDFs and EBDFs are not considered particularly toxic, they can have a low biodegradability and therefore may persist for long times in the environment.Few studies take into account An Acad Bras Cienc (2010) 82 (1) the increased toxicity of bioaccumulation and food chain transfer and its consequential effects at higher trophic levels (Neff et al. 1989, Frost et al. 2002).
In summary, there is a paucity of data on both the acute and chronic effects of drilling fluids and their associated wastes on many marine organisms and ecological processes.Whilst the impacts of drilling fluids on multiple ecological processes have been highlighted, the inherent variation in different natural environments coupled with variations in temporal and spatial scales have so far impeded the ability to predict the long-term impacts of drilling fluids on the marine environment (Holdway 2002).Design of suitable drilling fluid additives can lead to low environmental impact, waste minimization through recycling, and reuse of fluids and their components.

COMPUTER SIMULATION FOR STUDYING CLAY-ORGANIC SYSTEMS
In many papers concerning polymer-clay systems, the only form of analyses presented are powder X-ray diffraction (PXRD) and, to a lesser extent, thermo-gravimetric analysis (TGA), which reveal the interlayer spacing and organic/water content respectively.Fourier transform infra-red spectroscopy (FTIR) has also been used to study the adsorption of low M w amines by clays, for example to demonstrate the presence of Brønsted acid sites in the clay interlayer (Yariv and Cross 2002).These conventional analytical techniques, whilst useful, are by themselves incapable of elucidating the interlayer arrangement, let alone dynamical behaviour.However, the insight gained by these methods provides useful starting data for building models for computer simulation.Computer simulation is a valuable tool for gaining understanding of the arrangement and dynamics of intercalated molecules in clay systems (Greenwell et al. 2006b).The complexity and size of inhibitor-clay systems generally rules out the use of quantum mechanical (QM) based methods other than to understand the interactions between reactive groups and small components and clay sheets.Therefore, simulations are typically carried out using classical force-field based methods.
Here interactions between all of the component parts of the system are based upon combinations of empirical force-field terms, which describe bond vibrations, an-gles, torsions, non-bonded van der Waals and coulombic interactions.

MOLECULAR DYNAMICS
By applying initial velocities to a configuration of atoms and solving Newton's equations of motion the potential energy surface of a system may be traversed in a deterministic fashion and the temporal evolution of a system followed.This is known as molecular dynamics (MD).In this technique thermal energy is included using a thermostat, which allows potential energy barriers to be overcome, in a realistic manner.The main advantage of the method is that the dynamical evolution of a system, with time, may be followed, which allows comparison with additional experimental techniques such as NMR and quasi-elastic neutron scattering.The data from these simulations provides precise information regarding the coordinates of all atoms within the model at any point in time during the simulation period.This allows the interlayer arrangement and dynamics of organic and water molecules to be evaluated with equal precision.

MODELLING PERIODIC SYSTEMS
In order to model the bulk structure of materials (greater than 10 23 atoms) using relatively small models (generally less than 10 4 atoms), two methods are often employed: (i) The use of super-cells, where the original unit-cell, usually derived from a crystal structure, is replicated several times and then redefined as one larger simulation cell, (ii) periodic boundary conditions are applied to the simulation cell, where the super-cell is considered to be replicated infinitely in all three orthogonal space directions.It is assumed that such an approach will replicate the properties of the extended condensed phase system.

SIMULATION METHODOLOGY
Each of the model inhibitor-clay systems was prepared as described in our previous work in this area (Boulet et al. 2003, Greenwell et al. 2005, Thyveetil et al. 2008).
In this section we briefly recount some of the main simulation details.I. To investigate finite-size effects, one of the initial models was further replicated to form a super-cell containing 350840 atoms.This was simulated for 0.5 ns.A total of four Mmt based systems were simulated consisting of PEG (M w = 414 g mol -1 , see Fig. 1b) and PEO-DiAc (M w = 258 g mol -1 , see Fig. 1c), intercalated within Na + and K + Mmt galleries.In these simulations models were constructed in the order of 20,000 atoms.The unit cell composition of each model is given in Table I.In order to measure the materials properties of these inhibitor-clay systems, very large-scale simulations are needed to sample the long wavelength undulations required to extract the bending modulus.The undulations are only noticeable for large system sizes; for Mmt, the collective motion of the clay atoms only manifests as long-range undulations at length scales greater than 100Å.We have therefore constructed two very largescale Mmt models through replication of much smaller unit cells: i) intercalated with a double layer of water, with a system size of 1055000 atoms, ii) intercalated with a double layer of PEG (M w = 3978 g mol -1 ) with a system size of 1756020 atoms.

SIMULATION CODE AND PARAMETER DETAILS
The MD simulations were performed using the Largescale Atomic/Molecular Massively Parallel Simulator, LAMMPS (Plimpton 1995).LAMMPS uses algorithms and techniques that allow the code to exhibit a near linear relationship (scaling) between the number of processors used, or the size of the model system of interest, and the time taken for the simulation to be performed (Hein et al. 2005, Suter et al. 2007).To model the various systems two different parameter sets were required: (i) The Teppen force-field (Teppen et al. 1997), (an extension of the cff91 force-field) which has been specifically parameterized to account for the behaviour of clay minerals, as well as describing the organic molecules was employed to simulate all models > 500000 atoms in size.This force-field has been validated for the swelling behaviour of the Na + -Mmt clays (Hein et al. 2005) and used to study the interlayer arrangement and dynamics of all of the PPO-DiAm-Mmt and PEO-DiAc-Mmt inhibitor-clay systems, as well as the smaller (< 25K atoms) PEG-Mmt models (Boulet et al. 2003, Greenwell et al. 2005, 2006a, b).(ii) To calculate the elastic properties of inhibitor-clay systems a completely flexible force-field is required.To this end the more recently developed ClayFF force-field (Cygan et al. 2004) was used to describe the clay framework of the much larger models (> 1000000 atoms).This force-field treats most interatomic interactions as non-bonded Coulombic and Lennard-Jones (12-6) pair potentials, which allows complete translational freedom and permits simulation of complex disordered systems.Intercalated organic inhibitor molecules were described by CVFF force-field parameters which are compatible with the ClayFF forcefield.This combination of ClayFF and CVFF has been  (Liu 2007).
In all cases water was represented by the SPC model (Berendsen et al. 1981) which is consistent with the ClayFF, CVFF and Teppen force-fields.Short-range non-bonded interactions were described using a Lennard-Jones 9-6 potential where the Teppen force-field was used and a 12-6 Lennard-Jones potential for the larger systems simulated using the ClayFF and CVFF force-fields.In both cases these interactions were assigned a cut-off distance of 10Å.Calculation of the computationally demanding long range interactions was carried out using the PPPM algorithm (Hockney and Eastwood 1981, Luty et al. 1994, Toukmaji and Board 1996).The rRESPA algorithm was employed to economize simulation time, with an outer time step for electrostatics of 4 fs and an inner time step for bond terms of 0.5 fs.Three-dimensional periodic boundary conditions were applied to all models to represent the bulk material.Models were constrained to be orthogonal by the limits of the LAMMPS package.The inhibitor-clay structures were optimized to remove repulsive interactions and then simulated at 300K and 1 atm using an isobaric-isothermal (NPT) ensemble in which the temperature was controlled via a Langevin thermostat.

COMPUTATIONAL DETAILS
Simulations of the inhibitor-clay systems were performed on a variety of high-end computing resources.In typical recent simulations, a federation of three supercomputing grids was employed to perform the simulations as well as to visualize the data produced; since this is a somewhat novel computing infrastructure we also briefly describe pertinent aspects of it in this section.Simulations were performed on the following federated grid infrastructures: the UK's National Grid Service (www.ngs.ac.uk), including HPCx (www.hpcx.ac.uk), the TeraGrid (www.teragrid.org)and the EU Distributed European Infrastructure for Supercomputing Applications (www.deisa.org).Submission of jobs was facilitated by the Application Hosting Environment (AHE) (http://www.omii.ac.uk/).The AHE allows the submission of geographically distributed jobs through a single uniform interface which interoperates between Globus and Unicore grids and also retrieves output data automatically once a simulation has finished (Coveney et al. 2007a, b).

ANALYSIS AND VISUALIZATION
Analysis of the potential energy of each simulation showed that all of the model systems were equilibrated within the first 0.2 ns of simulation and data collection An Acad Bras Cienc (2010) 82 (1) was started after this period.Data was collected every 0.1 ps and statistical averages were evaluated over the last 800 ps of simulation where possible.Time averaged radial distribution functions (RDFs) and 1-D atom distribution plots were calculated for selected atom types.Visualization was achieved using the UCSF Chimera package for models containing less than 10,000 atoms.Larger models were visualized using the AtomEye software package (Li 2003), which has more efficient rendering capabilities for larger systems and can be run on parallel processors.Parallel rendering of the trajectories was performed on six processors of an SGI Prism, further reducing time spent on visualization.Fast networks are required in order to transfer data from supercomputing resources for post-processing and visualization.Switched optical networking was used in order to transfer the large quantities of data generated in this study, ∼ 250 GB. (Coveney et al. 2007b) Specifically, the JANET lightpath research network was utilized (http://www.ja.net/services/lightpath/).This service operates within the UK and connects to similar network infrastructures across Europe and the U.S.

RESULTS
In this section we first present the results obtained from our studies of Na + Mmt with PPO-DiAm swelling inhibitors, and examine the results of increasing protonation of the amine group.This is followed by the results from studies of PEG swelling inhibitors in Na + Mmt, the way PEG effects plasticity and other material behaviour of the clay-water-inhibitor system, and then the effect of functionalizing the PEG with DiAc.We next look at the effect of exchanging the Na + cation for K + cations as might be the case where KCl brines are used to help suppress swelling.By examining the interactions at a molecular level the subtle interactions between inhibitor and clay and effect on water arrangement can be deduced, vital for successful design of drilling fluid composition.

Interlayer arrangement and bonding in PPO-NH 2 inhibitor-clay systems
The PPO-NH 2 molecules in the simulated 7160 atom system formed a monolayer with a d-spacing of 13.98Å (±0.002Å).This is in good agreement with the experimental value of 14.3Å (Greenwell et al. 2005).Figure 3 shows the 1D atom density distribution across the interlayer regions for the amine/ammonium group N and H atoms.The PPO-NH 2 N atoms are located slightly offset either side of the mid-plane of the interlayer.Little or no H-bonding occurs between the amine groups and the siloxane surface oxygen atoms in the tetrahedral layer of the clay sheet.The other atoms in the PPO backbone are observed to lie along the mid-plane of the interlayer region.The interlayer water molecules in this system are arranged almost parallel with the clay sheets.The water O atoms are arranged near the interlayer midplane, but oriented towards the face of the clay sheets with the water H atoms residing close to, but slightly offset from, the mid-plane of the interlayer.As in all the systems we have investigated, Na + cations were arranged mainly adjacent to the face of the clay sheet, sometimes penetrating into the tetrahedral layer slightly.The RDF for this system showed that the Na + cations are closely coordinated by the O atoms of the PPO backbone (2.6Å), the interlayer water (2.4Å), the clay sheet (2.6Å), and by the N atoms of the amine groups (2.5Å).A PPO-aminosodium cation coordination has also been postulated by Lin et al. (2004), who suggested that the combined O and N atom coordination to Na + provides a driving force for adsorption of the organic molecules.

inhibitor-clay systems
The simulated d-spacing for the 7160 atom system ammonium system was 14.39Å (with a standard deviation of ±0.001Å), in reasonable agreement with the experimental work of both ourselves (Greenwell et al. 2005) and others (Lin and Chen 2004), but approximately 0.5Å more than for the simulated PPO-NH 2 system.The difference can be rationalized by considering the atom distribution within the interlayer (Fig. 3) and comparing it to the distribution for the PPO-NH 2 model.3 ) groups to locate adjacent to the negatively charged clay sheets, also increasing H-bonding interactions between the ammo-nium H atoms and the O atoms of the silicate sheet.This conformational change does not occur when only amine (-NH 2 ) groups are present.
The PPO-NH + 3 backbone C and O atoms, Figure 3b, are arranged along the mid-plane of the interlayer region, with the methyl groups and the O atoms both slightly offset either side of the mid-plane, due to the staggered nature of these groups in the PPO backbone.The Na + cations adopted positions along the face of the clay sheets (Fig. 3c), in some cases interpenetrating the rings formed by the siloxane surface oxygen atoms of the tetrahedral aluminosilicate layer.In the case of the ammonium-intercalated systems with very few Na + cations the interlayer water adopted an arrangement close to the faces of the clay sheets, with the water H-atoms oriented predominantly closer to the face of the clay sheet than the water O atoms.This arrangement of water molecules is in contradistinction to the PPO-NH 2 systems where the water O atoms were arranged near to the interlayer mid-plane and oriented co-planar with the clay sheets.
A more detailed analysis of the local environment about each atom type is given by the RDFs (not shown).The ammonium cation H atoms are coordinated strongly by interlayer water (1.8Å inter-nuclear separation), indicating that the ammonium ion behaves somewhat akin to a Na + cation.The ammonium H atoms also strongly interact with the siloxane surface O atoms of the clay sheet at a distance of some 2.4Å, as suggested by the 1-D atom distributions.This arrangement of the positively charged ammonium groups is again similar to the similarly charged Na + cations, which also interact closely with the O atoms in the PPO backbone.The ammonium group is oriented such that the H-bonding and favourable electrostatic interactions are maximized.In conditions of low interlayer water content, Na + cations, by comparison, tend to interpenetrate slightly into cavities within the tetrahedral layer of the clay sheet so as to maintain a full coordination shell of O atoms made up partly by water, partly by clay.
The large-scale, 350840 atom system, simulated for 0.5 ns, showed an average simulated d-spacing of 14.40Å (±0.001Å).The 1-D atom maps were found to differ slightly from the corresponding smaller model, the distribution of atom types retained the same broad features in both systems, but the atom density in the larger system was found to be less constrained than in the smaller model.Visualization showed that some undulations in the clay layers had formed (see Greenwell et al. 2005).It is likely that this effect is due to the large supercell employed; small periodic models are much more tightly constrained by symmetry to rigid clay sheets.In short, absence of such undulations is due to finite size effects.The observation of such clay sheet undulations are of significance as they allow the calculation of materials properties, vide infra.

Interlayer arrangement and bonding in PPO-NH + 3 /PPO-NH 2 inhibitor-clay systems
The simulated PPO-NH + 3 /NH 2 systems had average dspacings of 13.69Å (±0.001Å) and 15.15Å (±0.001Å) for the 33% ammonium and 66% ammonium systems respectively.The larger average d-spacing for the latter case arises due to one of the interlayers in the model being significantly more expanded (ca.15.9Å) than the other (ca.14.4Å), the reason for which is described in more detail below.
An examination of the atom density distribution across the interlayer shows that the ammonium groups in both PPO-NH + 3 /NH 2 systems (Fig. 3) behaved similarly to those in the PPO-NH + 3 system, being found adjacent to the clay sheet siloxane surface O atoms as can also be seen in Figure 2. In contrast to the PPO-NH 2 systems, the amine groups in the 66% ammonium system were further away from the mid-plane of the interlayer, with the bulk of the amine H atom density located towards the positively charged ammonium groups, suggesting possible H-bonding between these species, which is manifested in the RDF plot discussed later.The amine groups in the 33% ammonium system were distributed more in line with the amine groups in the PPO-NH 2 system (Fig. 3).The monomer backbone C and O atoms were distributed in the mid-plane of the interlayer, with the exception of the much wider interlayer noted above.In this larger interlayer phase, of d-spacing ca.16Å, the PPO appeared to form a bilayer, or pseudo-bilayer, arrangement.The pendant methyl groups formed a pseudo-trilayer and the O atoms were broadly distributed.Visualization of the monomers An Acad Bras Cienc (2010) 82 (1) showed that this was due to intra-molecular H-bonding causing a coiled monomer conformation to arise, resulting in an apparent bilayer arrangement.
The distributions of Na + cations and water across the interlayer are similar to those of the PPO-NH + 3 system (Fig. 3) for the 66% ammonium system, and followed the PPO-NH 2 system for the 33% ammonium system.This illustrates that the distribution of interlayer Na + and water varies according to the number of amine groups that have been protonated.So far as the intercalated organic molecules are concerned, in scenarios where there are predominantly ammonium groups the interlayer adopts an arrangement similar to the case where there are all ammonium groups, and vice versa for amine groups and compounds.The RDFs for the amine N atoms (not shown) indicate that the ammonium H atoms approach within 1.9Å on average, indicating strong H-bonding.The strong nature of such H-bonding in the simulated system would also be a plausible explanation for the high shift noted for the N-H bending mode absorption (indicative of H-bonding) in the PPO-NH 2 and PPO-NH + 3 /NH 2 systems FTIR spectra (Greenwell et al. 2005), suggesting that in the PPO-NH 2 system a mixture of H-bonded -NH 2 and NH + 3 groups exist in the interlayer, rather than just -NH 2 .
PEG BASED NA + -MMT INHIBITORS Simulations revealed average d-spacings of 1.70nm for the Na + -Mmt-PEG system.The 1-D atom density distribution for the PEO backbone C atoms shows two distinct peaks corresponding to a definite bilayer arrangement for PEG, concentrated within a relatively narrow region about the mid-plane.In the 1-D atom density distribution plot, the Na + cations are found embedded in the faces of the clay sheets, with a very small number of Na + cations further out into the interlayer region adjacent to the face of the clay sheets (Fig. 4).The water molecules were found to hydrate the exposed surfaces of the Na + cations studding the clay surface and pointing toward the interlayer space.As such the water molecule distribution follows the asymmetry present in the Na + 1-D atom density distribution and the formation of hydration spheres is apparent from the Na + /water RDF.Due to the presence of the Na + cations and associated water molecules at the face of the clay sheets, the organic molecules were further from the clay sheet resulting in the more compact bilayer.Na + ions interacted with the PEG hydroxide groups and the PEG backbone O atoms.

Bulk properties of PEG systems
Intercalation of polyethers such as PEG can not only affect the swelling properties of clay but can also change the elastic and viscoelastic properties of the clay material.However, experiment measurements of the elastic moduli of smectite clay platelets have, thus far, not been successful (Chen and Evans 2006), so it is not possible to calculate values of modulus to use in composite theory; inhibited clay systems can be considered as high clay fraction clay-polymer composites.To calculate the material properties of Mmt, in Figure 5 we show the spectral intensity per undulatory mode versus wave-vector in the x and y directions for Mmt with water and with 3978 g mol -1 PEG polymer, with system sizes of 1055000 and 1756020 atoms respectively (Table I).To calculate the elastic properties of the Mmt sheet, we relate the wavelength and amplitude of the thermal undulations of the clay sheet to its material properties, such as the bending modulus, k (Suter et al. 2007): where h is the height function of the clay sheet, q is the wavevector of the undulation, A is the area of the clay sheet, k B is the Boltzmann constant and T is the temperature.We find through fitting the long wavelength behaviour to a q −4 fit that the bending modulus in both x and y directions to be 1.7 × 10 -17 J, for both clay systems.The identical bending moduli indicates that the undulations of the clay are unaffected by the interfacial medium.This is likely to be due to the very large in-plane elastic modulus of the clay sheet, which is far greater than the intercalated medium.

DiAc functionalized PEO inhibitors
Simulations revealed average d-spacings of 1.59nm for the Na + -Mmt-PEO-DiAc systems respectively.The 1-D atom density distribution for the PEO backbone C An Acad Bras Cienc (2010) 82 (1) atoms has two distinct peaks corresponding to a definite bilayer arrangement for the Na + -Mmt-PEO-DiAc.The carbonyl O atoms on the PEO-DiAc were found, in part, to have a similar spatial distribution to the PEO C atoms, but also occupied much of the space along the mid-plane of the interlayer suggesting cross-linking between the acrylate end groups in different monolayers may be possible.A similar distribution occurred for the PEG hydroxyl groups.The RDFs also show that the Na + also interacted somewhat with the PEO-DiAc carbonyl O atoms and the PEO chain O atoms of the PEO-DiAc.
THE EFFECT OF EXCHANGING NA + FOR K + IN MMT SYSTEMS Simulations revealed average d-spacings of 1.69nm and 1.84nm, respectively, for the K + -Mmt-PEO-DiAc and the K + -Mmt-PEG system.The larger d-spacing of the K + -Mmt-PEG was found to correspond to a clear trilayer arrangement of organic molecules (Fig. 4) compared to the bilayer observed for the K + -Mmt-PEO-DiAc.In the K + -Mmt-PEO-DiAc the K + cations were distributed only adjacent to the faces of the clay sheet.This suggests that K + has better mobility within the interlayer, due to lower charge density resulting in lower interaction with water and monomers, and is able to migrate to the face of the clay sheet.For the K + -Mmt-PEG, a small amount of K + atom density was observed towards the mid-plane of the interlayer (Fig. 4).Analysis of the 1-D atom density distribution plots and RDFs for the water O atoms, the PEO C atoms and the PEO-DiAc carbonyl O atoms showed that the K + cation was approached as closely by the organic molecule as the water molecules, and in fact the carbonyl O atoms in the PEO-DiAc were closer to the face of the clay sheets than the water O atoms.The fact that the organic molecule approaches the cation as closely as the water supports the suggestion that the K + cation has less interaction with water due to its lower surface charge density relative to Na + .

DISCUSSION
In the PPO-DiAm systems, the interlayer spacing under the controlled conditions of computer simulations, for systems that differ only in the number of amine/ammonium groups and Na + cations, was found to be depen-dent on the number of intercalated ammonium groups, and showed good agreement with the experimentally observed monolayers at similar organic and water loadings.Computer simulation shows that the organic molecules studied in this work are generally arranged in a monolayer within the interlayer.The PPO backbones of the molecules are arranged along the mid-planes of the interlayer and the orientation of the headgroups depends on whether an ammonium or amine group is present.The arrangement of the PPO backbone contrasts with the behaviour of PEO based polymers (as can be seen in the following section), which are hydrophilic and arrange themselves along the face of the clay sheets to form bilayers (Boulet et al. 2003).
As discussed, in MD simulations the amphiphilic nature of the Na + -Mmt is influenced by the presence of organo-ammonium species, the orientation of water molecules within the interlayer depending upon whether the monomers were terminated with ammonium or amine groups, due to the changes in H-bonding networks and absence of Na + cations which might otherwise coordinate the water strongly.
The range of the calculated d-spacings for the Na + systems showed good agreement with experimental values, given the differences in composition and uncertainty in the degree of intercalated material, typically within about 5%.The ability of the Na + cations to diffuse into the tetrahedral pockets of the clay surface has already been observed both experimentally (Yang and Zax 1999) and theoretically (Hackett et al. 2000) for Li + and Na + cations.The absence of the cations and their associated hydration spheres, from the interlayer region, results in both a more organophilic interlayer region and more space for organic molecules for a given d-spacing.
The simulated K + -Mmt-PEO-DiAc systems showed particularly poor agreement to experimental d-spacing.However, it should be noted that the experimental XRD reflections are very broad for K + -Mmt-PEO-DiAc systems (Greenwell et al. 2006b).Paradoxically, these systems containing K + cations, which have been shown to be less able to migrate into the tetrahedral layers of the clay sheets, sometimes have very high organic loadings based on TGA data (Greenwell et al. 2006b).This can be accounted for through the experimental observation that these materials tend to exfoliate and hence exhibit domains where single clay sheets are not assembled in tactoids, but are rather dispersed throughout the polymer matrix.Exfoliation accounts for the lack of agreement between simulated and experimental d-spacings for the K + -Mmt systems, since exfoliation cannot be simulated in these strictly periodic models.The propensity for K + -Mmt to exfoliate is noteworthy as it is in contrast to the known resistance of the K + -Mmt to water swelling (Boek et al. 1995) and illustrates the need to understand mineral-organic interactions when designing drilling fluids; the use of K + salt inhibitors with certain organic inhibitors may lead to dispersion and destabilization of clay fractions.For all the simulated PEG systems no evidence was found for H-bond interactions between the protons of the PEG alcohol groups and the tetrahedral O atoms of the clay surface.Therefore, it seems that in the presence of water and cations, PEG is unlikely to adsorb strongly at the clay surface.The PEO chains for both PEO-DiAc and PEG tend to orientate with the O atoms towards the mid-plane for the Na + clays, away from the cations which reside at the clay sheet surface.This arrangement, which results in organic monomer C atoms adjacent to the organophilic silica surface, has been reported previously by others (Bujdak et al. 2000).
Fig. 5 -Spectral intensity per undulatory mode versus wave-vector in the x and y directions for MMT intercalated with water (1055000 atoms) and 3978 g mol -1 PEG (1756020 atoms) in the x and y directions.The dashed line is a fit to the undulatory q −4 modes for q < 0.04 with k = 1.7 × 10 −17 J.The local peak centred on q = 0.35Å −1 and 0.175Å −1 is due to artificial periodicity resulting from the way isomorphic substitutions are included in all the models.Isomorphic substitutions in real clays are randomly distributed so this peak can be disregarded.
The choice of functionalized PEO was also found to affect the cation distribution across the system interlayer.In the PEG systems hydroxyl (alcohol -OH) groups retained some of the cations and associated hydrations shells within the interlayer region.The magnitude of this effect is dependent upon the cation, with nearly all the K + migrated to the face of the clay sheet.Since the cations are retained in the interlayer region away from the clay sheet surface they are also closely associated with the monomer backbone O atoms.Therefore in the RDFs, the order of interaction for both the PEG hydroxyl O atoms, and the backbone O atoms, with the cations is: Na + > K + .
Conversely, the PEO-DiAc, having no alcohol OH groups, did not retain the cations in the mid-plane of the interlayer region; instead, the vast majority of the Na + cations migrate into the tetrahedral layer of the clay sheet, and the K + cations migrate to the face of the clay sheets.Therefore, comparing the RDFs for the inter-action between the different cations and the PEO-DiAc backbone O atoms, or the endgroup O atoms, showed the strength of the interaction with the low surface charge density cations decreasing in the order: K + > Na + .
Furthermore, both the resistance to swelling and the tendency of K + -Mmt to exfoliate can be rationalized in light of the fact that the low surface charge density K + cation is not particularly hydrophilic, and sheds its hydration shells (Boek et al. 1995).This both renders the clay more resistant to swelling in an aqueous environment and, if the driving force for intercalation is entropy-favoured displacement of interlayer water, facilitates the uptake of organic molecules, when present, by allowing water of hydration to be exchanged (Bains et al. 2001).

CONCLUSIONS AND FURTHER WORK
We have reported on the use of large-scale MD on high performance computing facilities, and using grid-computing methods, to undertake simulations of clay swelling inhibitor systems.The results allow interpretation of bulk experimental data of the mode of interaction of swelling inhibitors with clay minerals at an atomistic and molecular level.This highlights the use of computer simulations for the design of new and improved drilling fluids.The ultimate goal of computer aided design is to allow the oilfield industry to produce improved swelling inhibitors for reactive shales that maximize oil and gas recovery and reduce costs associated with wellbore instability.With increasing environmental awareness, the design of biodegradable fluid components, including swelling inhibitors, with similar efficiency to current technologies is a key priority.Furthermore, by using very large-scale simulations, enabled through grid computing, material properties can been calculated for hydrated clay systems, which are in agreement with the experimental data available for related systems.Future work will use the techniques and understanding of swelling inhibitor-clay interactions described here to design, in silico, new swelling inhibitors with improved environmental performance, i.e., increased biodegradability and consequently lower residence time in the ecosystem.
MODEL DETAILSThe model clay used was a Wyoming-like Mmt with chemical formula [Al 15 Mg 2 ][Si 32 O 96 H 16 ]M + (M + = Na + , K + ).Four possible scenarios of PPO-DiAm (M w = 248 g mol -1 , see Fig.1a) Na + -Mmt systems using experimental loadings of inhibitor and water have been simulated: (i) None of the Na + was exchanged and all the amine was unprotonated (-NH 2 ); (ii) 33% of the (-NH 2 ) and Na + was exchanged by ammonium species (-NH +3 ); (iii) 66% of the Na + was exchanged by ammonium species; (iv) all the monomer was protonated ammonium (-NH +3 ), equivalent to exchange of 83% of the Na + cations.The number of inhibitor intercalants and water molecules remained constant for each model, based upon experimental TGA data(Greenwell et al. 2005) for the PPO-NH 2 system and corresponded to 10% mass inhibitor content and 1% mass water.The super-cell composition of each model is given in Table

Fig. 2 -
Fig. 2 -Snapshots after 1 ns of molecular dynamics simulation of the 7160 atom 33% ammonium mix (PPO-NH 2 /PPO-NH + 3 ) system showing the interlayer arrangement of ammonium cations and adjacent clay sheets.The ammonium cation is arranged to maximise H-bond and electrostatic interactions.When compared to the Na + cation, the ammonium cation is unable to sit closer to the cavities in the tetrahedral layer of the clay sheet due to steric restrictions and strong H-bond interactions with surface siloxane surface O atoms.The colour scheme is: C gray, H white, O red, N blue, Si orange, Al green, Na + brown.
Figure 3a shows the distribution of the ammonium group N and H atoms in comparison to the siloxane surface O atoms of each clay layer.The ammonium groups are predominantly arranged adjacent to the face of the clay sheets, with the ammonium H atom density closer (ca.1.5Å) to the siloxane surface O atoms than the ammonium N atom density (ca.1.9Å), indicating the formation of a domain of ammonium groups H-bonded to the aluminosilicate sheets.The slightly expanded interlayer in the PPO-NH + 3 system arises due to conformational changes in the molecules, discussed below, to allow the positively charged ammonium (-NH +

Fig. 4 -
Fig. 4 -Snapshots after 1 ns of molecular dynamics simulation for (a) the PEG-K + Mmt composite and (c) the PEG-Na + Mmt composites.The one-dimensional atom density maps of the change in K + cation distribution for the simulated PEG and PEO-DiAc K + Mmt composites is shown in (b) and the change in Na + cation distribution for the simulated PEG and PEO-DiAc Na + -Mmt composites in (d).The x-axes show the distance from the centre plane of the interlayer in nm (×10).The y-axes show the arbitrary density.Dashed line = poly(ethylene glycol) and solid line = poly(ethylene oxide) diacrylate for comparison.

TABLE I Simulation cell composition of clay-polymer molecular models. For the mixed PPO-NH 2 /NH + 3 systems the percentage NH + 3 content is labelled.
(Padma Kumar et al. 2006cribe the intercalation of organolayered double hydroxides(Padma Kumar et al. 2006) and alkylammonium-intercalated smectites