From Organic Geochemistry to Statistical Thermodynamics: the Development of Simulation Methods

This work is a synthesis of researches to which the author contributed successively in three areas of applied research. The first of these areas is petroleum organic geochemistry, which is aimed to understand the formation of petroleum reservoirs from detailed analyses of the composition of fossil fuels, experimental simulations in the laboratory by pyrolysis, and reconstructions of the geological history of the sedimentary basins concerned. These reconstructions, which imply the mathematical simulation of many coupled mechanisms such as thermal cracking, fluid flow in porous media, sedimentation, sedimentary rock ;compaction and inter-phase equilibria, were the subject of the author's main contributions in this field. The results obtained at IFP in the 1980s, as part of a multidisciplinary team project, provided the basis for numerous industrial applications in connection with the new discipline of basin models. The second of these areas is petroleum thermodynamics, which involves measuring and predicting the equilibrium and transport properties in crude oil and natural gas production and processing. The contribution of the work presented in this vast area was to characterize petroleum fluids at high pressure, ;both by experimental identification of the phase diagrams specific to certain HP-HT fluids, and the prediction of their properties by means of equations of state. It also allowed a more accurate understanding of how the phase equilibria between oil and gas were likely to alter the composition of petroleum fluids. The third area of study is statistical thermodynamics, whose applications have recently soared with the generalization of molecular simulation. Monte Carlo and molecular dynamics methods not only offer an understanding at microscopic scale, but also a more reliable prediction of the equilibrium and transport properties in many cases where classical macroscopic models are limited by their empirical character. By relying on this occasion on close collaboration with a university team, it was possible to implement improved algorithms and to develop a new intermolecular potential to meet the needs of the petroleum industry. Thus effective applications of molecular simulation materialized in the processing of H2S gases, separations by adsorption in zeolites, and the calculation of properties of commercially unavailable hydrocarbons. These three areas, which appear very different at first glance, display significant similarities, like their dominant articulation around physical chemistry, their need of new software, and the empirical optimization of parameters whereof a physical meaning is nonetheless expected.


INTRODUCTION
Fossil fuels-oil, coal and natural gas-represent the principal energy source for humanity today.The rate at which mankind exploits these resources, which have taken millions of years to come into being, may seem irrational, but their massive replacement by other energies is not seriously conceivable in the short term.The energy supply of the planet will therefore depend for a long time to come of the renewal of petroleum resources by the discovery and production of new reservoirs.This renewal is contingent on the continuous improvement of the techniques employed.In exploration, the days are past when the petroleum geologist was content to prospect in the field to pick or position a well, and he relies today on an in-depth knowledge of the mechanisms of formation of the oil reservoirs.In production, no longer are we content to let the oil gush from the production wells, and waterflooding and gas injection in the reservoirs is now systematic to optimize recovery.In refining, the distillation of crude oil is no longer sufficient, and has for decades been supplemented by processing methods and catalytic cracking to supply more efficient and less pollutant fuels.
This article, which summarizes the research projects to which the author contributed from 1980 to 2002 at Institut français du pétrole and the University of Paris Sud, offers an opportunity to dwell on some of these areas of applied research.
The first of these areas is petroleum geochemistry, which underwent a scientific revolution like all of earth science, after the advent of plate tectonics in the 1960s.This discipline, which participated in the tremendous impetus aimed at the time to quantify the behavior of the earth's crust in physical and chemical terms, had succeeded in demonstrating the biogenic origin of oil.The presence of many "geochemical fossils", hydrocarbons directly derived from molecules known in living beings, served to demonstrate that oil and natural gas were primarily produced by a slow thermal degradation of the sedimentary organic matter (Durand, 1980;Tissot and Welte, 1984).
In this field, the 1980s saw the development at IFP of ambitious numerical simulation methods to understand the formation of oil reservoirs.This work was aimed to improve the kinetic models of thermal cracking proposed some years earlier (Tissot and Espitalié, 1975), and also to propose an original approach to describe the transport processes responsible for its accumulation in reservoir rocks.These achievements were the focus of a multidisciplinary project integrating geosciences, chemistry, thermodynamics, mechanics and applied mathematics.
The second area of this recapitulation is the thermodynamics of petroleum fluids, and chiefly the experimental characterization of fluids at high pressure, which in petroleum production technology corresponds to the interval between 50 and 120 MPa.These investigations, propelled by the wish to produce very deep reservoirs, helped to identify phase diagrams of a specific type when a natural gas contains heavy n-alkanes which crystallize at low temperature.They also culminated with the measurement of volumetric properties and viscosities in the same conditions.As to modeling, we will above all have to address the difficulty of predicting the volumetric properties and phase equilibria of crude oils and natural gases.This inadequacy is not only explained by the empirical character of routine equations of state like Redlich-Kwong-Soave (Soave, 1972) and Peng-Robinson (Peng and Robinson, 1976), but also by the improper consideration of the fluid compositions: above 7 or 10 carbon atoms, engineers generally ignore the measured composition (Pedersen et al., 1989).To improve the predictions of the properties of gas condensate, it was proposed to represent the heavy fraction by distributions of saturated and aromatic hydrocarbons.As we shall show, the results encouraged further work in this direction, which was unfortunately limited by the lack of experimental data on heavy hydrocarbons other than n-alkanes.This type of finding is also not specific to reservoir fluids: despite the progress offered by excess free energy models (Huron and Vidal, 1979), by contribution methods (Fredenslund et al., 1977), by the application of the principle of corresponding states (Lee and Kesler, 1975) and by equations of state based on concepts of statistical thermodynamics (Chien et al., 1983;Chapman et al., 1989), thermodynamic models were often inaccurate.In many cases, the acquisition of experimental data enabling their adjustment was lengthy and costly: heavy hydrocarbons not commercially available, highly toxic components, very high pressures, new molecules, etc.In the specific example of multicomponent adsorption equilibria relevant to separation processes, classical thermodynamic models had to rely on selectivity measurements, which encumbered the search for better adsorbents (Ruthven, 1984).
These observations logically encouraged a trend towards a third area of research from 1997, molecular simulation.By taking account of the molecular structure and the different forms of intermolecular energy, simulation can offer more reliable predictions in extended areas.The continued progress of computation capacity enabled this discipline to occupy a leading place in physical chemistry (Allen and Tildesley, 1987;Frenkel and Smit, 1996).Thanks to the Gibbs ensemble method (Panagiotopoulos, 1987) it is now possible to calculate phase equilibria, while molecular dynamics provides access to transport properties.In terms of grand canonical simulation, adsorption processes can be simulated at a pertinent scale in microporous materials used to perform certain separations that are difficult in chemical engineering (Lachet et al., 1997).For the effective use of molecular simulation in applied thermodynamics, many laboratories focus on algorithmic improvements to Monte Carlo methods employed to calculate the equilibrium properties (Smit et al., 1995) or methods of molecular dynamics (Rowley, 1994;Dysthe et al., 1999a;Dysthe et al., 1999b;Dysthe et al., 2000).In cooperation with the University of Paris-Sud, IFP attempted to contribute thereto by proposing an original method to calculate the bubble point of multicomponent mixtures and to calculate derivative properties like heat capacity or Joule-Thomson's coefficient.However, the accuracy of molecular simulation calculations always depends critically on the intermolecular potentials whose adjustment is relatively empirical.Continuing a project initiated by Toxvaerd (Toxvaerd, 1990;Toxvaerd, 1997), it was proposed to refine the determination of the anisotropic united atoms potential parameters for alkanes, on the basis of the thermodynamic properties which must be represented in chemical engineering.This work helped set the stage for a more general potential for other nonpolar hydrocarbons, and also for polar organic compounds, by also relying on the results of quantum chemistry.

Kinetics of Oil Formation
The pyrolysis of sedimentary organic matter with temperature programming under an inert gas stream (Espitalié et al., 1985) helps to experimentally simulate the formation of oil and to simply identify the kinetic control of this thermal cracking process.If pyrolysis is rapid (50°C/min) the maximum rate of liberation of hydrocarbon substances is observed around 500°C.If it takes place more slowly, at 0.3°C/min, this maximum release occurs around 400°C (Fig. 1).Based on this type of observation, it is possible to develop a kinetic model involving several parallel first order reactions obeying the Arrhenius' law, where the main activation energies range from 50 to 60 kJ/mole (Ungerer et al., 1985;Burnham et al., 1987;Ungerer and Pelet, 1987;Ungerer, 1990).The determination of the parameters of the model (the weight of each reaction and their common preexponential factor) relies on an optimization method by a least squares criterion.
It is no doubt difficult to provide a solid justification for this model, proposed by Tissot (Tissot and Espitalié, 1975), rather than another.For example, the approach of Burnham at the Lawrence Livermore Laboratory (Burnham and Braun, 1990) led independently to slightly different parameters by using a continuous spectrum of parallel reactions rather than discrete reactions.On the other hand, there is no doubt that control of the parameters by experiments of very different durations is an excellent key to ensure that the problem of optimizing the twenty free parameters of the model is properly posed, in other words, that the solution is stable.It is only from the time that this stability is achieved that the model can be extrapolated to very slow evolutions like those taking place in the subsurface.
In sedimentary basins, where the burial rate is a few tens to a few thousands of meters per million years and where the temperature rises by 20 to 40°C per kilometer of depth, the rate of temperature rise is about 1 to 100°C per million years.The kinetic models fitted to the pyrolysis results predict that the maximum oil formation must occur at around 120-150°C in these conditions.It is in fact in sediments which have reached these temperatures that manifestations of thermal cracking are systematically observed, such as higher content of volatile hydrocarbons and a decrease in the hydrogen content of the residual solid organic matter-since the hydrocarbons formed are richer in hydrogen than the initial solid.The model also correctly predicts the shift of the pyrolysis peak towards the high temperatures which is used to characterize the thermal alteration of organic matter (Ungerer et al., 1985).The validation of this approach required a patient analysis of natural samples not only involving IFP (Ungerer et al., 1985;Ungerer and Pelet, 1987;Forbes et al., 1991), but also the Lawrence Livermore Laboratory (Sweeney et al., 1987) and the KFA team in Jülich (Welte et al., 1997), among other leading laboratories active in this field.One of these validation examples is shown in Figure 2.
Admittedly, it is difficult to reconstruct to better than +/-10°C the temperature changes undergone by the rocks in their natural evolution, thereby limiting the accuracy of this validation.Despite these uncertainties, kinetic models describing the transformation of organic matter into oil are now intensively employed by the petroleum industry to delimit the zones of the sedimentary basins in which the oil has been formed.
In the industry, another major challenge is to predict the composition of the hydrocarbons formed, because it is much less advantageous economically to discover a gas reservoir than an oil reservoir, the second term designating a liquid form of petroleum in the reservoir conditions.To meet this Pyrolysis curves at different heating rates of a sample of sedimentary organic matter taken from a source rock in the Paris basin, and histogram of activation energies of the kinetic model determined from these curves (Ungerer and Pelet, 1987).
Figure 2 Trend as a function of depth of the transformation rate of organic matter into hydrocarbons, estimated from pyrolysis experiments on samples of wells of the Norwegian margin (rhombus) and calculated using a kinetic model like the one shown in Figure 1 (Forbes et al., 1991).demand, a first answer is to couple the inert gas flushing pyrolysis system with a cold trap system.Depending on the trap temperature, the curve of liberation of methane, hydrocarbons from C 1 to C 5 , and C 1 to C 15 can be recorded separately.Based on this additional experimental information, it was possible to make the kinetic model discussed above compositional, in other words, to supply at each stage the stoichiometry of the products formed in four main classes: methane, C 2 -C 5 , C 6 -C 15 and C 15+ (Espitalié et al., 1988).The studies conducted on organic matter of different origins (Fig. 3) helped to reveal systematic differences: coals of plant origin thus produce higher proportions of light hydrocarbons than matter of marine of lacustrine origin (Espitalié et al., 1988;Welte et al., 1997).However, the general predominance of the C 15+ fraction in the reaction products clearly indicated that pyrolysis under a vector gas stream represented an open environment in which the hydrocarbons generated were preserved from the secondary cracking reactions that yield light hydrocarbons.
Although the compositional kinetic model thus developed was improved, it was unable to predict by itself the quantities of gas found in numerous basins.In consequence, the modeling of the secondary reactions in the conditions of the sedimentary basins (pressures of several hundred bar) proved indispensable.

Kinetics of Secondary Degradation Reactions
At the time when the development of basin models demanded better consideration of secondary reactions, it turned out that the pyrolysis of oils under pressure was also the subject of detailed investigations for the needs of enhanced recovery by in situ combustion (Behar et al., 1988).From the theoretical standpoint, the excessively large number of elementary reactions involved made a detailed radical kinetic model inconceivable, and the kinetic models used in refining did not apply to complex systems such as crude oils.Furthermore, the few studies available attested that the influence of pressure was to significantly accelerate the cracking reactions up to 100 bar, and that the influence of Histogram of activation energies used to model the composition of products generated by primary cracking by means of a parallel reaction kinetic model for two types of organic matter (Espitalié et al., 1988).
pressure appeared to lessen beyond this limit.An empirical kinetic model was therefore proposed in which each of the unstable fractions (hydrocarbons other than methane and heavy products such as resins and asphaltenes) could be degraded by a first order reaction forming both the lighter and heavier fractions.In this system, the only final terms are methane and a polyaromatic carbonaceous residue called coke.This ultimate term also corresponds to the result of the thermodynamic stability analysis (Dayhoff et al., 1967).
The optimization of the parameters of the secondary cracking model is difficult because they are numerous (activation energies, stoichiometric coefficients, preexponential factor).By employing a constrained optimization method, we can impose that the stoichiometric coefficients respect the conservation of carbon and hydrogen by assuming a constant elementary composition for a given fraction.But beyond the choice of the optimization method, it is crucial to conduct experiments covering a broad range of heating times and temperatures, without which certain parameters would be undetermined.In particular, it is important that the same degree of advancement of secondary cracking be simulated both by short term pyrolyses at high temperature and by longer pyrolyses at low temperature.Based on these principles, which demanded fruitful mutual feedback between models and experiments, it was possible to develop kinetic models for two typical examples of oils (saturated and aromatic) based on five fractions: the four fractions used in primary cracking (C 1 , C 2 -C 5 , C 6 -C 15 and C 15+ ) and coke (Ungerer et al., 1988).As for the primary cracking reactions, the extrapolation of the model to geological evolutions yielded results in good agreement with the highest temperatures at which the oils had been found in the sedimentary basins, or around 180°C.Subsequently, this kinetic model was improved on several occasions by increasing the number of fractions (Behar et al., 1991) but by retaining the same general approach.
The coupling of primary cracking and secondary cracking models proved to accurately reflect the composition of the hydrocarbons generated by the main types of organic matter by pyrolysis in a closed medium, to the extent of the reproducibility of the experimental results (Ungerer, 1990).However, the application of this coupled model to sedimentary basins clashed with a major problem of principle, namely that neither pyrolysis in open medium, nor pyrolysis in closed medium, can be expected to fully simulate the thermal cracking mechanisms in a petroleum basin: the open medium would simulate the case of a system in which the hydrocarbons could travel to a colder reservoir as soon as they are formed, while the closed medium would simulate a basin in which these hydrocarbons cannot leave the place where they have been generated and would undergo the same thermal evolution as the solid organic matter.To simulate natural cases, which are intermediate, it is therefore necessary to consider the coupling between petroleum migration-name given to all its movements in the subsurface-and cracking mechanisms.By modeling the migration in a simplified manner, it was possible to compile a balance of the hydrocarbons expelled from the source rock in the drainage area of a known reservoir in the North Sea, and to confirm that their composition was compatible with that of the reservoir (Forbes et al., 1991).This verification would ideally have demanded a dynamic flow model like the one discussed below, but it was not advanced enough to consider this.
Before concluding the discussion of kinetic models, it is important to state what is needed to have a more reliable model.As demonstrated by the work of Doué and Guiochon (Doué and Guiochon, 1968), the thermal cracking of alkanes at high pressure is characterized, inter alia, by the virtual absence of olefins in the reaction products, and is roughly explained by a radical mechanism.To examine how these results could be extended to cracking processes involving organic geochemistry, a project conducted alongside the ENSIC of Nancy attempted to completely simulate the cracking of light hydrocarbons (hexane and hexane-benzene mixture) under pressure (Kressmann et al., 1990;Kressmann, 1991).Figure 4 shows the kinetic scheme obtained.
This work mainly offered the advantage of demonstrating the importance of bimolecular processes between the hexane molecules and the radicals (Fig. 4), so that the apparent Figure 4 Radical cracking of hexane at 350°C (Kressmann et al., 1990) activation energy of the cracking process (50-60 kcal/mole) is much lower than that of monomolecular scission (over 80 kcal/mole).The importance of bimolecular reactions also questions the strong hypothesis of secondary cracking models, namely that the degradation rate of a fraction is independent of the concentrations of the other species.This being said, the radical mechanisms involved already include several hundred elementary processes for very simple systems with a low degree of advancement.Their resolution demanded solving systems of ordinary nonlinear and poorly conditioned differential equations, which demanded several hours of computer time on the computers available at the time (Cray II).The extension of the radical approach to real systems requires both the implementation of lumping methods in order to reduce the number of equations, as well as reliance on theoretical chemistry methods to calculate the kinetic parameters of the many elementary processes involved.This research path has advanced considerably since then (Duin et al., 2001), but its application to complex systems is still an objective of present research paths.

Modeling of Oil Migration
The migration of oil and gas is a very general phenomenon since hydrocarbon accumulations are generally found in sedimentary rocks different from those in which oil is formed: the former, called reservoir rocks, are permeable rocks (sandstone, porous carbonates), while the latter, called source rocks, are generally very fine grained rocks (shale, marl, coal).Even in basins which are known to have generated large quantities of oil, an understanding of migration processes is essential to locate the costly exploration boreholes on the most favorable geological structures.
It was recognized very early that migration in permeable rocks was explained by a gravitational segregation and that this could be described by the classical formalism of flows in porous media.Similarly, it was acknowledged in the 1970s that the entrapment of hydrocarbons in the reservoirs was explained, depending on each case, by the impermeability of the overerlying rock-salt, anhydrite-or by capillary forces at the interface with the reservoir rock (Berg, 1975).On the contrary, the mechanisms at play in the source rocks were still the subject of debate in the early 1980s, because capillary forces were considered to pose an obstacle to migration in these rocks, whose pores were submicronic in size: some researchers felt that the hydrocarbons migrated in solution in the water expelled by the sediments during their compaction (Price, 1976), while others considered that molecular diffusion was the dominant mechanism (Leythaeuser and Yükler, 1982), and yet others suggested a two-phase flow in which the hydrocarbons formed a separate phase from the water (McAuliffe, 1980).
Following a critical analysis of the literature, it emerged that the main mechanism responsible for the expulsion of oil from the source rocks was the last of the three mentioned above, namely migration in a separate phase (Durand, 1983;Ungerer et al., 1983).The underlying reasons can be summarized by simply stating that the solubility of liquid hydrocarbons in an aqueous phase is much lower and that molecular diffusion is a mechanism that would have excessively favored light hydrocarbons over the heavy fractions.On the contrary, expulsion in a separate phase better reflected the high electrical resistivity found in the source rocks, and a certain similarity of composition between the crude oils and the products remaining in the source rocks.
Once migration in a separate phase in source rocks such as reservoirs was accepted, it remained to find an adequate formalism to approach its modeling.This was a particularly ambitious task, because at the time, the only team that had tried to model the migration of oil at the scale of a basin had limited consideration to its movement in the reservoir rocks (Welte and Yükler, 1981).The approach pursued (Ungerer et al., 1984) was characterized by consideration of the following processes: -the geometric evolution of the basin was represented by a mobile mesh in which the rows of meshes are added successively to simulate the deposition process of new sediments; -the thermal evolution of a level is inferred from the evolution of its depth, based on the geothermal gradient observed or from the solution of the Fourier's heat transfer equation; -the formation of oil in the source rocks is described by the kinetic model of primary cracking discussed above; -the fluid flow is described by a two-phase Darcy's law in which the permeability varies by several orders of magnitude according to the particle size distribution, porosity, hydrocarbon saturation (this term designating the fraction of pore volume occupied by the phase concerned) and the prevailing stress; -compaction, i.e. the irreversible reduction of porosity of the sedimentary rocks, is represented by a law linking the porosity and the effective stress.
As complex as this formalism may appear, it was impossible to discard a single one of these elementary processes to represent the formation of oil reservoirs realistically.On the contrary, the model tended to sideline a vitally important mechanism with the phase equilibria between oil and gas, whose important implications during migration were the subject at the same time of the thesis of E. Nogaret at LIMHP (Nogaret, 1983).
To concretize this mathematical model by a numerical simulation, it was illusory to make a mere vertical unidimensional study as done by the authors dealing with the kinetics of oil formation.In fact, the formation of many reservoirs cannot be explained without migration with a strong horizontal component.This is why the route of twodimensional numerical simulation (vertical + horizontal) was forthwith selected, despite the numerical problems associated with the resolution of the system of equations.Given the total duration to be simulated, which ran into millions of years, it was in fact necessary to apply time steps of several thousand years to remain within reasonable computer time.This resulted in a numerical instability which demanded the development of implicit numerical schemes.Besides, the use of deformable and nonrectangular meshes was necessary to adapt to the limits of the geological discontinuities (Fig. 5).This was an additional difficulty, because this type of mesh had barely been addressed in the literature.
Despite the shortcomings of the early simulations, the formalism proposed for migration rapidly proved capable of simulating all the mechanisms governing the formation of oil reservoirs in various cases: expulsion from source rocks, segregation by density in reservoir rocks, entrapment by and leakage through caprocks.Simultaneously, the interbedding of permeable and impermeable levels resulted in a very strong permeability anisotropy at the basin scale, which reconciled the requirement of a low vertical permeability-to explain the entrapment of the oil-and a high horizontal permeability, to supply reservoirs located far from the zones where the cracking was active.Beyond the petroleum aspects, the model accurately described the pressure and fluid flow regime at the scale of the basins.For example, the introduction of low permeabilities for shales (about one nanodarcy, or 10 -21 m 2 ) helped to explain that despite the very large timescale of the mechanisms simulated, the pressure profiles often displayed very wide deviations from equilibrium at great depth (Fig. 6).
Once these preliminary results were obtained, many improvements were made.A first alternative was to consider sedimentation mechanisms in greater detail and explicitly simulate heat transfer mechanisms, whose transient character is significant in certain basins (Doligez et al., 1987).A second alternative concerned numerical resolution: test of a finite element method in connection with the thesis of V. Bouvier (Bouvier, 1989), implementation of a compositional scheme in that of I. Faille (Faille, 1992).At the same time, this coupled model of fluid migration and heat transfer saw its first application to real exploration cases, albeit for synthesis rather than prediction (Ungerer et al., 1987) since it was limited to two dimensions (Fig. 7).
The basic principles of the basin models set forth in the 1980's by the IFP team have not been fundamentally questioned, either in the field of thermal cracking or of migration.However, many improvements continued at the scientific level in teh 1990s: three dimensional modeling, explicit consideration of oil-gas phase equilibria during migration, representation of faults.At the industrial level, the use of basin models was progressively extended and generalized in oil exploration, and the IFP group is still one of the leading software suppliers in this area.Variation in fluid pressure as a function of depth observed in two wells of the Norwegian margin (q) and calculated using the twodimensional migration model (continuous line).The hydrostatic equilibrium profile is shown as a reference (Ungerer et al., 1990).Geological section across a section of the Norwegian margin indicating the fluid migration directions calculated using a two-dimensional migration model.The hydrocarbon saturation is the fraction of pore volume occupied by the hydrocarbon phase (Ungerer et al., 1990).

Phase Diagram of Natural Gases Containing Heavy Hydrocarbons
The natural gases discovered by Elf in the deepest boreholes of the North Sea and the Aquitaine basin in the late 1980s sometimes displayed unusual properties compared with the fluids normally encountered: particularly high dew point pressures, sometimes over 50 MPa, crystallization of paraffins at low temperature.These fluids also displayed a characteristic composition, rich in methane and containing heavy hydrocarbons up to C 30 or more.To start producing these reservoirs, it was indispensable to have a clear understanding of the phase diagram of these fluids and to relate it to the composition.This raised specific experimental difficulties due to the high pressures and temperatures (190°C and 110 MPa).In cooperation with Elf, the phase diagrams of natural fluids (Fig. 8) and synthetic fluids (Fig. 9) were determined by means of high pressure cells (Ungerer et al., 1995;Ungerer et al., 1998).
The composition of the synthetic fluids had been selected in order to approximate that of natural fluids, methane being the majority compound (70% or more) with the remainder consisting of light liquid hydrocarbons.Their study revealed that the addition of heavy aromatic or paraffinic hydrocarbons in a small molar fraction (1% or less) resulted in the appearance of (gas + solid) and (gas + solid = liquid) domains at temperatures close to ambient.It also showed that the heavy hydrocarbons had the effect of considerably increasing the saturation pressures.These observations were not really surprising, because they fitted into the logic of the behavior of methane-heavy hydrocarbon binary mixtures investigated experimentally at the University of Delft (van der Kooi, 1981;Glaser, Peters et al., 1985;Flöter, 1999), at IFP (Arnaud et al., 1996;LeRoy et al., 1997) and at the CNRS (Marteau et al., 1998).For example, we found the negative slope of the curve of appearance of a solid phase at high pressure in the (P, T), diagram, as it had been observed on methane rich binary mixtures (Fig. 8).
Nevertheless, it was impossible to observe a diagram comprising a three-phase region in the P, T plane for a binary system due to the phase rule, whereas it was possible for our synthetic fluids since they had three components or more.Our experiments therefore provided a conceptual guide for the interpretation of the phase diagrams of natural fluids.These measurements were then used in the literature to validate thermodynamic models describing thermodynamic models describing crystallization in complex systems (Nichita et al., 1999).In a different perpspective, the simultaneous use of experiments on model fluids and equations of state proved a very efficient approach to predict the distribution of light liquid hydrocarbons in oil and gas phases at depth.It was thus possible to provide organic geochemists with well-tested tools to understand subtle composition variations in oil reservoirs (Carpentier et al., 1996).Figure 9 Phase diagram of a synthetic fluid consisting of methane (79%), pentane (20%) and α-methylnaphthalene (1%), showing the same type of behavior as gas condensate in the previous Figure 8 ( Ungerer et al., 1998).

Measurement and Correlation of Volumetric Properties at High Pressure
The prediction of volumetric properties by cubic equations of state is obviously insufficient, and to correct this defect, Péneloux et al. (1982) proposed the volume translation technique.The great advantage of this technique is to avoid modifying the representation of the liquid-vapor equilibria on which the cubic equations of state are fitted.Based on the literature data, it is possible to show that volume translation allows systematic errors to subsist (de Hemptinne and Ungerer, 1995).In particular, it was found that for n-alkanes from pentane to decane, the error committed at high pressure on the molar volume of a pure substance increases linearly with temperature.This finding led to the recommendation of a variable volume translation with temperature, whose parameters are adjusted to liquid volumes at 50 MPa (Ungerer and Batut, 1997;Sant'Ana et al., 1999).Figure 10 gives an example of the improvement made.This approach helps avoid the inconsistencies observed at high pressure when volume translation is adjusted to accurately describe liquid volumes at saturation.It also offers the advantage of being compatible with the mixing rules classically used during the lumping of several pure substances into pseudocomponents.
In the case of a mixture, an adequate volume translation is not sufficient to clearly represent the molar volume, because the excess volume at constant pressure must also be represented.This quantity is generally small in the case of mixtures of liquids, for which it is less than 1% of the volume (Katzenski and Schneider, 1982;Ashcroft et al., 1992).In the case of contrasted mixtures (gas-liquid) at high pressure, the measurements taken at IFP by J.F. Arnaud during his thesis project showed that the scale of the excess volume is significantly greater, raising the problem of its modeling.Volume translation has no influence on the excess volume in a cubic equation of state, but it is possible to introduce a quadratic mixing rule on the covolume to faithfully represent this quantity (Arnaud et al., 1996).Measurements of excess volume thereby provide independent information of the phase equilibria measurements, the former being rather influenced by the covolume mixing rule and the latter by the attractive term of the equation of state.

Representation of the Phase Behavior of Gas Condensates
A condensate gas is defined by the fact that a decrease in pressure at the reservoir temperature leads to the condensation of a liquid phase (Pedersen et al., 1989).The composition of these gases is dominated by methane, of which the molar fraction is 60 to 80%, and by light hydrocarbon gases.The liquid hydrocarbons of these fluids are mainly hydrocarbons with C 5 to C 10 , but may include heavier hydrocarbons: up to C 20 or even C 30 .As shown by the study of binary systems and synthetic mixtures discussed earlier, the phase behavior of methane rich mixtures is very sensitive to the composition of the heavy hydrocarbons.Hence it is not surprising that the prediction of the phase equilibria of gas condensate by equations of state is very sensitive to the representation of the heavy fraction of the gas.Unfortunately, the precise composition is generally unavailable, since chromatography cannot provide a molecular analysis above C 10 due to the many isomers present (Durand et al., 1989).It is hence routine practice to optimize the thermodynamic properties of the heavy fraction to best represent the liquid deposit curve.The resulting model is difficult to extrapolate to other temperatures.
To improve this procedure, distribution functions can be used to represent the hydrocarbons of the heavy fraction: this is the track explored by M. Sportisse in her thesis (Sportisse, 1996;Sportisse et al., 1997).In this work, continuous distributions were employed to reduce the number of parameters to be optimized and not to reduce the number of equations of the phase equilibrium calculation.This approach helped to test several types of functions (exponential, gamma) to describe the heavy fraction by distributions of alkanes or aromatic hydrocarbons.It thus appeared that phase equilibria are just as sensitive to the chemical family as to the analytical form of the distribution.For example, the exclusive use of a distribution of n-alkanes leads to a systematic underestimation of the dew point pressure, and it Curves of liquid deposition of a gas condensate observed ( ) and calculated by representing the fraction by distribution functions (Sportisse et al., 1997).The distribution of heavy hydrocarbons is determined here in order to represent the behavior observed at 104.4°C, and the curves corresponding to other temperatures are predicted without adjustment.
was necessary to include distributions of alkylbenzenes and polyaromatic hydrocarbons to correct this effect.For gases of which the detailed composition is known, the heavy fraction thus modeled proves to be fairly representative of the real composition.It also allows an improved prediction of the density of the condensate and a better extrapolation of the curves of liquid deposition towards the low temperatures (Fig. 11).

What Prospects for Petroleum Thermodynamics?
The petroleum industry would like to benefit from more predictive methods for a broad range of applications, including liquid-vapor equilibria, the crystallization of solid phases (paraffins, gas hydrates), the flocculation of asphaltenes, and adsorption, etc.However, progress is certainly difficult to achieve in this area, and the few results set forth in the previous chapter may appear modest when seen against the efforts made.Thermodynamicists are sometimes tempted to blame the lack of experimental data for a good share of these difficulties.It cannot be denied that in the field of crude oil and natural gas, the rarity of measured equilibrium properties of pure substances is a serious hindrance for hydrocarbons with more than ten carbon atoms other than n-alkanes, and this is no doubt one factor that limits potential progress.It is also true that the equilibrium data of binary mixtures at high pressure are not substantial enough, in the same way as crystallization measures under pressure, flocculation under pressure and adsorption selectivity on well-known adsorbents.The systematic acquisition of these data would be very expensive, and the industry is not prepared to finance this effort for profits that seem too far in the distant future.
Yet if the analysis is taken further, this situation could be explained as much by the inadequacy of present theories as by industry decisions.
Thus in the example of equations of state, it is questionable whether really decisive progress has been achieved since the work of Peng-Robinson and Soave.While many equations of state have improved the description of a particular system, up to now this has always occurred to the detriment of the generality of the equation.An eloquent example in this respect is the Lee-Kesler's method of corresponding states (Lee and Kesler, 1975), very accurate for components similar to the reference substances, but unrealistic for substances that are too different (Hemptinne and Ungerer, 1995).The SAFT equation also helps to clearly represent the equilibria of certain systems (Chapman et al., 1989;Gil-Villegas et al., 1997), but by remaining in the general scheme of linear molecules and, for each system, relying on a specific parameterization of the distribution function.Moreover, like any analytical equation of state, it does not reproduce the scaling laws near the critical point.On other equations of state, the addition of "crossover" terms helps to satisfy these scaling laws (Albright et al., 1986;Rainwater and Williamson, 1986), but this procedure requires knowledge of the critical point, which one would in fact have preferred to calculate from the equation of state in the case of mixtures!Finally, another prohibitive weakness of existing theories is their difficulty in predicting differences in properties between near isomers (for example, between isoalkanes with the same number of branches).The consequence of these inadequacies is a very strong demand for experimental measurements to adjust model parameters to specific cases.
A more severe criticism could be leveled concerning the prediction of adsorption selectivities which the industry needs to design adsorption separation processes on crystalline microporous solids.The use of thermodynamic models similar to those employed in liquids or the reliance on statistical models (Ruthven, 1984;Tournier, 2000) does not help predict selectivities, but only to correlate them on the basis of experimental measurements.
Faced with these findings, it is likely that the community of thermodynamicists will have to completely change the problematics to provide general answers.Without much chance of being wrong, we can claim that this will entail reliance on theories wherein the molecular structure and the nature of the molecular interactions will be described in greater detail.Given the insufficient results of the models inspired by statistical thermodynamics, one can even claim that each molecule will have to be represented individually, as well as any changes in conformation, at least in delicate cases (adsorption, differences in properties between isomers).This is why the future of thermodynamic models is rather to be sought in terms of molecular simulation techniques, the subject of the last part of this article.

Computer Algorithms Based on Monte Carlo Methods
The first numerical simulations using Monte Carlo methods were nearly contemporary with the first appearance of computers, dating half a century back (Metropolis et al., 1953).The first statistical bias algorithms (Rosenbluth and Rosenbluth, 1955) appeared soon thereafter, followed by the Widom insertion test which helps to calculate the chemical potential (Widom, 1963), mentioning only the best known of the pioneering work at the time.In the 1960-70s, we therefore had the theoretical concepts necessary to address complex systems, as attested by the still very real utility of the reference works written at the time (Münster, 1969;McQuarrie, 1976).However, the applications remained limited to small size molecules for many years.
Only in the 1980s did the continuous increase in computation capacities significantly broaden the field of possibilities of Monte Carlo methods.This is no doubt why many essential algorithms for calculating thermodynamic properties were developed quite late.Among them is the calculation of phase equilibria in the Gibbs ensemble (Panagiotopoulos, 1987) , the simulation of flexible alkanes by a configurational statistical bias (de Pablo et al., 1992;Smit et al., 1995), the simulation of water by an orientational bias (Cracknell et al., 1990), the calculation of saturation vapor pressures at low temperature by the integration of the Clapeyron's equation (Kofke, 1993), the simulation of polymer materials by means of Monte Carlo movements acting on a limited segment of the molecule (Dodd et al., 1993) and the calculation of phase equilibria not corresponding to a statistical ensemble stricto sensu thanks to pseudoensembles (Mehta and Kofke, 1994).These achievements allowed the development of the Gibbs code as part of a collaborative project between IFP and the Physical Chemistry Laboratory of the University of Paris Sud (Mackie et al., 1997).Despite these breakthroughs, Monte Carlo methods cannot today meet all the needs of applied thermodynamics.This is largely due to the lack of appropriate intermolecular potentials, as we shall discuss below.But it is also due to the lack of algorithms in two directions: the calculation of phase equilibria under particular constraints and statistical bias methods adapted to the calculation of complex molecules.

Calculation of Phase Equilibria Under Particular Constraints
The Gibbs ensemble only partly covers the needs connected with phase equilibria, because the quantities imposed are the temperature, the total number of molecules, and a third quantity which is the volume or pressure.In view of the applications, algorithms are also needed in which the quantities are imposed.When classical thermodynamic models are used, it is for example routine to impose the composition of the liquid phase for a bubble point calculation or that of the vapor phase for a dew point calculation.
Generally speaking, it is important to be able to plot in the (P, T) plane the limits of the domain wherein a complex mixture of a given composition becomes two-phase, i.e. its phase envelope (Michelsen, 1980).Among the other types of specifications are calculations at imposed enthalpy which serve to calculate the variation of temperature or pressure resulting from an isenthalpic expansion or from the input of a given quantity of energy.The direct calculation of azeotropes and critical points is another important need.
In this area, an original Monte Carlo algorithm has been proposed for the direct calculation of the bubble point (Ungerer et al., 1999;Ungerer et al., 2001) using two simulation boxes without interface but nonetheless coupled via changes in volume (Fig. 12).This algorithm can be summarized by stating that it consists in calculating the probability of acceptance of a molecule transfer from one phase to the other as if we worked with the Gibbs' ensemble, but by refraining from updating the configuration of the liquid phase, whose composition thus remains rigorously constant.The vapor phase, in which the updates are effectively made, undergoes insertions and destructions comparable to those characterizing the grand canonical ensemble.Its mean composition and mean pressure supply the desired results in calculating the bubble point.
Compared with the simulation techniques previously available (Vrabec and Fischer, 1995;Escobedo, 1998;Miyano, 1998), this work helped to avoid dependence on an iterative process to determine the liquid phase chemical potentials and hence offered undeniable progress.However, it appeared that the direct algorithm for calculating the 284 Widom test Algorithm of direct calculation of the bubble point at imposed temperature.The formulas shown for acceptance rates correspond to a configurational bias algorithm (Ungerer et al., 2001).
Figure 13 Liquid-vapor equilibrium constants in bubble point conditions calculated using the proposed algorithm (Ungerer et al., 2001).
bubble point had been proposed independently at the University of Wisconsin by Escobedo under the term of faked move (Escobedo, 1999).This author having also proposed an alternative method adapted to the calculation of dew pointsfor which the direct algorithm is unstable-it was decided not to further pursue the development of this type of algorithm.In fact, Escobedo in the meantime proposed an algorithm for isenthalpic calculations (Escobedo and Chen, 2001) and the azeotrope calculation algorithms progressed (Pandit and Kofke, 1999).It would have been difficult to stand up to this competition.The sidelining of this research subject appears a posteriori to be a wise decision, especially since the determination of intermolecular potentials demanded major efforts.However, this work provided many lessons.Inter alia, it served to identify the approximate nature of the factorization of the partition function as classically employed.Strictly speaking, the partial integration of the partition function reveals a factor whose dependence as a function of the conformation is presumed to be negligible in Monte Carlo techniques (Ungerer et al., 2001).While this hypothesis seems justified for molecules with strong intramolecular potentials (Go and Sheraga, 1976) , it is not the case of long chain polymers, whose rigidity is weak at the scale of the chain.
The practical limitations of the Gibbs ensemble to a relatively restricted temperature interval (from the normal boiling point to close to the critical temperature) is another reason to employ complementary methods.For molecules of substantial size (C 10 to C 30 ) the temperatures at which thermodynamic properties are to be calculated are in fact often lower than the boiling point.In the case of pure substances, thermodynamic integration has proved to be an ideal method for calculating the saturation vapor pressures at the same time as the densities and vaporization enthalpies (Ungerer et al., 2000).To accelerate this type of calculation in the future, it is necessary to employ parallel tempering techniques (Yan and de Pablo, 1999) which, by exchanging configurations between simulations made at different temperatures, offer the advantage of effectively using the capacity of parallel machines, without necessarily demanding extreme speed of communication between processors.However, in using these methods, one must be sure to allow the configurations to decorrelate sufficiently before the exchanges.

Improved Statistical Bias
Thanks to configurational statistical bias techniques, it is possible to calculate the phase equilibria of linear or branched alkanes, which would be otherwise untractable even at the cost of very long computer time.However, the possibilities of simulation are extremely limited concerning large size molecules.Thus in the example of squalane-a C 30 isoalkane-the use of configurational statistical bias in the Gibbs ensemble did not permit access to sufficiently low temperatures for comparison with experimental measurements (Neubauer et al., 1999).To extend these calculations to temperatures closer to process requirements, it is undoubtedly necessary to improve the statistical bias algorithms.
In the case of long chain linear alkanes (C 20 to C 35 ), E. Bourasseau, in his DEA diploma course at IFP, implanted a local rotational move that serves to relax the central portion of the chains, which is untouched by growth movements.This algorithm, in which a carbon atom rotates about the axis formed by its two neighbors, allows single-phase simulations at temperatures far below the boiling point (Bourasseau et al., 2002).The advantage of this local move is that it applies to any molecule possessing a linear chain, regardless of its structure (alkylbenzenes, ethers, cycloalkanes).For very long chain alkanes (polyethylene), the introduction of a reptation move will be useful for improving effectiveness, but this move cannot be extended to all molecules.
In the case of isoalkanes with many branches, the local rotational movement is not sufficient to explore all the conformations of the middle chain because the branches all represent fixed points.To correctly simulate these compounds at low temperature, complementary Monte Carlo movements are necessary.One possibility is the pivot of a moiety of the molecule about one the bonds of the main chain, or a concerted rotational movement of several consecutive centers (Dodd et al., 1993), or the simultaneous placement of two atoms at branches (Macedonia and Maginn, 1999).Cyclic alkanes represent another family of molecules for which recent research has encountered difficulty in calculating the phase equilibria while taking account of the different possible internal conformations (Neubauer et al., 1999).In his first thesis year, E. Bourasseau developed a particular transfer movement using a canonical reservoir of conformations, periodically updated by local rotational movements (Bourasseau et al., 2002).To obtain a significant proportion of accepted movements, the transfer is made in two steps.In the first step, designed to find cavities of substantial size in the liquid, several insertion tests are performed on a Lennard-Jones' particle at places selected at random in the receiver phase, and one of them is selected for its Boltzmann's probability.In the second step, many insertion tests of the complete molecule are performed at the place selected in random orientations, the internal conformation being sorted at random in the reservoir.This statistical bias algorithm helped to reduce computer time by a factor of three per accepted transfer.
With a view to calculating larger and less flexible molecules, it is important to further improve the effectiveness of the transfers.How to do this?To start with, analyze the source of the effectiveness of the bias discussed above with respect to an unbiased transfer movement.Its success mainly stems from the fact that a cheap selection step (the test of several positions by a Lennard-Jones potential) helps to limit to a favorable position the energy calculation of the second step, which is more expensive because involving all the force centers of the molecule.In this way, the complete computation is avoided in sites where the overlap with the other molecules would in any case prevent the transfer from being accepted.In the case of inflexible molecules displaying a significant elongation, like polycyclic hydrocarbons (decaline, naphthalene), an interesting research direction would be to design a bias in three steps instead of two: the first step would be unchanged (insertion tests of Lennard-Jones' particles at several places), the second step would comprise insertion tests of a second Lennard-Jones particle at a fixed distance from the position selected for the first particle, and the third step would consist in testing several conformations and orientations of the molecule about the axis of the positions selected.In this movement, the aim of the first two steps would be to find elongated cavities in the liquid, so as to reserve the costly calculations of the final step for the most favorable sites.
A final point worth examining concerns the calculation of the chemical potential by insertion tests.Although the problem has been the subject of much work in the literature, no simple method appears to provide a definitive solution in the case of dense liquids (Kofke and Cummings, 1997).Reliance on progressive insertions is certainly possible, but perhaps difficult to insert in the normal course of a simulation, where other qualities than the chemical potential need to be calculated.A recent research direction consists in coupling a destruction test and an insertion test (Boulougouris et al., 1999).This approach was successfully applied to calculate the solubility of alkanes in water, which is a difficult test (Boulougouris et al., 2001).It could conceivably provide a fairly general solution to the direct calculation of the chemical potential at low temperature, when the Gibbs ensemble is inapplicable.In the case of pure substances, this could offer an alternative to thermodynamic integration.

Why New Intermolecular Potentials?
Why is the development of intermolecular potentials still an important problem, despite the abundant reporting in the literature?There are many reasons for this, and they may be worth recalling.As a consequence of the present inability to accurately calculate the dispersion energy by ab initio calculation methods, the intermolecular potentials are of empirical form and their parameterization can assume various forms according to the final objective.Many potentials so far developed (Jorgensen and Madura, 1984;Cornell et al., 1995;Sun, 1998) were not intended to represent the saturation vapor pressure, but only the density of the liquid phase and the enthalpy of vaporization in ambient conditions.Logically, these potentials offer a poor prediction of the saturation vapor pressure (Chen et al., 1998).In fact, this quantity is essential for predicting phase equilibria because it guarantees the correctness of the chemical potential in a liquid phase for pure substances.
If we consider investigations aimed to represent liquidvapor equilibria, the "all atoms" potentials in which each atom is explicitly represented (Chen and Siepmann, 1999) require very long computer time.The united atom potentials, in which the hydrogen atoms are not individually represented, offer the advantage of reducing computer time by about one order of magnitude, while allowing a good plot of the coexistence curve of pure substances.However, the united atom potentials stricto sensu which place the force center on the nucleus of the main atom of the group (usually carbon) are somewhat lacking in physical meaning and generality.This is demonstrated by the fact that in these potentials (Smit et al., 1995;Martin and Siepmann, 1998;Nath et al., 1998) the Lennard-Jones diameters attributed to the CH 3 , CH 2 and CH groups are higher than that of methane (Möller et al., 1992) and that in the case of iso-alkanes, the energy parameter attributed to the CH 3 group varies from 50 to 80 K depending on the position of the group and the size of the molecule (Martin and Siepmann, 1999;Nath and de Pablo, 2000).On the contrary, the anisotropic united atom potential proposed by Toxvaerd (Toxvaerd, 1990;Toxvaerd, 1997) seems to offer greater physical meaning, and this why the Physical Chemistry Laboratory decided to investigate this route from the start of its cooperation with IFP.

Figure 14
Calculation by simulation of the saturated vapor pressure of n-alkanes: comparison of TRAPPE united atom intermolecular potential (Martin and Siepmann, 1998), the AUA3 potential (Toxvaerd, 1997) and the optimized AUA4 potential (Ungerer et al., 2000).
It was therefore proposed to optimize the parameters of the Toxvaerd's potential for n-alkanes, attempting to represent several properties in a wide range of temperatures: saturated vapor pressure, enthalpy of vaporization and liquid density (Fig. 14 ).
Since then, E. Bourasseau has demonstrated that this potential correctly represented the properties of long chain alkanes (Bourasseau et al., 2002b) and M. Lagache obtained good predictions of thermodynamic quantities derived for ethane and butane (Lagache et al., 2001).This performance should be credited to the physical meaning of the parameters obtained by optimization.In particular, the fact that optimization can locate the positions of the force centers near the geometric center of the atoms of the group is highly encouraging.As to iso-alkanes, it has been demonstrated that the AUA potential optimized on isobutane could reasonably predict the difference in properties between n-heptane and its branched isomers.

Families of Components Considered
Since this work, E. Bourasseau continued the development of the AUA potential to adapt it to other hydrocarbon families.In the case of cycloalkanes, his work allowed the simultaneous representation of the properties of cyclopentane and cyclohexane with parameters of the CH 2 group very close to those of n-alkanes (Fig. 15).The good representation of cyclooctane, which is not part of the molecules to which the potential was fitted, shows that this potential is more Results obtained by E. Bourasseau (Bourasseau, et al., 2002a) using an AUA type intermolecular potential for cycloalkanes.The potential is fitted to the properties of cyclopentane and cyclohexane, but not to those of cyclooctane, which is also well represented.
general for cyclic hydrocarbons than the united atom potentials stricto sensu (Errington and Panagiotopoulos, 1999;Neubauer et al., 1999).In the case of olefins, the influence of the branches and the unsaturation position is clearly indicated by the AUA potential recently optimized by E. Bourasseau.Compared with the united atom potentials proposed by Spyriouni et al. (1999) or Wick et al. (2000) for olefins, this potential serves to represent more equilibrium properties (density, saturation pressures and enthalpies of vaporization) by means of parameters consistent with those employed for saturated hydrocarbons.A precise and transferable AUA potential has also been developed for aromatics, objects of close cooperation with A. Mackie's group of the University of Tarragona.This research direction therefore offers highly promising prospects for mixed families encountered in petroleum fluids: alkylcycloalkanes, alkylaromatics, naphthenoaromatics, etc.Another interesting prospect is the development of potentials for organic molecules possessing polar functions, as initiated by J. Delhommelle successfully for organic sulfides and thiols (Delhommelle et al., 2000).This approach consisted in relying on the results of quantum chemistry to determine the partial charges while preserving the anisotropic united atom potential of the alkanes for the nonpolar portion of the molecule.In cooperation with the TotalFinaElf group, S. Kranias carried out a work of the same type to represent other polar functions (ketones, aldehydes), in the course of publication.M. Lagache is currently working on a comparable development for the highly toxic organomercury compounds present in traces in some natural gases.Her work presents an additional difficulty insofar as the experimental data are less abundant and where the molecular structure of these compounds is poorly known.Hence there is every reason to believe that the family of potentials of anisotropic united atoms will be significantly enriched in the coming years.

Form of Interaction Potential and Mixing Rule
By thus proceeding into the development of intermolecular potentials, it is worth raising the question of the form of the interaction potential and the mixing rule.One could, for example, consider the replacement of the Lennard-Jones's potential by a Buckingham's potential, in which the exponential form of the repulsion is known to better represent the behavior of noble gases (Wu and Sadus, 2000).It is nonetheless possible that the larger number of parameters to be determined of the Buckingham's potential would mean indeterminations in the optimization process, and it is uncertain whether this potential will provide the same improvement for anisotropic united atoms as the one procured for monatomic gases.
In the same range of ideas, one could question whether the choice of the Lorentz-Berthelot's mixing rule is the best, insofar as the work done by J. Delhommelle on noble gases (Delhommelle and Millié, 2001) or by other researchers on more complex mixtures (Potoff et al., 1999;Errington et al., 2000) tend to show that the Kong and Waldmann-Hagler's rules could yield better results, at least for interactions between groups of very different size.These rules rely in particular on a geometric mean on the product εσ 6 and not on ε, which conforms more to the theoretical form of the dispersion energy.
The tests performed on methane-alkane mixtures in connection with the Master degree work of E. Soldi-Lose in 2000 failed to identify any tangible effect of the mixing rule, so that it did not seem advisable to question the present choice so far.However, the recent study by A. Wender of H 2 S -CO 2 mixtures showed that for temperatures slightly higher than the critical temperature of the mixture, its volumetric properties are particularly sensitive to the choice of the mixing rule.On this system, the results speak clearly in favor of the Kong and Waldmann-Hagler's rules, and could even help to decide between the two.An interesting research direction is hence to extend this type of comparison to other systems to determine the degree of generality of these conclusions.

Methodology of Optimization of Intermolecular Potentials
Another point worth examining is the optimization method developed during the work on n-alkanes, on which sufficient data is now available to consider improvements.Undeniably, its strong point is the definition of the composite error function which integrates three quantities (saturation vapor pressure, enthalpy of vaporization and liquid density) at very different temperatures, on several pure substances.Provided allowance is made for necessary adjustments (for example, the consideration of volumetric data measured in supercritical conditions), this type of error function proves to be a powerful tool because it helps to use numerous independent data while demanding a moderate number of simulations.One could possibly consider entering derived quantities like those discussed in the following chapter in the error criteria.However, the slowness of convergence of these quantities obtained by fluctuation methods and the lesser availability of the corresponding experimental data, tend to imply that this would not be a major improvement to the optimization process.
On the contrary, the gradient method employed to determine the minimum of the error function is sometimes a weak point.Besides, the gradient method as proposed may display wide instabilities when the initial parameters are too far from the optimum.If this is the case, we cannot hope to find the optimum directly, and it must be accepted to decrease the error function step by step in the gradient direction.Moreover, the calculation of the gradient component by finite differences is a lengthy, inaccurate and difficult process.To avoid it, one can directly calculate the partial derivatives of a mean quantity X with respect to any parameter λ of the potential from the known fluctuation formula (Gray and Gubbins, 1984) where U is the total potential energy.
In exchange for the calculation of a few additional terms during the simulation, this direct calculation of the gradient components is likely to considerably accelerate the empirical portion of the development of intermolecular potentials.The first tests of this method conducted recently by E. Bourasseau effectively yielded absolutely positive results for representing the properties of olefins (Bourasseau et al., 2003).

The Use of Anisotropic United Atoms for Adsorption
In terms of applications, adsorption in zeolites is a field in which the industry would be happy to have a united atom potential that displays transferability qualities comparable to those observed in liquid-vapor.On the scientific level, it is also interesting to examine whether the displacement of the force centers with respect to the carbon nuclei provides additional predictive power for simulating adsorption in zeolites bearing the compensation cationic charges.It is one of the major objectives assigned to the thesis P. Pascual, begun in autumn 2001 under the direction of A. Boutin on the adsorption of alkanes, to develop this new approach.
Based on the approach which yielded good results in liquid-vapor, one could propose to develop the zeolitealkanes potential on the basis of adsorption isotherms covering a wide range of numbers of carbons and temperatures.In a first step, the consideration of zeolites devoid of extra-structural cations like silicalite would help to parameterize the oxygen-alkane and T-atom-alkane interactions.In the second step, the study of zeolites bearing cations in known crystallographic positions (like the NaY, KY and BaX faujasites) would help to derive the dispersionrepulsion parameters of cation-alkane interactions.And finally, if this succeeds, the approach could be extended to the adsorption of aromatics, olefins and cycloalkanes, so as to obtain a coherent interaction potential for adsorption and fluid phases.It is likely that this process will be confronted with a lack of data-if only of adsorption isotherms of pure substances measured at different temperatures-and it is therefore logical to propose an experimental program designed to fill these gaps.

Towards the Industrial Use of Monte Carlos Methods
A first concern of the IFP and CNRS investigations was to assess the accuracy with which it was possible to describe the properties of hydrocarbon gases at high pressure.The work of B. Neubauer had already established the feasibility for calculating one of the properties of natural gases containing several components, with highly encouraging results (Neubauer et al., 1999).It was therefore tempting to grapple with a more difficult problem, namely the prediction of derivative thermodynamic quantities: heat capacity, thermal expansion coefficient, compressibility factor, Joule Thomson's coefficient.This topic is the first part of the thesis of M. Lagache, carried out at the LCP in connection with a Concerted Research Project involving six other French university laboratories and oil companies.To calculate the derivative quantities more directly, M. Lagache applied a fluctuation technique discussed briefly by Jorgensen (Jorgensen, 1986;Lagache et al., 2001).It is in fact possible to calculate the derivative of a mean quantity X with respect to an intensive variable Y by considering the covariance of X and the extensive quantity Z associated with Y: where λ is a factor depending on the variables considered.
Since Monte Carlo methods do not calculate the kinetic energy, it is necessary to use them only to calculate residual properties (deviation from ideal gas properties) and to calculate the ideal contributions from experimental measurements.Once this somewhat delicate division has been made, the results obtained are excellent, better in fact than the initial likelihoods.Without altering the intermolecular potentials which were nonetheless fitted to the first order derivatives of the free energy, it is in fact possible nearly quantitatively to represent the variations in quantities that are especially difficult to predict like the heat capacity at constant volume (Fig. 16).Note that this test is particularly difficult since the C v (ρ) curve at constant temperature shows a peak which the equations of state do not reproduce, even qualitatively.The correct representation of the peak of C v by simulation offers hope that the slight divergence of this quantity at the critical point can also be reproduced if we use sufficiently large simulation boxes.
Another field of application of delicate phase equilibrium calculation is the behavior of acid gases, i.e. natural gases containing significant contents of H 2 S and CO 2 .In general, phase equilibria involving H 2 S have been the subject of fewer experimental measurements as those exclusively involving hydrocarbons, due to the toxicity and corrosive nature of this gas.It is therefore extremely interesting to employ molecular simulation to predict these equilibria, as initiated by J. Delhommelle during his thesis work.To supplement these investigations, it was decided to simulate ternary phase diagrams in the water-methane-H 2 S system, for which few data are available.In this study, it was important first to validate the choice of the intermolecular potentials by comparison with the properties of the pure substances of the mixture.It thus turned out that with the potentials considered for water and H 2 S (Delhommelle, 2000), the properties of the pure substances were less well described when the polarization energy was taken into account.While this was not surprising in the case of water, since the TIP4P potential had been optimized without taking account of polarization et al., 1983), this finding was difficult to explain in the case of the polarizable model of H 2 S proposed by Delhommelle.Finally, nonpolarizable potentials were employed for H 2 S (Kristof and J. Liszi, 1997) and for water (TIP4P).Once this stage of the choice of the potential was terminated, the analysis of the binary water-H 2 S system (Fig. 17) provided the opportunity to employ a liquid-liquidvapor equilibrium calculation at imposed total volume to calculate the three-phase equilibrium pressure.The prediction of the phase diagram is very good, and some details are extremely well rendered.The simulation helped to explain why the solubility of water is higher in the H 2 S rich liquid phase than in the vapor phase: this behavior is probably connected with the self-combination of water molecules to form oligomers (dimers, trimers, etc.) in this phase (Fig. 18).This mechanism, rather unexpected for solubility levels in the range of 2-3%, shows that we are not in the regime of infinite dilution for water.
Now that simulation has demonstrated its possibilities for acid gases, we can consider applying it more systematically for reinjection projects of these gases into the subsoil, after having separated them from the hydrocarbons: prediction of volumetric properties at high pressure, determination of the heat of mixing, phase equilibria of H 2 S-hydrocarbon mixtures.The training student work of A. Wender recently completed at IFP helped to confirm these hopes by rapidly supplying a large number of high quality predictions.By resorting to molecular dynamics, it should be also possible to determine transport properties such as viscosity and thermal conductivity in the future.
Another type of application of Monte Carlo methods consists in predicting the properties of hydrocarbon components with more than ten carbon atoms, i.e. those making up most of diesels and kerosenes.It is in fact surprising to find that in this carbon range, only a few families, like n-alkanes, alkylbenzenes and condensed polyaromatics (naphthalene, phenanthrene, anthracene) have been the subject of sufficient measurements to accurately determine their equilibrium properties.In the absence of pure substances available commercially at reasonable prices, measurements are difficult to take on many isomers of little known families such as isoalkanes and cycloalkanes, with more than ten carbon atoms.The same applies to mixed compounds like alkylcyloalkanes and naphthenoaromatics, which are also well represented in crude oils (Petrov, 1984;Tissot and Welte, 1984).Since molecular simulation reasonably describes the differences in properties between isomers with the potentials recently developed and possesses good properties of extrapolation to high carbon numbers, it could supply more reliable predictions than the group contribution methods used today.This type of prediction proves feasible (Fig. 19) and is now seriously considered for better use of the detailed analytical results on petroleum products obtained by analytical methods coupling chromatography with mass spectrometry.
In the same range of ideas, the prediction of the behavior of polymers is a field in which molecular simulation can make a large contribution.For example, the solubility of gases at high pressure in amorphous polymer matrices is an important point that must be controlled for the design and Figure 17 Phase diagram of the water-H 2 S as a binary system at 343 K. operation of flexible lines which used these materials considerably.The convergence of this type of system demands reliance on specific algorithms such as reptation, recoil growth (Consta et al., 1999) and gradual insertion (Kaminski, 1994) to effectively sample all the conformations.It is to be hoped that the use of these algorithms will help to address industrial problems successfully.
A final industrial challenge, and certainly not the least, concerns the prediction of adsorption in zeolites.This prediction is essential to design adsorption separation processes used to separate isomers whose liquid-vapor equilibrium properties do not permit the economic use of conventional separation by distillation.Thus meta is separated from paraxylene, the latter being an important intermediate in petrochemicals.It is also by adsorption that iso-alkanes are separated from n-alkanes, the latter being undesirable in premium gasolines due to their poor octane number, and must therefore be recycled to the reactor, where they undergo conversion to iso-alkanes.These adsorption processes, as well as many other industrial processes using the same principle, are based on the specificity of moleculeadsorbent interaction at molecular scale.In this respect, it is not surprising that classic thermodynamic models fail to predict the adsorption selectivities in the absence of selectivity measurements on the absorbent concerned, even if these models incorporate concepts of statistical mechanics.This implies a real opportunity for the development of molecular simulation methods capable of supplying preliminary indications before experimental measurements become available.Prediction of the saturation vapor pressure of n-pentadecane and of two C 15 alkanes by means of the AUA4 potential.The three simulations corresponding to the highest pressures were obtained directly by two-phase simulation (Gibbs ensemble), while the lower pressures were obtained by thermodynamic integration using single-phase simulations.
Compared with cases of fluid phase equilibria, the simulation of adsorption equilibria is a more difficult problem in many respects.We are still far from having a family of potentials transferable from one zeolite to another and from one adsorbate to another: so far, the potentials used are in fact specific to a particular system.In his first thesis year at LCP, P. Pascual succeeded recently in producing a first improved anisotropic united atom potential.This potential displays transferability properties much better than those of classic united atom models (Vlugt et al., 1998).Beyond this result, it must be kept in mind that the frequent presence of cations in zeolites makes electrostatic interactions and polarization more important, and the possibility that these cations could shift under the influence of adsorbed molecules highly complicates the calculations if these effects are to be considered.In fact, the simulation of adsorption in zeolites demands a bridging of two relatively distant disciplines, namely chemistry of solids and molecular simulation of phase equilibria.For the equilibria, the local electrostatic charges of the molecules in vacuum (Delhommelle et al., 2000) may be determined by 2-meC 14 -simul.
Similarly, an effective research direction could be to use methods of theoretical chemistry used in chemistry of solids (DFT with breakdown into plane waves) to calculate the local charges electrostatic charges in the zeolite pores.These methods are in fact the most widely used in the field of materials with a crystalline structure (Wimmer, 1993).

CONCLUSION
Although numerical simulation at the scale of a million years was met at its inception by the skepticism of many geologists marked by the naturalistic tradition of their field, this discipline has progressively gained credibility thanks to many regional studies that have helped to compare observations and simulation results.Designated today by the term of basin models, the tools derived from this approach are routinely used in industry and in the academic world.It is extremely satisfying to have been able to contribute to the advance of this discipline through basic research that is still recognized today.With the perspective of time, we can state that this success rewards the large risks assumed by confronting difficult problems.It must be acknowledged that our team had two big chances: to have had the backing of a hierarchy motivated by the development of basin models, and to have been operational at the right time to participate fully in this scientific adventure.
If we consider thermodynamics and molecular simulation, the situation is very different in many ways: the discipline is based on solid foundations, in which many highly competent teams have been active for many years and whose potential applications cover a broader field than the oil industry.Beyond these differences, there are some significant similarities with basin models: long computer time, the need for versatile softwares to analyze realistic systems, empirical optimization of parameters from which a physical meaning is nonetheless expected.Yet the most striking common feature is the difficulty that a new technology encounters at all levels to establish itself.In terms of applications, the confidence of industry can only be secured very gradually, once many examples have been favorably investigated, and this process is especially slow because the teams employing the new methods are made of a few people at the start!No doubt there is no panacea for success in establishing a new technology, except a strong dose of perseverance!Having taken this look at the past and present, what does the future offer?In the field of basin models, the more detailed consideration of reaction mechanisms, based inter alia on a theoretical calculation of rate constants, is one path for progress that is highly interesting for better control of the kinetics of oil formation, including in cases where thermal alteration is not the only factor (thermoreduction of sulfates, for example).As to migration modeling, no doubt considerable work still needs to be done before realistic three-dimensional geometries (faults, etc.) can integrate flows involving equilibria between the water, oil and gas phases, not to speak of the physicochemical aspects which play a role in the isotopic composition of the fluids.
In the field of applied thermodynamics, interesting prospects exist in the use of new techniques, including inter alia, vibrating wire viscodensimetry used at the University of Clermont-Ferrand for fluids of high pressure.The measurement of low saturation vapor pressures with very small quantities of samples developed at the University of Lyon, the systematic analysis of petroleum fluids by GC-MS coupling.Yet it is no doubt in modeling, and particularly in molecular simulation, that the greatest prospects for progress exist.Generally speaking, the algorithms and intermolecular potentials developed in recent years have only been partly used for industrial purposes, and their practical application alone is a very fruitful pathway.And still concerning the petroleum industry and chemical engineering, three research directions can be distinguished: the calculation of properties of pure toxic substances, unstable or unavailable commercially; the simulation of complex systems in the fluid state (equilibria involving polar compounds, gas condensate, supercritical fluids) and amorphous state (polymers); and adsorption equilibria in industrial adsorbents.
Beyond the industrial application of already validated techniques, the improvement of algorithms and molecular simulation softwares also offers rich prospects.Considering Monte Carlo methods and equilibrium properties, worth noting is the establishment of parallel tempering techniques, which will ultimately help to improve a factor that clearly limits the uses of molecular simulation, namely the productivity of the engineer's work.Many other improvements are worth testing, such as the non-Boltzman sampling and the gradual insertion technique.With the explicit simulation of interface, it will be possible to predict interfacial tensions and achieve a detailed understanding of the role of surfactants (Goujon et al., 2001;Chen et al., 2002).In a nearby area, network calculations help to consider much larger systems (a few tens of nanometers) thereby helping to simulate micellar systems that are inaccessible with classical atomistic techniques.The simulation of crystalline molecular solids, such as crystals of paraffins or hydrates, could also become an attractive tool.
It would be incomplete to fail to mention the techniques of molecular dynamics among these prospects.Whether for free fluids or for confined fluids, the prediction of the transport properties by classic models is much less satisfactory than that of equilibrium properties.In this area, simulation is already capable of surprisingly accurate predictions when relaxation times are not too long.By contrast, the prediction of transport properties remains an unresolved problem when relaxation is very slow-except for the case in which the path is known and where kinetic Monte Carlo techniques applies.Increased computer power will probably not suffice this problem, and its solution rather entails the development of new algorithms.
It will be useful to terminate by stressing the importance that molecular simulation can and must assume.In quantum chemistry, the success of this discipline is such that some quantities are already obtained more reliably by calculation than by experiment.However, this must not obscure the fact that for most applications, simulation techniques are still based on empirical parameters and approximate numerical solutions.This is particularly the case of statistical thermodynamics discussed here.This means that simulation remains dependent on the experimental approach and can only substitute for it partially, after intensive validation tests.
Nor can simulation substitute for classical models, which will remain indispensable for rapidly calculating properties in process simulators where the calculations must be repeated millions if not billions of times.On the contrary, it is legitimate that simulation claims a more important place in research and development programs of the petroleum industry and chemical industries.It is striking in this respect to observe the difference in involvement of the chemical engineering laboratories in Europe and the United States, where modeling at the molecular scale is clearly felt as the most promising "frontier", in scientific as well as industrial terms.The Old Continent still has time to seize the opportunities, but it cannot procrastinate. Figure1 Figure3 Figure 5Mesh employed to simulate migration in the Mahakam river delta in Indonesia.The shale/sandstone intercalations are shown by rows of meshes touching the stratigraphic horizons(Ungerer et al., 1990). Figure7 Figure 8Phase diagram of a North Sea HP/HT gas condensate showing the appearance of a solid paraffin phase at low temperature.

Figure 10
Figure 10Relative error committed by the Peng-Robinson's equation on the density of liquid hexane at 50 MPa by means of a constant volume translation and a temperature dependent volume translation(Ungerer and Batut, 1997). Figure11

Figure 15
Figure15 Figure 16 Calculation made by M. Lagache of the heat capacity at constant volume of methane from statistical fluctuations (T = 210 K , i.e.T/Tc = 1.1).

Figure 18
Figure 18Example of configuration of a simulation box of the H 2 S rich liquid phase at 343 K and 70 bar in equilibrium with an aqueous phase.The water molecules (mauve) are either isolated, or in the form of aggregates (dimers, trimers).