Applications of Molecular Simulation in Oil and Gas Production and Processing

Applications of Molecular Simulation in Oil and Gas Production and Processing — Molecular simulation is an emerging technique which consists in performing a detailed simulation of microscopic systems involving typically a few hundreds of molecules. On the basis of these simulations, appropriate statistical averages are performed to derive fluid properties that can be compared with experimental measurements. The purpose of the article is first to provide the reader with basic notions of molecular simulation. Then, application examples are given in several fields of the oil and gas industry. In reservoir engineering, the examples relate to the properties of poorly known hydrocarbons, the thermal properties of natural gases, and the phase equilibria and volumetric properties of acid gas hydrocarbon mixtures. The application to gas production and processing is illustrated by the phase equilibria involving methanol (a common solvent or hydrate inhibitor) and the solubility at high pressure of gases in polymer materials. In hydrocarbon processing, the solubility of hydrogen sulfide in hydrocarbons at high temperature is discussed. Petroleum Industry Applications of Thermodynamics Applications de la thermodynamique dans l'industrie pétrolière D o s s i e r Oil & Gas Science and Technology – Rev. IFP, Vol. 61 (2006), No. 3


INTRODUCTION
The essence of molecular simulation is to consider explicitly all the molecules of a small size system and to determine the behaviour of the latter from a careful computation of the interactions between its components.Moreover, the position of most atoms must be explicitly computed, in order to account for the difference in geometry between individual molecules (Fig. 1).These computations take much more computer time than classical thermodynamic models: from a few hours to several weeks of computing time, depending on system complexity.Yet, these methods are increasingly used in the oil and gas industry to compute thermodynamic properties.Why is this so?There is no short answer, and several kinds of reason must be quoted.Firstly, molecular simulation is the unique way to gather the prediction of all thermophysical properties in a single theoretical framework, which is statistical thermodynamics [1-3].It is sufficient to provide a receipe to compute the potential energy of the system from the positions of its atoms or molecules, and all the system's properties can be derived by statistical thermodynamics.This holds for equilibrium properties like vapour pressure, as well as for transport properties like viscosity [4][5][6].Of course, finding the adequate receipe for computing the potential energy without exceeding reasonable computing times is not an easy task, and the existing receipes will certainly be improved in the years to come.Nevertheless, this field is sufficiently mature that it can provide reliable predictions in many instances, as will be illustrated in the present article.In these cases, molecular simulation has shown an unprecedented ability to extrapolate prediction from small molecules to large ones, from low pressures to high pressures, and from low temperatures to high temperatures.
The second major reason to consider molecular simulation is the difficulty of modelling many systems with classical thermodynamic models (equations of state, group contribution methods, activity coefficient models, etc.) in absence of experimental data.Despite the large amount of experimental data in the literature, the oil and gas industry faces several problems where this kind of difficulty is encountered: -systems containing toxic or corrosive components like hydrogen sulfide or organo-mercuric compounds; -heavy hydrocarbons, which are not commercially available as pure chemicals to perform measurements; -unstable systems at process temperature; -properties of microporous adsorbents prior to their synthesis; -systems at extreme pressures.
In such cases, the extrapolation capacity of molecular simulation makes it the safer way to predict thermophysical properties.This is especially the case when specific interactions occur at the molecular scale, as in zeolitic microporous adsorbents [7], or when polar systems are considered.
In the present article, we will first give a brief introduction to molecular simulation methods.Then we will provide an overview of application examples in reservoir engineering, gas production and hydrocarbon processing.However, we will restrict these examples to fluid phase properties.It must be stressed that molecular simulation of fluid phases is a vast field, that is impossible to treat in detail within such an article.We will therefore frequently refer the readers to articles or manuals where they will find complementary information.

INTRODUCTION TO MOLECULAR SIMULATION METHODS
Molecular simulation comprises mainly quantum chemistry methods and statistical thermodynamic methods, the first being based on the fact that real systems can access only 388 For each type of problem, molecular simulation has provided useful predictions together with a better understanding of the relation between properties and chemical structure.It provides an intermediate way between experiments and classical thermodynamic models.Cheaper and more rapid than new measurements, it nonetheless allows an improved determination of the parameters of routine models.
discrete (or quantized) energy levels.For the purpose of the present article, we will focus on the methods based on statistical thermodynamics, in which the classical approximation allows to consider that all energy levels can be accessed by the system [1].

Statistical Ensembles
The concept of statistical ensemble is the cornerstone of statistical mechanics.To put it simply, it consists in a collection of all possible states of a system compatible with some set of constraints, such as imposed temperature or imposed number of molecules [1].From this collection of snapshots, average properties can be determined.Depending on the application desired, different statistical ensembles will be used (Table 1).Each statistical ensemble is characterized by its probability density, i.e. the probability of occurrence of each system state in the collection.An important specific aspect of statistical ensembles is the fluctuations of the variables that have not been imposed.Let us consider the example of the canonical ensemble, where temperature, number of molecules and volume are imposed.Its probability density is proportional to exp(-β E i ) where β = 1/kT and E i is the energy of the system state i, k = 1.381 10 -23 J.K -1 being the Boltzmann constant.In classical thermodynamics, energy would be imposed to the system.In the canonical ensemble, the energy is fluctuating, i.e. the collection contains system states at different energies.However, we can compute the average energy, noted 〈E〉: (1) where m is the number of system states in the ensemble.The energy of a particular state of the system does not teach us anything meaningful, but the average energy can be compared with experimental results, for instance through the molar energy: (2) where N A = 6.02 10 23 is Avogadro's number and N the total number of molecules of the system.
Another useful statistical ensemble is the isothermalisobaric ensemble, or NPT ensemble, where pressure is imposed instead of volume.This ensemble is used for instance when the density of a gas mixture is to be determined at known pressure and temperature.Its probability density is proportional to exp(-βE i -βPV i ) where V i is the volume of system state i.As volume fluctuates, the volume of a given state does not obey an equation of state, but the average volume does.
The ensemble that is most adapted to adsorption in microporous solids is the Grand Canonical ensemble, where the temperature T, the volume, and the chemical potential μ j of each species j are the imposed variables.Its probability density is proportional to where N ji is the number of molecules of type j in system state i.For applications to the adsorption on solids, the chemical potentials can be determined from the partial pressure of the components (among other possible ways) and the energy computation must account for the interaction energy between each molecule and the microporous adsorbent.The most desired output is the average number of molecules of every species in the system, i.e. the amount adsorbed.The ensemble which is the most widely used for phase equilibrium calculations is the Gibbs ensemble [8] in which two phases are introduced in separate simulation boxes without an explicit interface.This may be done either by imposing the global volume of the two phases, or the pressure.The temperature and the total number of molecules are also imposed in the Gibbs ensemble.When applied to pure compounds, this ensemble is applied at constant volume to respect the phase rule.Its outputs are then the equilibrium properties: saturated vapour pressure, vaporization enthalpy, and the coexisting densities of the liquid and vapour.In the case of mixtures, the Gibbs ensemble may be applied either at imposed volume or at imposed pressure.In both cases, it provides the equilibrium compositions and densities of coexisting phases.
Besides common average properties like density or energy, the analysis of fluctuations allows one to determine thermodynamic properties like the heat capacity, the The two ways of building a statistical ensemble: Molecular dynamics (MD) and Monte Carlo (MC) simulation.Average properties can be determined from time averages in MD and from ensemble averages in MC.Both averages are equivalent, as a consequence of the ergodicity theorem.
compressibility, the thermal expansion coefficient or the Joule-Thomson coefficient [9, 10].Practically, there are two main methods to simulate statistical ensembles (Fig. 2).The first is through molecular dynamics, which consists in solving the equations of motion, and the second is Monte Carlo simulation, in which a statistical method is used.

Energy of Molecular Systems
The total energy of a molecular system is the sum of its kinetic energy, K, and of its potential energy U: Due to the theorem of equipartition, the kinetic energy is linked with temperature.In molecular dynamics, the solution of the equations of motion allows to compute the kinetic energy at each time, and average kinetic energy may be used to determine the temperature of the system -provided it has reached equilibrium-.Solving the equations of motion makes use of the interatomic forces, or equivalently to the potential energy.In Monte Carlo simulation, the temperature is imposed and the equipartition theorem is used to evaluate the contribution of the kinetic energy to ensemble averages by an appropriate analytical integration.In both cases it is essential to compute the potential energy from the molecular coordinates.It is decomposed classically into intermolecular (or external) energy and intramolecular energy: The intermolecular energy arises from the interaction between different molecules.In the applications considered here, it is decomposed into the following contributions: (5) where U LJ is the Lennard-Jones intermolecular energy arising from dispersion and repulsion interactions, U el is the electrostatic (or coulombic) interaction energy, and U pol is the polarisation energy.The Lennard-Jones energy, which is dominant in low polarity systems like alkanes, is obtained by a summation over all the pairs of force centers: where r ij is the separation distance separating force centers.
Potential energy models considering individual atoms as Lennard-Jones force centers are referred to as "All Atoms" [11][12][13].If the force center represents a group of atoms and is located on the main atom of the group, the model is said "United Atoms" [14][15][16].Alternatively, the force center may be located at an intermediate position between the atoms belonging to the group (Fig. 3), and the model is then known as "Anisotropic United Atoms" [17,18].Most of the applications shown in the present article are based on the latter potential.The electrostatic energy is computed directly from the Coulomb law, assuming that the molecules bear electrostatic point charges: where the sum runs over all possible pairs i, j of partial charges q i , q j , r ij is the separation distance and ε o = 8.8541910 -12 C 2 N -1 m -2 .
The polarization energy originates from the deformation of the electronic cloud of a molecule under the influence of a surrounding polar molecule or solid.We did not introduce this term in the following calculations.This means that the polarization energy is implicitly accounted for in the Lennard-Jones energy or in the electrostatic energy.
The intramolecular energy is assumed to be the sum of the following terms: where U bend is the bending energy arising from the variation of the angle between two successive chemical bonds, and U tors is the torsion energy arising from the variation of the dihedral angle defined by three successive chemical bonds.For an individual bending angle (resp.torsion angle) these potential energies are computed according to the following functional forms: (9) (10) where θ and ϕ are the bending angle and the torsion angle.
U LJ int represents the distant neighbor energy between groups separated by more than three bonds.It is computed with the same Lennard-Jones potential as intermolecular interactions.
In flexible molecules, these contributions are summed over all bending angles, all torsion angles and all pairs of force centers separated by more than three bonds.If a molecule (or a part of molecule) is considered rigid, these contributions are not taken into account.This is the reason why the streching energy, which is linked with the variations of bond lengths, is not accounted for in Equation ( 8).These variations are indeed sufficientily small in the conditions where the following simulations are applied.It must be stressed that alternate expressions can be used instead of Equations ( 9) and (10), which are just empirical correlations.For instance, the bending energy is often written [11,[14][15][16] and the main reason for preferring Equation ( 9) is that its evaluation from atomic coordinates takes less computer time [17,19].

Monte Carlo Methods
The basis of Monte Carlo methods is to generate successive configurations of the simulated system by a statistical method to respect the probability distribution of the desired statistical ensemble.Each elementary move consists in generating a new test configuration by a random change (for instance a random translation of one molecule).Then the energy Principle of the Anisotropic United Atoms (AUA) intermolecular potential.United Atoms refer to the representation of one main atom (here carbon) and bounded hydrogens as a single particle.In classical United Atoms (UA), the force centre is located on the nucleus of the main atom while in Anisotropic United Atoms (AUA), it is located close to the geometrical centre of the constituting atoms.In AUA as in UA models, the Lennard-Jones intramolecular interactions are not counted between sites separated by three bonds or less, because they are accounted in specific terms (bending and torsion).

Anisotropic United atom force center
United atom force center difference ΔU between the test configuration and the last accepted configuration is computed, and the test configuration is either rejected or accepted, based on a specific criterion involving the energy difference.If we consider the example of a translation move in the canonical ensemble, the test configuration is accepted with the following acceptance probability: (11) This means that the new test configuration is always accepted when ΔU is negative, i.e. when its energy is lower than the last accepted configuration.It is then added to the ensemble.If ΔU is positive, a random number is generated between 0 and 1, and the new test configuration is accepted if this random number is lower than exp(-β ΔU).If the random number is larger than exp(-β ΔU), the test configuration is rejected, i.e. the last accepted configuration is duplicated and appears an additional time in the ensemble.
This algorithm, known as Metropolis sampling after the first author of a precursor work [20], provides the basis of most Monte Carlo methods.It is applicable to numerous types of Monte Carlo moves that are necessary to explore as much as possible all system configurations, and to the various possible statistical ensembles.It may be also adapted to bias the generation of the new test configuration, i.e. to generate more frequently favorable configurations that have an increased probability to be accepted.These statistical bias methods are now standard in Monte Carlo simulation to handle either systems with preferred orientation [21], long flexible chains [14,22] or branched flexible molecules [23].Indeed, it is a common problem when complex molecules or dense phases are considered that the acceptance probability of some Monte Carlo moves is very low.
The various types of Monte Carlo moves that are considered in our applications are the following: This variety of moves is needed to sample all the configurations in a given ensemble as well as all possible molecular orientations, locations, and internal conformations.As an example, the moves involved in a typical Gibbs ensemble simulation are indicated in Figure 4.
In addition to these moves, it is also possible to perform test insertions of a molecule in the simulation box.The analysis of the interaction energy of the test particle with the rest of the system allows one to determine the chemical potential of the species [24].This is useful because chemical potential cannot be obtained by simple thermodynamic averages like volume or energy.
A characteristic feature of molecular simulation is that it cannot be used in the close vicinity of critical points, because it would be necessary to simulate systems of excessively large size that are incompatible with reasonable computing times.This is because the characteristic size of density or energy fluctuations increases when the critical point is approached, and larger systems should be considered.In the case of a pure compound, the critical coordinates are obtained by assuming the following scaling law: (12) where β ≈ 0.325 is a universal critical exponent, ρ V and ρ L are the densities of the coexisting vapour and liquid phases, and T c is the critical temperature.When this relationship is complemented by the law of rectilinear diameters: (13) the critical temperature T c and the critical density ρ c can be regressed from Gibbs ensemble calculations performed at lower temperatures.
In the case of binary mixtures, the determination of critical coordinates is somewhat different because we are more interested in the critical pressure at constant temperature.Also, the critical composition must be determined in addition to the critical density and the critical pressure.Nevertheless, expressions comparable to pure compounds can be used, and the critical locus can be also reliably determined [25].
Finally, an interesting Monte Carlo technique to investigate a system at several state points (for instance m different temperatures) is parallel tempering [26].In this technique, m systems are simulated simultaneously.Periodically, they can exchange their temperatures with an appropriate acceptance probability.As a result, one configuration visits successively all the temperatures (equivalently, one temperature is successively represented by distinct configurations).This allows a more efficient sampling of all possible configurations than processing independent simulations.It also takes advantage of the capacity of parallel machines.

Molecular Dynamics
Molecular dynamics consist in solving the equations of motion for every atom in every molecule, accounting for intramolecular and intermolecular interactions inside the system.The Newton equations of motion are expressed in the form: (14) where r i is atom position, i.e. a three dimensional vector (x i , y i , z i ), and F i is the force acting on the atom from the other atoms of the system.The force vector F i can be computed by deriving the potential energy U versus coordinates: (15) In this expression, U includes the interaction energy both with other atoms of the same molecule and with other molecules.The solution of the equations of motion is such that total energy (i.e. the sum of kinetic and potential energy) is constant.Therefore, the solution of the equations of motion at constant volume simulates the microcanonical or NVE ensemble, but not the canonical ensemble, in which total energy fluctuates.It may be shown that sampling the solution of Equation ( 14) at regular time intervals is adequate to sample the NVE ensemble.
Molecular dynamics can be used to obtain equilibrium properties, such as density or enthalpy.However, it is generally less efficient than Monte Carlo simulation.In fact, molecular dynamics is most useful to simulate the dynamic behaviour of a system, particularly its transport properties.By transport properties, we mean here the coefficients which govern the rate of mass transfer (diffusion coefficients), of heat transfer (thermal conductivity) or of momentum transfer (viscosity).For instance, the self-diffusion coefficient of a component may be expressed from the mean squared displacement, according to: Translations, rotations and partial regrowth (internal equilibrium) Volume changes (pressure equilibrium) Principle of the Gibbs ensemble Monte Carlo simulation.Two simulation boxes without an explicit interface are used to represent the phases in equilibrium.Internal Monte Carlo moves (translations, rotations, regrowth , etc.) are used to impose thermal equilibrium at the temperature considered.Volume changes are used to maintain pressure equilibrium, and molecule transfers between phases allow chemical equilibrium between both phases to be imposed.
In a molecular dynamics simulation, it may be checked that the mean squared displacement is indeed a linear function of time, as is the case when Fick's law is obeyed.Thus, the Fickian diffusion coefficient can be rather easily determined from equilibrium molecular dynamics.
Through the Green-Kubo equations, it is also possible to derive shear viscosity and thermal conductivity from appropriate time correlation functions [27][28][29].
In order to perform a molecular dynamics calculation at imposed temperature, some means must be used to change the total energy of the system during the simulation, so that the Boltzmann distribution of energies is respected.This is the purpose of thermostat methods [30,31], in which system energy is changed by altering particle velocities (or center of mass velocities for molecules) at regular intervals.
If a simulation is desired at imposed temperature and pressure, i.e. in the NPT ensemble, a specific technique has to be used to allow for volume variations.The aim of such techniques, known as pressostat or barostat methods [32], is to obtain the right probability distribution (Eq.4).They are equivalent to equilibrating the system with an infinite pressure reservoir.

Software and Hardware Considerations
With current workstations equipped with PC processors, molecular simulation no longer requires expensive supercomputers to be operated on realistic systems.Moreover, the application of molecular simulation techniques will be made easier and easier by the continuous increase of computer memory and speed.Over the last thirty years, the increase of individual processor speed has been observed to follow Moore's law, which states that computer speed doubles every eighteen months at constant cost.Although this trend will probably slow down within the next decade for fundamental reasons, there will be still a substantial reduction of computer time.For difficult problems (such as those involving polymers, for instance), the availability of parallel computers at moderate cost, coupled with adapted techniques like parallel tempering [26], will certainly be a way to generalize the use of molecular simulation in a near future.
Another factor that will favour the development of molecular simulation is the availability of well-tested, multipurpose software.In molecular dynamics, both open source software, such as DLPOLY, and commercial software, are available.Open source Monte Carlo softwares are generally restricted in the range of applications they can address.In contrast, the GIBBS software developed jointly by IFP, CNRS and Université de Paris Sud, which has been used for the applications shown in this article, has a larger coverage since it addresses phase equilibria, phase properties and adsorption in zeolites [25].

Properties of Heavy Hydrocarbons
The first application presented in this article is related to the prediction of thermodynamic properties of hydrocarbons representative of diesel range fractions.Indeed, in this range of carbon atom numbers, i.e. from 10 to 20 carbon atoms, the properties of pure compounds, except for linear alkanes or linear alkylbenzenes, are seldom known.Properties of isoalkanes, cycloalkanes and mixed families are thus poorly documented, although they are present in significant amounts in crude oils and in manufactured products.To face the lack of experimental data for these hydrocarbons, Monte Carlo simulations have been performed in order to obtain useful information for petroleum engineers such as saturation pressures, vaporization enthalpies and liquid densities at different temperatures.For this purpose, two types of Monte Carlo algorithm have been used: the Gibbs ensemble algorithm to predict equilibrium properties at high temperatures, and the NPT algorithm followed by the thermodynamics integration to extend prediction to lower temperatures.Techniques increasing the efficiency of the algorithms such as configuration bias, reservoir bias and parallel tempering were also implemented in the algorithms.All calculations have been performed using the AUA potential optimized on a few reference alkanes [18,33].The main parameters of this potential are given for hydrocarbons in Table 2.It may be observed that these parameters form logical trends.For instance, the Lennard-Jones diameter σ decreases with decreasing number of attached hydrogens in the sequence CH 3 , CH 2 and CH of the alkyl groups.The results obtained on several branched alkanes having from 7 to 19 carbon atoms are discussed in this section.The predicted properties on these hydrocarbons are compared to experimental data obtained from the DIPPR Data Bank [34] when available.Figure 5 shows the saturated vapor pressure curves obtained from simulations.The comparison with available experimental data indicates that good predictions are obtained for the whole range of carbon number studied.It is particularly satisfactory to see that the small difference in the vapor pressure curves of 2,3-dimethylpentane and 2,4dimethylpentane is well predicted, because these two hydrocarbons have the same chemical formula, the same degree of branching and differ only by the position of one methyl group.The comparison between 2-methylhexane and 2,3-dimethylpentane or 2,4-dimethylpentane indicates that simulation is also able to reproduce the slight increase in vapor pressure with increasing branching.The calculations performed on the three C 15 branched alkanes studied are pure predictions since there is no experimental measurement on these compounds.The variation of vapour pressure with the number of branched methyl groups is qualitatively similar to the trend observed for the C 7 branched alkanes, i.e. an increase of the vapour pressure with increasing branching.Such consistency demonstrates the reliability of the calculations even though no experimental data are available.Calculations performed on 2,6,10,14tetramethylpentadecane (pristane molecule) are also found in good agreement with the experimental measurements in the temperature range for which they are available [35].This type of result shows that the intermolecular potential used, optimized on the basis of a short isoalkane only (isobutane), can be used directly to predict equilibrium properties of long branched alkanes.The calculated critical temperature and normal boiling temperature of the several isoalkanes studied are presented in Table 3 along with available experimental data.The results show that the error between experimental data and simulation does not exceed 10 K even near the critical point.

Properties of Natural Gases
Natural gases are mixtures in which gaseous hydrocarbons (mainly methane) are most abundant, but they can contain significant amounts of liquid hydrocarbons, which influence strongly their thermodynamic properties.The advantage of molecular simulation to address these fluids is that it provides excellent predictions of the volumetric behaviour at high pressure, without needing to calibrate specific parameters as it is generally required with cubic equations of state.For  instance, it has been shown that the compressibility factor (Z = PV/RT), the heat capacity, the thermal expansivity, the isothermal compressibility are accurately predicted for light hydrocarbon gases and their mixtures [9].This has been obtained with the same intermolecular potentials as those used for liquid-vapour equilibria, i.e.Möllers potential for methane [49] and the AUA potential [18] for n-alkanes.
If a full analysis is available, as is the case with current techniques of gas chromatography, molecular simulation is expected to show a better predictive capacity than classical equations of state.
An example where molecular simulation was applied to reservoir engineering problems is provided by the production of deep gas condensate fields of the North Sea where pressure may be as high as 110 MPa and temperature may rise to 460 K (HPHT gas).Unexpectedly, the production tests of these fields have shown that the fluid tends to heat up in the production tubing instead of cooling down as it is generally observed.Qualitatively, this is explained by the behaviour of the Joule-Thomson coefficient (i.e. the derivative of temperature versus pressure at constant enthalpy).Although enthalpy is not exactly constant for the gas in real conditions because it exchanges heat with the rocks surrounding the well, the knowledge of the Joule-Thomson coefficient is useful to know the approximate amplitude of the heating effect.The Joule-Thomson coefficient of gases is known to change from positive values (i.e.cooling upon expansion) to negative values (i.e.heating) for pressures of the order of 30-60 MPa, which explains qualitatively the observed behaviour of natural gases.However the Joule-Thomson coefficient is extremely difficult to measure in representative conditions because the temperature effect is very small with the low amount of gas sample available (a few hundreds of cm 3 ).Equation of state calculations are possible [36] but they have not been validated against measured Joule-Thomson coefficients, so that molecular simulation is a good candidate method.
In order to perform the simulation of a real HPHT gas, Lagache et al. [10] obtained its complete composition from different types of chromatographic analysis.The fluid was represented by 19 components, the total number of molecules in the simulated system being 500.This allowed the authors to represent both the distribution of hydrocarbons by molecular weight between methane and C 35 and also the existence of cyclic and aromatic hydrocarbons which were present in the gas (Table 4).
Simulation results are shown in Figure 6, where the computed Joule-Thomson coefficient of the HPHT gas is plotted versus pressure at reservoir temperature (463 K).The trends predicted for methane and ethane, which are the most abundant compounds of the gas, are given for comparison.This prediction could be partly confirmed by a measurement of the inversion pressure, i.e. the pressure corresponding to the change of sign of the Joule-Thomson coefficient (it can be obtained from volumetric measurements).The measured inversion pressure was consistent with the prediction within statistical uncertainties.

TABLE 4
Composition of the HP-HT fluid and representation by a discrete number of molecules [10] Lumped Joule-Thomson coefficient of hydrocarbon fluids versus pressure at 463 K.The lines represent smoothed simulation results for methane (dashed line), ethane (dotted line) and the HP-HT natural fluid (full line), the composition of which is given by Table 4. Experimental Joule-Thomson coefficients are shown for methane (dots).The simulation results on the HP-HT fluid could be checked by the measurement of the Joule-Thomson inversion temperature (see text).

Reinjection of Acid Gases in Deep Reservoirs
In order to improve the economics of producing H 2 S-rich gas fields, one option being considered is to reinject H 2 S in deep reservoir rocks, instead of converting it into sulfur.In these reservoirs, H 2 S may contact oil deposits under high-pressure conditions.The design of such reinjection processes requires numerous data about the phase equilibrium, volumetric and transport properties of H 2 S and H 2 S-hydrocarbons mixtures.Due to the toxicity and corrosivity of H 2 S, experimental data are difficult to obtain, and molecular simulation can be applied to provide predictions of phase diagrams, phase densities [37] and phase viscosities [38].Phase equilibria have been calculated using the Gibbs ensemble Monte Carlo while transport properties have been obtained using equilibrium molecular dynamics in the NPT ensemble.The intermolecular potential of Kristof and Liszi [39] was employed for H 2 S while hydrocarbons were represented by the AUA intermolecular potential previously described.
For three systems, in which the liquid hydrocarbon was a linear alkane (n-pentane, n-hexane and n-pentadecane), Figure 7 shows the pressure-composition phase diagram for a temperature range between 344 and 423 K. Whatever the temperature (above or below the critical temperature of H 2 S), an excellent agreement is found with available literature data [57,58].For the three mixtures, the maximum absolute deviation on vapor phase and liquid phase compositions does not exceed 0.05.The H 2 S/nC 5 mixture has been studied at a temperature below the critical temperature of H 2 S. For this reason, the phase diagram ends up at the vapor pressure of pure H 2 S at the imposed temperature.For the two other mixtures (H 2 S/nC 6 and H 2 S/nC 15 ), that have been studied at 423 K, i.e. a temperature above the critical temperature of H 2 S, no attempt to extend the calculations close to the critical point of the phase diagrams has been performed in the present work.In the case of the H 2 S/n-pentadecane mixture, the solubility of n-pentadecane in the vapor phase is so small that the right side boundary of the two-phase domain is almost confused with the vertical axis corresponding to pure H 2 S.
Similar calculations on binary systems of H 2 S and other hydrocarbon families (naphtenics and aromatic compounds) are now under progress.The generated data will be used in order to develop specific calibrations of the binary interaction parameters needed in engineering equations of state.
Figure 7 (P,x) phase diagrams for three H 2 S/n-alkane mixtures.The simulated liquid phases (filled triangles) and vapor phases (open triangles) are compared to the corresponding experimental liquid compositions (full lines) and vapour compositions (dashed lines) from ref. [57] on the pentane-H 2 S system and from ref. [58] on the hexane-H 2 S and pentadecane-H 2 S systems.

Gas Processing
The development of improved processes for the removal of hydrogen sulphide, carbon dioxide and water from natural gases requires a good knowledge of the thermodynamic properties of the mixtures involved.In addition to the compounds to be removed, these mixtures comprise the hydrocarbons of the gas (mainly methane) and the solvents that are used in the process.Physical solvents like methanol and chemical solvents like ethanolamines are widely used for such separations.However, experimental measurements are very costly and lengthy because of the toxic and corrosive nature of hydrogen sulphide and high pressure conditions.Thus the classical thermodynamic models are often lacking the data that are normally needed to calibrate their interaction parameters.This is especially true when polar compounds like water or methanol are involved, because equations of state do not predict their phase equilibria as reliably as with non-polar hydrocarbons.Through its explicit account of electrostatic interactions between polar molecules, molecular simulation can be of great help.
In Figure 8, we show the results of monophasic simulations of a ternary mixture of methane, hydrogen sulphide and carbon dioxide in the NPT ensemble, that aim to determine the volumetric properties of the gas.With the use of the Kong combining rules, it has been shown that volumetric properties could be predicted with an excellent accuracy compared to state-of-the-art equations of state [37].The purpose of these simulations is to check the behaviour of the equation of state which will be used repeatedly in process design on comparable mixtures.In order to perform the test on a wide range of conditions, 72 state points were investigated.Thanks to the parallel tempering technique, nine simulations (one per pressure level) were sufficient to investigate all state points.For each pressure level, seven temperatures were run in parallel, attempting to exchange configurations between neighbouring temperatures every 5000 steps.As it can be inferred from the figure, the statistical uncertainty on density is very low (approximately 0.2%).This example shows that Monte Carlo simulation can now produce rapidly a fair amount of accurate volumetric data on mixtures.However, care must be taken to test whether the system is monophasic or diphasic.Indeed, a monophasic simulation may compute the density in metastable conditions.For instance, if the true equilibrium corresponds to a biphasic equilibrium, the NPT simulation may compute the density of a metastable homogeneous liquid or of a metastable homogeneous gas phase, depending on the initialization.In simulations, the fluid does not have the time to undergo phase separation inside the simulation box, in the same way as a bottle of sparkling water does not undergo phase separation immediately after its opening although this brings the liquid into a two-phase region.
Phase diagram calculations in the Gibbs ensemble are exemplified by the case of the H 2 S-methanol system at the temperature of 348 K, for which data are available [40].As shown in Figure 9, the agreement of simulation results and experimental data is excellent.This has been obtained by using an intermolecular potential from the literature for methanol [41] and the same potential as above for hydrogen sulphide.Similarly, it has been shown that the phase diagram of the propane-methanol system, which presents an Phase diagram of the H 2 S-methanol mixture at the temperature of 348.15K from Gibbs ensemble simulations (symbols) and experimental measurements (lines) [40].
Figure 10 Size distribution of methanol clusters in the liquid phase of a propane-methanol mixture at 2 MPa and 343 K.The straight line, which corresponds to the exponential formula shown, is a least squares fit of simulation results in this semilogarithmic plot.azeotrope, is also very well predicted [25].Encouraging results have been also obtained on the water-methanol system by us and by other investigators [42].
When a polar fluid like hydrogen sulphide is mixed with a very polar fluid like water, it might be thought that the system is rather miscible at high pressure.In fact, this is not the case [43] as phase split occurs between an aqueous phase with 3-4% of dissolved H 2 S and a H 2 S-rich liquid phase with approximately 3% water.As shown in a previous work [37], the phase bahaviour is well predicted by simulation with ordinary intermolecular potentials for water [44] and H 2 S [39].An interesting outcome of simulation is that it helps to understand fluid structure.For instance, the snapshots of the H 2 S-rich liquid at high pressure showed that water molecules are not dispersed, but tend to form small clusters.This phenomenon contributes to increase the solubility of water in this phase compared with non-associating fluids like hydrocarbons, where such clusters are not observed.Moreover, the distribution of cluster sizes can be obtained in some cases, as illustrated by the propane-methanol system at 343 K and 2 MPa, shown in Figure 10.The linear trend in this semi-logarithmic plot indicates an exponentially decreasing distribution, as is expected from classical theories of associating fluids [45] in which the association constant (K i = x n+1 /x n ) is related to the free energy of the association reaction A n + A → A n+1 .The main difference with simulation is that the association parameters are not fitted to explain phase equilibrium, but they are obtained as a result of the intermolecular interactions.The departure of simulation results from the classical behaviour (i.e. from the straight line) for small cluster size is due to specific effects which are ignored in classical association models.

Solubility of High Pressure Gases in Polymers
The extension of gas solubility calculations, as shown in Section 2.3, from alkanes to very long chain alkanes, such as polyethylene, is of considerable industrial interest.Indeed, polymers are used in oil production to ensure the internal and external pipe leakproofness.Knowledge of the factors influencing permeation of small molecules, like carbon dioxide, hydrogen sulfide or water, through these materials is required to evaluate blistering and corrosion risks.
In the present work, the solubility of carbon dioxide in molten polyethylene at 433 K has been investigated by means of Monte Carlo simulations.These calculations have been performed in the Gibbs ensemble at imposed pressure, using reptation moves to relax the polymer matrix and rotational bias [21] for the gas molecules.The optimized AUA potential [18] has been used for polyethylene and an All Atoms potential model [46] for carbon dioxide.The polyethylene matrix was modeled by 15 linear chains of 100 carbon atoms.Slow convergence, associated with large fluctuations, has been encountered during these simulations, especially at high pressure.At 15 MPa for instance, more than 300 million Monte Carlo iterations were required and two independant runs have been performed to check the validity of the results.The calculated solubilities, expressed as the mass of carbon dioxide dissolved in 100 g of empty polymer, are presented in Figure 11.The predicted values obtained with this model system correspond fairly well with experimental measurements of Sato et al. [47] and of Chaudhary and Johns [48] on polyethylene.Another information that is very important to characterize is the volumetric behavior of the system.The swelling of the polymer, i.e. the relative difference in volume between the saturated polymer and the empty polymer in the same conditions of pressure and temperature, has also been calculated from simulation data (Fig. 11).For instance, the predicted swelling is ∼8% at 10 MPa and ∼15% at 15 MPa.These significant values show that there is a substantial rearrangement of the polymer as it dissolves CO 2 .

APPLICATIONS IN REFINING
Processing of heavy crude oils is increasingly used to provide liquid fuels (gasoline, diesel) for which consumption is much higher than heavy fractions.Among the various elementary processes, hydrotreating is very important as it allows cracking of heavy compounds without generation of excessive amounts of coke, thanks to the use of high pressure hydrogen.Another important role of hydrotreating processes is to remove sulphur-bearing compounds such as thiols and thiophene derivatives.Sulphur is thus converted in H 2 S, which is distributed between the liquid and vapour phases.In order to design such units, it is essential to use well-tested phase equilibrium models.This design is traditionally based on correlations calibrated on available measurements of H 2 Salkanes phase equilibria.However, the lack of experimental data at typical process temperature (430°C) obliges one to use these models far above the investigated range of temperatures.There is little probability that such data will be available soon, for the simple reason that hydrocarbons are highly unstable at such high temperatures.Molecular simulation is therefore very useful to provide the safest possible extrapolation from known phase equilibria (below 200°C) to process conditions.As an example of the results that can be obtained, we simulated the phase behaviour of the H 2 S-nC 20 system at 400°C, using the same potentials as the reservoir engineering applications discussed above.As shown in Figure 12, the original parametrization of the Grayson-Streed model shows systematic deviations, and the more recent parametrization by Erbar shows a wrong curvature of the bubble point curve in the (P,x) diagram.In these systems, molecular simulation has allowed development of an improved correlation that accounts consistently for the influence of pressure, temperature and carbon number on equilibrium compositions.Process design is now benefiting from simulation-based improvements without requiring significant changes in process engineering software.
High temperature solubility of H 2 S in eicosane at 400°C.Molecular simulation results are used to calibrate an improved correlation following the Grayson-Streed scheme.
Molecular simulation is also very useful to understand and predict the adsorption properties of zeolites, which are often used for separations in hydrocarbon processing.The reader is referred to a recently published book in which these applications are exposed [25].

CONCLUSION
From a fundamental standpoint, molecular simulation presents numerous advantages versus classical thermodynamic models, in compensation for the important computing time it requires.Among these advantages, we can quote the ability to account for the detailed molecular structure-which is important to determine property differences between isomers-, the capacity of decomposing the energy according to repulsion, dispersion and electrostatic contributions, or the explicit consideration of fluid structure when association is significant.Although molecular simulation still requires empirically determined parameters to evaluate intermolecular interactions, it shows unprecedented predictive capacity to investigate new compounds, new properties, new mixtures or new temperature and pressure conditions.
These capacities of molecular simulation have been illustrated in the present article by a few application examples in the oil and gas industry.All branches of the oil and gas industry can benefit now from the predictions of pure component properties by molecular simulation.This is particularly interesting for many families of heavy hydrocarbons that are poorly known (e.g.isoalkanes, napthenoaromatics or alkylpolyaromatics).In reservoir engineering, molecular simulation is used to predict the properties of natural gases at high pressure (such as the Joule-Thomson coefficient, which is particularly difficult to measure by experimental means) or the behaviour of hydrogen sulphide-hydrocarbon systems (for which the few data available are not sufficient for a reliable design of acid gas reinjection processes in deep reservoirs).In oil and gas transport and processing, molecular simulation is used to predict the behaviour of hydrocarbon mixtures with water and polar solvents like methanol or amines.It is also a useful tool to model the solubility of gases in pipe-coating polymers.Finally, molecular simulation is very useful in hydrocarbon processing to predict phase equilibria in process conditions (e.g.around 400°C), where experiments would be impossible due to the fast degradation of hydrocarbons.
Is molecular simulation operational?Although we must recognize that the above achievements are pioneering applications, there are several indications supporting a positive answer.Firstly, molecular simulation no longer requires expensive supercomputers, and it can be operated efficiently on workstations or Linux clusters that are not more expensive per processor than ordinary PCs.More importantly, multipurpose software is available.For instance, all the examples shown can be processed with GIBBS, a Monte Carlo software developed jointly by IFP, CNRS and Université de Paris Sud which is commercially available.Finally, learning to use this technique for efficient applications is made easier by the availability of reference textbooks oriented either toward fundamentals [2, 3] or toward applications [25].Additional fluid properties which are expected to be predicted by these methods in the forthcoming years are interfacial properties (such as interfacial tension), chemical equilibria [42,50], and viscosity, among others.Therefore, it can be conjectured that molecular simulation will soon play a significant role in applied thermodynamics, with an intermediate position between experimental measurements and existing thermodynamic models.It will be most useful in cases when few data are available, being then cheaper and faster than experiments and more reliable than classical models.

Figure 1
Figure 1 Example of a simulation box of a mixture of cyclohexane and H 2 S in the liquid phase.Cyclohexane molecules are represented by six anisotropic united atoms and H 2 S molecules by a single united atom.
1. translation of a molecule; 2. rotation of a molecule around a randomly selected axis; 3. volume change of the simulation box (NPT ensemble, Gibbs ensemble); 4. transfer of a molecule from one box to the other (Gibbs ensemble); 5. destruction of an existing molecule or insertion of a new one (Grand Canonical ensemble); 6. regrowth of a part of a molecule (flexible chains); 7. flip, i.e. rotation of an atom around the axis formed by its immediate neighbours (flexible molecules); 8. reptation, i.e. destruction of one end of the chain and growth of the other end (flexible chains); 9. displacement, i.e. destruction of an existing molecule and insertion of a new one in the same box; 10. pivot, i.e. rotation of the part of a molecule located beyond some atom, around this atom (flexible chains).

Figure 11 Concentration
Figure 11 Concentration of carbon dioxide dissolved in polyethylene at 433 K. Simulation results (filled diamonds) are compared to experimental data (open diamonds).Associated polymer swelling calculated from simulation results is also shown (filled triangles).

TABLE 3
Estimated critical temperature and normal boiling temperature of the branched alkanes considered.