Experimental and Computational Modeling of Microemulsion Phase Behavior

The phase behavior of microemulsions formed in a surfactant-brine-oil system for a chemical Enhanced Oil Recovery (EOR) application is complex and depends on a range of parameters. Phase behavior indicates a surfactant solubilization. Phase behavior tests are simple but time-consuming especially when it involves a wide range of surfactant choices at various concentrations. An efficient and insightful microemulsion formulation via computational simulation can complement phase behavior laboratory test. Computational simulation can predict various surfactant properties, including microemulsion phase behavior. Microemulsion phase behavior can be predicted predominantly using Quantitative Structure-Property Relationship (QSPR) model. QSPR models are empirical and limited to simple pure oil system. Its application domain is limited due to the model cannot be extrapolated beyond reference condition. Meanwhile, there are theoretical models based on physical chemistry of microemulsion that can predict microemulsion phase behavior. These models use microemulsion surface tension and torque concepts as well as with solution of bending rigidity of microemulsion interface with relation to surface solubilization and interface energy.


Introduction
With growing global energy demand and depleting reserves, enhanced oil recovery (EOR) has become more important. Among various EOR processes, chemical EOR has been labeled an expensive process, with overwhelming parameters needed to describe a chemical EOR process and not practical to measure every one of them [1]. The chemical formulations for chemical EOR process consist of single or a combination of alkaline, surfactant and polymer. The traditional chemical EOR processes are polymer flooding, surfactant and alkaline flooding. Over the years, different modes of chemical flood injections were devised. There are the binary mix of alkali-surfactant (AS), surfactant-polymer (SP), alkaline-polymer (AP), and alkaline-surfactant-polymer (ASP) slug [2].
The research and development effort to design a robust chemical EOR formulation tailored to a specific field is challenging and laborious. To design a successful surfactant related chemical EOR formulation, Pope [3] highlighted 8 surfactant

Microemulsion and Winsor phase behavior
Microemulsions are dispersions of oil and water stabilized by surfactant molecules [11][12][13]. They can take on many structures such as water droplets in oil, oil droplets in water, sponge like, bicontinuous structures, and lamellar phase [14]. Unlike emulsions, they are thermodynamically stable. This because of the oil-water IFT is low enough (below 10 ˗2 nM/M) to compensate the dispersion entropy. The interfacial energy is balanced by the dispersion entropy when the dispersion sizes are small enough, which is below 100 Å [15]. Emulsion has much larger dispersion sizes at approximately 1 μm. The different in size between microemulsion and emulsion explains their difference in properties and appearance; however, their fundamental difference is thermodynamic stability [9].
For a given set of conditions (temperature, composition), microemulsion displays well defined structures. The physical properties of microemulsion often undergo an abrupt change over a narrow concentration range [16]. It is generally accepted that this rapid change in the property over the range of concentration is due to the formation of surfactant aggregates or micelle in solution [9]. The nucleation of micelles is spontaneous, forming a structure that can vary between spherical and cylindrical depending on the surfactant molecular structure, solution composition and temperature. Figure 1 depicted the intermicellar equilibrium and the associated phase changes. Micellar structure, S 1 is formed when the hydrophilic group of a surfactant are in contact with water while the hydrophobic group are gathered within the interiors of micelles to create small regions from which water is essentially excluded. When aggregates of surfactant form in apolar solvent, it is called inverted micelles, S 2 . Inverted micelles promote the solubility of water in apolar solvents. Micellar aggregates have a finite mean lifetime, where their structures are mobile, the interfaces are flexible with a rapid exchange of molecules between the neighboring region and the aggregate [9]. Therefore, both S 1 and S 2 do form isotropic solutions with a bicontinuous and fluctuating structure in a given optimal condition. The microstructure of the sequence of states progressing from S 1 to S 2 is difficult to study and hence, these systems were large ignored until recent years, when it has been recognized that the isotopic solutions between S 1 and S 2 are agents which can significantly enhance oil recovery [9].
In less frequent cases, the microemulsion structure is made of elongated cylinders, eventually interconnected, or of distorted lamellar (sponge like), as depicted as M 1 and M 2 in Figure 1. These structures are encountered when the spontaneous curvature of the surfactant layer is small and approaches zero. On the contrary, sponge like structure can be significantly swollen by both oil and water. The G phase in Figure 1 is liquid crystal with lamellar structure [15].
A microemulsion can exist in three types of systems. Winsor's introduced three types of simple phase diagram that are characterized by the nature of the polyphasic zone at low and moderate concentration. They are Winsor I (WI), Winsor II (WII) and Winsor III (WIII) [17]. Winsor I and II are also commonly known as II˗ and II+. Below a certain salinity, C seu , the system is WI. Above a certain salinity, C sel , the system is WII. If the salinity is between C seu and C sel , the system is WIII. In a WIII system, the IFT is lower than WI and WII [4]. WI and WII diagrams shows a characteristic of 2-phase behavior with water or oil microemulsion in equilibrium with the excess phase. WIII diagram shows a bicontinuous microemulsion in equilibrium with both the excess water and oil phases. WIII is considered to have the best probability of recovering additional oil. WII is considered to have the second- Intermicellar equilibrium and associated phase changes [9]. best chance to recover additional oil because it shows interaction between the aqueous phase and crude oil. Even though WI demonstrates interaction between the crude oil and the aqueous phase, it is considered to have poorer oil recovery potential than WII.
Winsor's phase behavior studies in the early 1950s introduced the R ratio of interactions of the adsorbed surfactant at the interface with the neighboring oil and water molecules as a criterion to take into account along with the effects of all formulation variables found in a ternary surfactant-oil-water (SOW) system, i.e., the surfactant head and tail characteristics, the nature of the oil, the aqueous phase salinity, as well as temperature and pressure [18].
In micellar solutions, 3 distinct regions can be identified: an aqueous region, W, an oil or organic region, O, and a surfactant region, C. The variation of the dispersing tendencies at the O and W faces of the C region is expressed qualitatively by Winsor [17] as: where A CO and A CW are the interaction of surfactant molecules per unit area at the interface with oil and water respectively, A OO is the interaction between two oil molecules, and A WW is the interaction between two water molecules. Winsor described that an optimal microemulsion (Winsor III) is formed when the microstructure surface is flat, i.e., R = 1. When R < 1, there is a tendency to form oil-inwater emulsion (Winsor II or II+), whereas when R > 1, the tendency is to form water-in-oil emulsion (Winsor I or II˗). Salager [18] presented a diagram (Figure 2) to link the Winsor R ratio with the observe phase type. R < 1, R = 1 and R > 1 correspond to WI (II˗), WIII, and WII (II+) diagrams, respectively. This shows that any formulation change that alters one of the interactions indicated in the ratio can increase or decrease R. When the formulation variation is properly selected to change R from R < 1 to R > 1 or vice versa, it changes the phase behavior from WI to WII or vice versa, with an intermediate WIII three-phase behavior at R = 1 [18].

Laboratory determination of phase behavior
Phase behavior laboratory test is conducted to assess surfactant performance in generating optimal microemulsion. The test involves the use of graduated cylinders with stopper and an explosion proof oven. A typical range of surfactant concentration and salinity for the phase behavior test are 0.25 to 2.5 wt.% and 5000 to 40,000 ppm respectively. The procedure starts with adding 50% of surfactant solution in brine into the graduated cylinder with 50% of crude oil of interest. The cylinder is then mixed vigorously and aged in the explosion oven for 14 days at reservoir temperature. The phase type and volume of the mixture are observed and measured on day-14 to identify the optimum microemulsion formulation and condition. The complete phase type descriptions are given in Table 1.
Type III is considered to have the best probability of recovering additional oil. Type II is considered to have the poorest chance to recover additional oil. Type II˗ is considered to have the second-best chance to recover additional oil because it shows interaction between the aqueous phase and crude oil. Even though Type II+ demonstrates interaction between the crude oil and the aqueous phase, it is considered to have poorer oil recovery potential than Type II˗. Figure 3 illustrates the relationships among phase behavior laboratory test observation, microemulsion structure and Winsor's R ratio.

Computational simulation of phase behavior
Apart from the commonly use of phase behavior laboratory test to assess microemulsion of a particular surfactant for chemical EOR application, it is viable to Phase type Phase type descriptions

II
Two fluid envelopes exist-A bottom aqueous phase and a top oil phase. No color is visible in the aqueous phase. The crude oil and aqueous phase volumes are equal to the volumes placed in the tube. The surfactant has been driven into the crude oil and no crude oil swelling has taken place (Type II+ phase behavior).

II˗
Two fluid envelopes exist-A bottom aqueous phase and an oil phase. The bottom aqueous phase is colored due to surfactant carrying oil into the aqueous phase. The crude volume can be swollen due to the interaction with the surfactant (added and in-situ), but this is not a requirement for this designation.

III
Three or more fluid envelopes exist-A bottom aqueous phase, one or more middle emulsion phases, and a top crude oil phase. The aqueous phase can be colored with saponified acids (if alkali is presence) from the crude oil; however, this does not necessarily have to be the case.

II+
Two fluid envelopes exist-A bottom aqueous phase and a top crude oil phase. The bottom aqueous phase is clear because the surfactant (added and in-situ) resides in the crude oil phase. The crude oil phase is swollen due to surfactant carrying water into the crude oil phase. properties and finally to optimize structures and mixtures of surfactants [19]. Molecular modeling tools, in combination with recently developed intermolecular potentials, can provide precise information about microscopic phenomena and lead to accurate estimation of thermophysical properties [20][21][22][23]. Meanwhile, QSPR is an analytical method for breaking down a molecule into a series of numerical values describing its relevant chemical and physical properties. It remains as the focus of many studies aimed at the modeling and prediction of physicochemical and biological properties of molecules [24].

Quantitative structure-property relationship (QSPR)
QSPR is an approach to relate molecular descriptors with experimental values of properties based on statistical method. Its prediction accuracy is dependent on the size and quality of database and calculation of the relevant descriptors. There are many QSPR-like terms being used for more specific situations, such as Quantitative Structure-Activity Relationships (QSAR), Quantitative Structure-Toxicity Relationships (QSTR), Quantitative Property-Property Relationships (QPPR), Quantitative Sequence-Action Model (QSAM) and Quantitative Structure-Reactivity Relationships (QSRR) [25]. QSPR models have been developed to predict properties of pure surfactants only. Development of QSPR models for mixtures of surfactants is still a challenge [26]. Salager et al. [18] presented the use of Hydrophilic-Lipophilic Deviation (HLD) equations to attain optimum chemical EOR formulation for simple surfactant, oil and water system. The use of HLD concept to predict optimum surfactant formulation is a hybrid approach that combine HLD equations with experiments data, which is a QSPR approach. It was demonstrated that the phase behavior and optimum formulation can be manipulated with four main independent variables: brine salinity, oil alkane carbon number (ACN), surfactant parameter and temperature. However, these models are limited to simple system. Jin et al. [27] predicts the optimum surfactant salinity using HLD equation and measured parameters including the equivalent alkane carbon number (EACN), salinity, surfactant head dependent parameter, K surf value and surfactant characteristic curvature, C c . The work is extended to predict the IFT behavior using the Hydrophilic Lipophilic Difference-Net Average Curvature (HLD-NAC) equation of state.
Moreau et al. [28] applied the QSPR method to predict surfactant optimal salinity based on its correlation with surfactant structures. The QSPR models have been proven in reference conditions but they cannot be extrapolated to other conditions outside the application domain.
Budhathoki et al. [29] use the HLD equation to design the ratio of surfactant mixtures to form optimal microemulsion at reservoir condition. Correct surfactant head dependent parameter, K surf and the surfactant temperature dependent parameter, αT values play a crucial role in the accuracy of the HLD method. Both K and αT values can be obtained via a combination use of HLD equation and a series of phase behavior laboratory work for each individual surfactant. This approach can reduce the experimental test matrix to design optimal surfactant mixture, but it is limited to surfactants with known K surf and αT values through extensive phase behavior laboratory work.

Dissipative particle dynamics (DPD)
Many of the interesting phenomena that occur in complex fluids occur at the mesoscale, which is roughly defined as the spatio-temporal scales ranging from 10-10 4 nm to 1-10 6 ns [30]. These scales of simulation are not feasible using MD simulation. DPD is a coarse-grained type of molecular simulation technique which could reduce the length and timescale of molecular dynamics (MD) simulations, allowing simulation of large and complex system. DPD modeling method can reach large length scales by combining molecule groups into particles or beads, and longtime scales by replacing atomistic forces with soft effective forces [31]. DPD is widely popular simulation approach due to its algorithmic simplicity and huge versatility. By varying the conservative forces between coarse-grained beads, one can readily model complex fluids such as polymers, colloids, surfactants, membranes, vesicles and phase separating fluids [30]. DPD can give insights on spatial organization of surfactants, interesting mechanistic information for films evolution or trends on surface tensions regarding structure of the adsorbed tensioactive molecules at an interface [19]. However, the challenges for a successful DPD simulation is finding robust and general methods for parameterization of the simulation system [32][33]. This is an active research area with recent approach to apply machine learning for DPD parameterization [34].
A breakthrough approach by Fraaije et al. [35] demonstrated the use of surface torque analysis in simulating surfactant phase behavior with DPD, to determine the optimal brine salinity specifically. Prior to this, QSPR statistical approach have been the only known approach for decades in determining surfactant phase behavior. Rekvig et al. [14] varies the surfactant chain length and topology to investigate the effect of surfactant structure and composition of the monolayer on the bending rigidity. This work of Rekvig et al. is of particular interest, where the linking of bending rigidity to surfactant structure in predicting the stability of microemulsions is demonstrated. This is important because it is agreed that the bending rigidity is a key parameter in understanding structure and phase behavior of microemulsion [37]. Further details of the approach by Rekvig et al. in determining bending rigidity is discussed further in the next Section 4.2.2.

Surface tension analysis
The work of Fraaije et al. [35] is the only found published work for direct determination of surfactant phase behavior with theoretical foundation based on physical chemistry of microemulsion. The calculation principal is based on an observation by Helfrich [38][39] that one can calculate the surface torque density (torque) from the first moment of the molecular stress profile, provided the surface is tensionless. A positive torque implies a tendency to bend toward the oil phase and form a microemulsion with oil droplets dispersed in an aqueous phase, while a negative torque a tendency to form a microemulsion that has water droplets dispersed in an oil phase. A surface with zero torque has an indifferent tendency where the system will form the optimal or balanced microemulsion with average zero curvature. Fraaije et al. run the DPD simulation including electrostatics and ion interactions with added salt, surfactants, and oil. The relationship between mechanical coefficients and the stress profile is expressed in Eq. (2) [39][40]: where M n is the stress moment, σ is the stress tensor and z is the coordinate perpendicular to the surface.
In mathematics, a moment of a function is a specific quantitative measure, used in both mechanics and statistics, of the shape of a set of points. For example, for a set of data points representing mass, the 0-th moment is the total mass, the first moment divided by the total mass is the center of mass, and the second moment is the rotational inertia. Similarly, for a set of data points representing probability density, the zeroth moment is the total probability (i.e., one), the first moment is the mean, the second moment is the variance, and the third moment is the skewness. The n-th moment, M n , of a real-valued continuous function f(x) of a real variable about a value c is given is Eq. (3) below: Fraaije et al. presented a surface tension analysis namely Method of Moments to shows that torque can be calculated based on Eq. (3). This is done by calculating the torque of an microemulsion interface from the first moment, M 1 , only if the zeroth moment, M 0 (the interfacial tension) is zero exactly. Otherwise, both values of the neutral surface position and the interfacial tension has to be known. Meanwhile, in the tensionless limit, the value of z s is inconsequential. The simulation is run across various salinities to find the optimal salinity at which the microemulsion surface torque is zero. Note that there is no clear-cut boundary on when a positive or negative torque transitions in to a small value (around zero). Therefore, the boundaries are somewhat gradual. Fraaije et al. demonstrated that the Method of Moments is principally correct by successfully deriving known empirical coefficients of decade-old QSPR models.
It was noted in Fraaije et al. work that the torque is related to the bending rigidity and spontaneous curvature through Eq. (4): where τ is the torque, k is the bending rigidity and C 0 is the spontaneous curvature. However, their approach does not allow for direct calculation of bending rigidity or the spontaneous curvature to deduce the stiffness of the interface. Furthermore, they have yet to attempt a full treatment of the compound mixture for application in actual crude oil system.

Interface fluctuation analysis
Microemulsion structure is governed by the elastic constants, the bending rigidity, k and the saddle splay modulus, k s [41]. Bending rigidity characterizes the resistance of the interface toward bending. A low bending rigidity means large thermal undulations and low stability.
Rekvig et al. [14] used DPD to simulate surfactant monolayers on the interface between oil and water to calculate the bending rigidity by analyzing the undulation spectrum. The effect of surfactant density, chain length, adding co-surfactant and linear versus branched surfactant on bending rigidity are investigated. The results show that increase of the monolayer thickness has a larger effect on the bending rigidity than increasing the density of the layer. The bending rigidity also increases with surfactant chain length and is larger for linear than branched surfactants. Bending rigidity decrease linearly with mole fraction of short surfactants. Mixed film has a lower bending rigidity than the corresponding pure film for all mole fractions.
The work by Rekvig et al. [14] is reference with an earlier work by Goetz et al. [42] for lipid bilayer. Goetz et al. was the first to compute bending rigidity in molecular dynamics simulations. Rekvig et al. used a simple model of head, tail, water, and oil beads in DPD to capture the essential properties of ternary systems such as phase separation and adsorption. During the simulation at very low IFT, the interface is not strictly flat and undulatory waves can be observed (Figure 4).
These fluctuations of the interface are analyzed to compute bending rigidity. It is firstly done by characterizing the interface based on continuum theory [43]. This is then followed by the adoption of Helfrich [44] free energy of the interface [Eq. (5)] with local displacement from the average position of the interface to enable easier monitoring in simulations.
where H is the principal curvature of the surface at M (unique point), K is the gaussian curvature, depending on its sign, the surface will be curved like a sphere or like a saddle-splay, these number give us an idea of the surface's shape, k is the bending rigidity, k s is the splay modulus which is related to the shape of the interface (K) and c 0 is the spontaneous curvature.
Bending rigidity is obtained by analyzing the undulation spectrum of the interface. Based on hypotheses from Rekvig et al. [14] using equipartition principal and fast fourier transform (FFT) application to decomposes the undulation signal into different wave lengths, Eq. (5) is transformed into Eq. (6) below: whereh is the approximation of local displacement from the average position of the interface, q is the wave vector, 2π λ , λ is the corresponding wavelength, k B is Boltzmann constant (1.38 Â 10 23 J/K), T is absolute temperature in Kelvin, A is the interface area and γ is the interfacial tension.
Given the spectral intensity, S(q): Combining Eqs. (6) and (7), Eq. (8) is devised: Eq. (8) is used to fit the interface undulation spectrum analysis' results to estimate the bending rigidity, k. In DPD simulation, all units are normalized where k B T is unity.
Published experiment data on bending rigidity values for microemulsion system may be used as a reference for the model predicted bending rigidity values. However, such published data is scarce or not related to monolayers. Majority of the publications are generally focused on theory. Zvelindovsky [45] mentioned that the bending rigidity for a surfactant monolayer between water and oil is usually in the range of 1-20 kBT. Martínez et al. [46] performed MD simulation for SDS surfactant in dodecane and brine system at zero salinity. The associate bending rigidity is 1.3 kBT. SDS is sodium dodecyl sulfate or sodium lauryl sulfate, sometimes written as sodium lauril sulfate. Kegel et al. [47] found that the bending rigidity for SDS surfactant with alcohol in cyclohexane and brine system is around 1 kBT. Binks et al. [48] found a bending rigidity of around 1 kBT for AOT surfactant in nonane and a brine system at optimum salinity and surfactant concentration. AOT is a twin tailed, anionic surfactant with a sulfosuccinate head group stabilized as a salt by a sodium cation. It was also reported that the bending rigidity value depends on the alkane length, where bending rigidity decreases with increase in alkane length.

Monte Carlo (MC) and molecular dynamics (MD)
Both molecular modeling methods of MD and MC studies have been carried out for academic considerations, mainly on the surfactant aggregation process, instead of for industry application. This is mostly due to the use of atomistic description of the system that requires extensive computational power, which is not practical for industry application. As mention in Section 4.2.2, Goetz et al. [42] provides the first explicit connection between computer simulations with molecular resolution and elastic membrane models based on differential geometry. Goetz's method demonstrates a relationship between bending rigidity and the IFT.

Conclusion
Phase behavior of microemulsion is commonly assessed via laboratory study. These studies are straightforward but laborious especially when it involves a huge range of surfactant choices. Computational simulation is an alternative approach to provide insights into microemulsion phase behavior. There are limited computational simulation studies to predict surfactant phase behavior, whereby the widely used method since the beginning is empirical correlations as in QSPR approach.
There are very few non-empirical approaches to predict surfactant phase behavior. These approaches are based on combination of physical chemistry of microemulsion surface tension, torque and bending rigidity concepts.