David Steinmetz

Benoît Creton

Véronique Lachet

Bernard Rousseau

Carlos Nieto-Draghi

Benoit Creton

Carlos Nieto- Draghi

Simulations of interfacial tension of liquid-liquid ternary mixtures using optimized parameterization for coarse-grained models

Keywords: Dissipative Particle Dynamics, interfacial tension, liquid-liquid equilibrium, ternary diagram, parameterization, conservative force

Abstract: In this work, liquid-liquid systems are studied by means of coarse-grained Monte Carlo simulations (CG-MC) and Dissipative Particle Dynamics (DPD). A methodology is proposed to reproduce liquid-liquid equilibrium (LLE) and to provide variation of interfacial tension (IFT) as a function of the solute concentration. A key step is the parameterization method based on the use of Flory-Huggins parameter between DPD beads to calculate solute/solvent interactions.

Parameters are determined using a set of experimental compositional data of LLE following four different approaches. These approaches are evaluated and obtained results are compared to analyze advantages/disadvantages of each one. These methodologies have been compared through their application on six systems: water/1,4-dioxane/benzene, water/acetone/chloroform, water/benzene/acetic acid, water/benzene/2-propanol, water/hexane/acetone and water/hexane/2-propanol. Gibbs (NVT) ensemble CG-MC simulations have been used to check the validity of parameterization approaches for LLE reproduction. Then, CG-MC simulations in the osmotic (µ solute N solvent P zz T) ensemble were carried out considering the two liquid phases with an explicit interface. This step allows to work at the same bulk concentrations as the experimental data by imposing the precise bulk phase compositions and predicting the interface composition. Finally, DPD simulations were used to predict IFT values for different solute concentrations. Our results on variation of IFT with solute concentration in bulk phases are in good agreement with experimental data but some deviations can be observed for systems containing hexane molecules.

INTRODUCTION

Formation of a liquid-liquid interface results from the mixing of immiscible or partially miscible fluids. This physical phenomenon is characterized by the interfacial tension (IFT) which quantifies the imbalance of intermolecular forces between molecules leading to an accumulation of free energy at the interface. Characterization of liquid-liquid interfaces is a major issue for many applications involving emulsion stabilization or liquid-liquid extraction [START_REF] Crowell | Pharmaceutical applications of liquid-liquid extraction[END_REF][START_REF] Fathi Azarbayjani | Impact of Surface Tension in Pharmaceutical Sciences[END_REF][START_REF] Thomas | Enhanced Oil Recovery -An Overview[END_REF] . For instance, the chemical Enhanced Oil Recovery (cEOR) is an application in the oil and gas industry whose effectiveness depends on the water/oil IFT reduction. It consists in the injection of Alkaline/Surfactant/Polymer (ASP) combinations to mobilize the oil trapped in the reservoir by acting on capillary forces [START_REF] Thomas | Enhanced Oil Recovery -An Overview[END_REF] . The capillary number decreases as viscosity increases and/or interfacial tension decreases. The ASP formulation aims at reaching an ultra-low IFT for the brine/surfactant/crude oil systems. Due to complex involved phenomena, identification and selection of relevant surfactants is challenging, and it requires a large number of trial and error tests. Modeling tools such as molecular simulations are adapted to improve the efficiency of this process by providing information about phenomena occurring at the molecular level and at the interface [START_REF] Creton | Prediction of Surfactants' Properties using Multiscale Molecular Modeling Tools: A Review[END_REF] .

IFT for water/oil mixtures have already been intensively studied with atomistic scale simulation methods such as Molecular Dynamics (MD) [START_REF] Buhn | A molecular dynamics study of a liquid-liquid interface: Structure and dynamics[END_REF][START_REF] Natália | Interfacial Tension Behaviour of Water/Hydrocarbon Liquid-Liquid Interfaces: A Molecular Dynamics Simulation[END_REF][START_REF] Neyt | Quantitative Predictions of the Interfacial Tensions of Liquid-Liquid Interfaces through Atomistic and Coarse Grained Models[END_REF][START_REF] Papavasileiou | Predictions of water/oil interfacial tension at elevated temperatures and pressures: A molecular dynamics simulation study with biomolecular force fields[END_REF][START_REF] Wang | Molecular Dynamics Simulation of the Liquid-liquid Interface for Immiscible and Partially Miscible Mixtures[END_REF] and Monte Carlo (MC) [START_REF] Gao | Theoretical examination of hexanol-water interfaces[END_REF][START_REF] Linse | Monte Carlo simulation of liquid-liquid benzene-water interface[END_REF] . Many works [START_REF] Buhn | A molecular dynamics study of a liquid-liquid interface: Structure and dynamics[END_REF][START_REF] Natália | Interfacial Tension Behaviour of Water/Hydrocarbon Liquid-Liquid Interfaces: A Molecular Dynamics Simulation[END_REF][START_REF] Wang | Molecular Dynamics Simulation of the Liquid-liquid Interface for Immiscible and Partially Miscible Mixtures[END_REF] have shown that atomistic simulation methods provide reasonable values of IFT for binary systems in comparison to experiments. Furthermore, additional information that is difficult to access experimentally can be extracted from simulations such as interface thickness, position and orientation of molecules at interfaces, and molecular compositions. However, calculations of IFT values in agreement with experimental data required a reliable model of intra-and inter-molecular interactions. A recent study conducted by Papavasileiou et al. [START_REF] Papavasileiou | Predictions of water/oil interfacial tension at elevated temperatures and pressures: A molecular dynamics simulation study with biomolecular force fields[END_REF] have

shown that IFT values calculated from atomistic simulations are very sensitive to the chosen molecular model. They studied IFT temperature dependence in the range from 383.15 to 443.15 K at 1.83 MPa for water/oil systems using different atomistic force fields that have been primarily developed for biomolecules. They studied binary mixtures (water/toluene and water/ndodecane) and a ternary system (water/toluene/n-dodecane) where the organic phase is described with the General AMBER Force Field (GAFF) and the Lipid14 force field for toluene and ndodecane, respectively. Water molecules were modeled with three force fields: TIP4P/2005, TIP3P and SPC/E models. They showed that, TIP4P/2005 and SPC/E are responsible for a overestimation of IFT while TIP3P model leads to slightly underestimated IFT values. In addition, they also proposed to model water/hydrocarbon interactions with modified Lorentz-Berthelot combining rules by introducing a correction term, , fitted on experimental data.

They showed that this correction term improve the accuracy of the IFT predictions using all water models for binary systems. However, they were not able to provide accurate reproduction of experimental IFT values for the ternary system. Recently, Ghoufi et al. [START_REF] Ghoufi | Computer modelling of the surface tension of the gas-liquid and liquid-liquid interface[END_REF] presented a review about molecular simulations in which different methodological factors (size-effects, truncation procedures, long-range corrections and potential models) were analyzed to bring out their effects on interfacial systems and IFT values.

Atomistic scale simulation methods have been used many times to represent systems that are more complex than binary or ternary water/hydrocarbon mixtures such as crude oil fractions 13- 19 . For example, systems including asphaltene compounds which have very complex molecular structure have been modeled in order to study the nature of interactions and to calculate interaction energies of asphaltene/asphaltene systems [START_REF] Costa | Aggregation of model asphaltenes: a molecular dynamics study[END_REF][START_REF] Gao | Molecular Dynamics Simulation: The Behavior of Asphaltene in Crude Oil and at the Oil/Water Interface[END_REF][START_REF] Greenfield | Molecular modelling and simulation of asphaltenes and bituminous materials[END_REF][START_REF] Murgich | Molecular Simulation and the Aggregation of the Heavy Fractions in Crude Oils[END_REF][START_REF] Sedghi | Atomistic Molecular Dynamics Simulations of Crude Oil/Brine Displacement in Calcite Mesopores[END_REF] . However, representation at the atomistic scale of a liquid-liquid system containing an explicit interface and bulk phases requires a large length scale and, therefore, a high computational cost. Furthermore, time scales accessible to atomistic simulations are too short to observe phenomena that are often studied in liquid-liquid systems such as formation of micelles or micelle reorganization and their diffusion to the interfaces when surfactant molecules are added [START_REF] Rekvig | Investigation of Surfactant Efficiency Using Dissipative Particle Dynamics[END_REF][START_REF] Marrink | Molecular Dynamics Simulation of the Kinetics of Spontaneous Micelle Formation[END_REF] .

In the case of large and complex systems, the numerical study of interfaces requires the use of mesoscopic simulation techniques based on a coarse-grained model (CG). CG model consists of grouping atoms or molecules into particles in order to reduce the number of degrees of freedom of the system and, therefore, the computational time. Several works in the literature deal with liquid-liquid interfaces described using the MARTINI force field [START_REF] Neyt | Quantitative Predictions of the Interfacial Tensions of Liquid-Liquid Interfaces through Atomistic and Coarse Grained Models[END_REF][START_REF] Ndao | Coarse-Graining the Liquid-Liquid Interfaces with the MARTINI Force Field: How Is the Interfacial Tension Reproduced?[END_REF] . This force field is based mainly on a four-to-one mapping, it means that four heavy atoms plus associated hydrogen atoms are grouped into a single particle. In the case of water, four water molecules are grouped into a particle. Alternatively to the MARTINI model, the Dissipative Particle Dynamics (DPD) seems to be a promising molecular simulation method. In DPD simulations, the coarse-grained model can be easily adapted according to the phenomena studied. For example, the study of interactions between asphaltenes [START_REF] Zhang | Aggregate Structure in Heavy Crude Oil: Using a Dissipative Particle Dynamics Based Mesoscale Platform[END_REF] or orientation [START_REF] Ruiz-Morales | Coarse-Grained Molecular Simulations to Investigate Asphaltenes at the Oil-Water Interface[END_REF] of these molecules at interfaces requires small coarse-grained levels (i.e. water particle corresponds to three water molecules) while emulsion phenomena [START_REF] Alvarez | Dissipative Particle Dynamics (DPD) Study of Crude Oil-Water Emulsions in the Presence of a Functionalized Co-polymer[END_REF][START_REF] Song | Investigation of demulsification efficiency in water-in-crude oil emulsions using dissipative particle dynamics[END_REF] are simulated using large coarse-grained levels (i.e. water particle corresponds at least to eighteen water molecules).

Although the CG model decreases the level of detail on molecular structure compared to atomistic model, CG simulations have already been used to predict quantitative values of the IFT for binary mixtures [START_REF] Neyt | Quantitative Predictions of the Interfacial Tensions of Liquid-Liquid Interfaces through Atomistic and Coarse Grained Models[END_REF][START_REF] Ndao | Coarse-Graining the Liquid-Liquid Interfaces with the MARTINI Force Field: How Is the Interfacial Tension Reproduced?[END_REF][START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF][START_REF] Goel | Estimation of interfacial tension for immiscible and partially miscible liquid systems by Dissipative Particle Dynamics[END_REF][START_REF] Maiti | Bead-bead interaction parameters in dissipative particle dynamics: relation to bead-size, solubility parameter, and surface tension[END_REF] . It has been shown that MARTINI model allows a prediction of the IFT within ±10 mN/m with respect to experiments for organic-water systems such as alkanes/water, benzene/water or chloroform/water [START_REF] Ndao | Coarse-Graining the Liquid-Liquid Interfaces with the MARTINI Force Field: How Is the Interfacial Tension Reproduced?[END_REF] . However, additional calibration of the parameters on the basis of liquid-liquid properties seems necessary for more quantitative results [START_REF] Neyt | Quantitative Predictions of the Interfacial Tensions of Liquid-Liquid Interfaces through Atomistic and Coarse Grained Models[END_REF] .

Using DPD simulations, Goel et al. [START_REF] Goel | Estimation of interfacial tension for immiscible and partially miscible liquid systems by Dissipative Particle Dynamics[END_REF] have estimated the interfacial tension for a large number of immiscible and partially miscible systems. The percentage of error obtained on the value of IFT is usually a few percent (<10% However, the representation of multi-constituent systems within molecular simulations remains a difficult task. Indeed, IFT values calculated with molecular simulations can be compared with experimental data only if compositions of bulk phases are well reproduced in the simulation boxes. Typically, when a solute is partially miscible in the two bulk phases of a liquid-liquid equilibrium, partition coefficient of this compound must be perfectly reproduced, thus, the calculated value of IFT can be compared to the experimental value for a given concentration of solute. Therefore, parameterization of CG simulations for multi-constituent systems requires that chemical interactions must be precisely described to model forces governing the IFT and to take into account solubility and miscibility of chemical species. To the best of our knowledge, a reliable method to fulfill these requirements is not available in the literature.

The key point of this work is to develop a parameterization strategy for DPD simulations involving multi-constituent systems including partially miscible solutes and providing quantitative predictions of the interfacial tension. We propose four parameterization approaches to reproduce ternary liquid-liquid systems. These approaches are compared with each other using different simulation techniques based on DPD model to reproduce miscibility, IFT and to investigate interfacial composition. The manuscript is organized as follows: in section 2, we describe theoretical aspects about the DPD simulation method, and we present existing approaches for DPD parameterization and discuss their limits. Then, alternative parameterization methods are proposed to take into account the miscibility and solubility. In section 3, the simulation results obtained using different parameterization approaches are presented for six ternary systems and the ability of these approaches to reproduce liquid-liquid equilibrium and predict IFT is discussed. Finally, section 4 gives our conclusions and perspectives. Hoogerbrugge and Koelman to model the hydrodynamic behavior of complex fluids at the mesoscopic scale [START_REF] Hoogerbrugge | Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics[END_REF][START_REF] Koelman | Dynamic Simulations of Hard-Sphere Suspensions Under Steady Shear[END_REF] . DPD method is based on a coarse-grained model in which atoms or molecules are grouped within DPD "beads" in order to reduce the number of degrees of freedom.

Beads interact with each other, and time evolution is governed by the Newton's laws:

,

,

where , and are the position, velocity, and force applied to the ith bead, respectively. It should be noted that in the DPD model, all masses, m i , are set to unity (reduced unit). The total force exerted on a bead i is defined as the sum of conservative , dissipative and random forces, as expressed in equation (3). To simulate molecular skeleton, typically long hydrocarbons or amphiphilic molecules, an intramolecular force is added to bind two neighboring beads:

.

The conservative force, , describes pairwise soft repulsive interactions between beads:

1 , (4)
where is the interaction parameter representing the maximum repulsive magnitude between ith bead and jth beads. is the cut-off radius and represents the maximum range of interactions.

The unit vector pointing from the jth bead to the ith bead is defined by /| | with and .

The dissipative force, , can be interpreted as a frictional force representing the viscosity effects and is defined as:

• , (5)
where is the dissipative weight function, is the friction coefficient and is the velocity difference between the ith bead and the jth bead, . To avoid the system freezing, energy is injected through the random force, , as follows:

, ( 6
)
where is the random weight function, is the random force amplitude and is a random number with zero mean and unity variance when averaged over time.

Español and Warren [START_REF] Español | Statistical Mechanics of Dissipative Particle Dynamics[END_REF] have shown that to ensure the simulation obeys a canonical ensemble, the weight functions and as well as the friction coefficient and the random force amplitude are related by the fluctuation-dissipation theorem as shown in equations (7) and (8).

2 , (7) 1
, 0, , (8)
where is the Boltzmann constant and the temperature. Finally, the intramolecular force, , can be expressed as an ideal spring: ,

where is the spring constant and is the equilibrium spring distance. In this work, intramolecular forces are used for alkane molecules and parameters are 100 (DPD unit) and 0.7 . It should be noted that for amphiphilic molecules, Goicochea et al. [START_REF] Goicochea | Dependence of thermodynamic properties of model systems on some dissipative particle dynamics parameters[END_REF] and more recently Deguillard et al. [START_REF] Deguillard | Interfacial tension in oil-watersurfactant systems: on the role of intra-molecular forces on interfacial tension values using DPD simulations[END_REF] have shown that parameters of intramolecular force may influence the values of water/hydrocarbon interfacial tension.

b. State of the art of parameterization of the conservative force Parameterization of the conservative force, , has been extensively studied to properly mimic interactions between beads in order to reproduce some specific properties such as interfacial tension [START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF][START_REF] Goel | Estimation of interfacial tension for immiscible and partially miscible liquid systems by Dissipative Particle Dynamics[END_REF][START_REF] Maiti | Bead-bead interaction parameters in dissipative particle dynamics: relation to bead-size, solubility parameter, and surface tension[END_REF] , critical micelle concentration [START_REF] Vishnyakov | Prediction of the Critical Micelle Concentration of Nonionic Surfactants by Dissipative Particle Dynamics Simulations[END_REF] or radial distribution functions [START_REF] Sepehr | Dissipative Particle Dynamics interaction parameters from ab initio calculations[END_REF] .

Parameterization can be performed using two different approaches. Approaches gather under the term top-down derive parameters from macroscopic properties (isothermal compressibility, interfacial tension and so on) whereas in approaches called bottom-up parameters are determined on the basis of atomistic configurations (using potential of mean force [START_REF] Trément | Conservative and dissipative force field for simulation of coarse-grained alkane molecules: a bottom-up approach[END_REF][START_REF] Trément | Conservative and dissipative force field for simulation of coarse-grained alkane molecules: A bottom-up approach[END_REF] , iterative Boltzmann inversion [START_REF] Lahmar | Onset of Entanglements Revisited. Dynamical Analysis[END_REF] or ab initio calculations 38 methods). As far as we know, bottom-up approaches to parameterize DPD simulations have never been used to reproduce liquid-liquid systems or to compute interfacial tension but several works can be found in the literature where hydrodynamics properties for polymeric systems are studied [START_REF] Trément | Conservative and dissipative force field for simulation of coarse-grained alkane molecules: a bottom-up approach[END_REF][START_REF] Trément | Conservative and dissipative force field for simulation of coarse-grained alkane molecules: A bottom-up approach[END_REF][START_REF] Lahmar | Onset of Entanglements Revisited. Dynamical Analysis[END_REF][START_REF] Lahmar | Influence of the adjustable parameters of the DPD on the global and local dynamics of a polymer melt[END_REF] . Therefore, bottom-up approaches are not presented hereafter and, we will thus focus only on top-down approaches.

The most widely used method to determine interaction parameters between two identical beads (like beads) has been developed in 1997 by Groot and Warren [START_REF] Groot | Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation[END_REF] . They proposed a relationship between the isothermal compressibility and the interaction parameter between like beads as shown in equations ( 10) and (11).

1 2 ̅ , (10)
with

1 , (11)
where is the dimensionless isothermal compressibility, ̅ is the DPD number density, is the number density of molecules and is a constant equal to 0.101 as proposed by Groot and

Warren [START_REF] Groot | Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation[END_REF] . Note that in the DPD model, interaction parameters are temperature dependent. The dimensionless isothermal compressibility of water at ambient conditions with ̅ 3 is approximately equal to 16 which leads to the value 25 commonly used in the literature. However, alternative formulations introduced the degree of coarse-graining, , which is the number of water molecules grouped within a DPD bead, to express the interaction parameter between like beads [START_REF] Mayoral | Modeling the temperature dependent interfacial tension between organic solvents and water using dissipative particle dynamics[END_REF][START_REF] Mayoral | Study of interfacial tension between an organic solvent and aqueous electrolyte solutions using electrostatic dissipative particle dynamics simulations[END_REF] :

1 2 ̅ . ( 12
)
For example using 3, the interaction parameter between two water beads is

78.3

. In this work, the two approaches will be compared in order to highlight the effects of parameters on the reproduction of the miscibility of compounds and the variation of IFT on liquid-liquid equilibria.

Interaction parameters between different particles (unlike beads) have been derived from several macroscopic properties in the literature. Groot and Warren [START_REF] Groot | Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation[END_REF] have related interaction parameters with the Flory-Huggins parameters (χ) from the Flory-Huggins (FH) theory of polymers solutions [START_REF] Flory | Thermodynamics of High Polymer Solutions[END_REF][START_REF] Huggins | Solutions of Long Chain Compounds[END_REF] . However, this approach is based on an important approximation, it requires that all liquids have approximatively the same isothermal compressibility ( ). In this way, interaction between like beads can be considered as reference energy, so that interactions between unlike beads can be expressed in terms of excess energy ∆ compared to the reference:

∆ . (13)
Then, Groot and Warren have shown that ∆ is proportional to χ, leading to the following expression [START_REF] Groot | Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation[END_REF] :

• , (14)
where is a constant that depends on the type of system. The value of is approximately equal to 3.5 at the density ̅ 3 and for systems consisting of beads without any intramolecular force (monomeric mixture) [START_REF] Groot | Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation[END_REF][START_REF] Wijmans | Phase behavior of monomeric mixtures and polymer solutions with soft interaction potentials[END_REF] . The connection between interaction parameter with the FH theory is convenient because χ parameter values are available in the literature for many systems or it can be experimentally measured [START_REF] Clarke | Measurements of the Flory-Huggins Interaction Parameter for Polystyrene-Poly(4-vinylpyridine) Blends[END_REF][START_REF] Mumby | Determination of .chi. from liquid-liquid phase data and the computation of phase diagrams for quasi-binary polymer solutions and blends[END_REF] or calculated [START_REF] Ruiz-Morales | Coarse-Grained Molecular Simulations to Investigate Asphaltenes at the Oil-Water Interface[END_REF][START_REF] Maiti | Bead-bead interaction parameters in dissipative particle dynamics: relation to bead-size, solubility parameter, and surface tension[END_REF][START_REF] Alasiri | Dissipative particle dynamics (DPD) study of the interfacial tension for alkane/water systems by using COSMO-RS to calculate interaction parameters[END_REF] or derived from atomistic simulations [START_REF] Ryjkina | Molecular Dynamic Computer Simulations of Phase Behavior of Non-Ionic Surfactants[END_REF][START_REF] Sun | COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase ApplicationsOverview with Details on Alkane and Benzene Compounds[END_REF] .

Travis et al. [START_REF] Travis | New parametrization method for dissipative particle dynamics[END_REF] have related interaction parameters between unlike beads with solubility parameters (δ) using the Regular Solution Theory (RST) introduced by Scatchard-Hildebrand [START_REF] Scatchard | Equilibria in Non-electrolyte Solutions in Relation to the Vapor Pressures and Densities of the Components[END_REF][START_REF] Hildebrand | The Derivation of Equations for Regular Solutions[END_REF] :

2 , (15)
where and are the solubility parameters of the ith and jth beads, respectively. The dimensionless equation can be expressed as:

̅ ̅ ̅ 2 , (16)
where ̅ and ̅ are the dimensionless solubility parameters of the ith and jth beads, respectively, and , and are the dimensionless interaction parameters ( ). This approach has the advantage of removing the assumption of identical repulsions between like beads.

Recently, Vishnyakov et al. [START_REF] Vishnyakov | Prediction of the Critical Micelle Concentration of Nonionic Surfactants by Dissipative Particle Dynamics Simulations[END_REF] proposed a thermodynamic approach in which interaction parameters between unlike beads are obtained by fitting infinite dilution activity coefficients for binary solutions. It should be noted that this approach requires the same approximation used by Alasiri and Chapman 51 have shown that infinite dilution activity coefficients can be related to the χ parameters, thus, interaction parameters are obtained using equation ( 14). This approach has been validated on interfacial tension for water/alkane systems.

Anderson et al. [START_REF] Anderson | Dissipative particle dynamics: Systematic parametrization using water-octanol partition coefficients[END_REF] proposed a parameterization scheme to determine interaction parameters based on water-octanol partition coefficients and liquid phase densities. Interaction parameters between unlike beads are tuned, step by step, using a brute-force approach, in order to reproduce the experimental partition coefficient. In addition, at each step, the interaction parameters between like beads and the cutoff radius for each interaction pair are optimized to fit phase densities. This parameterization method is noteworthy because, to our knowledge, cut-off radius optimization has never been discussed to parameterize DPD simulations in order to reproduce experimental data.

Quantitative data of IFT are generally obtained from DPD simulations by parameterizing interactions between unlike beads using Hildebrand solubility parameters [START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF][START_REF] Goel | Estimation of interfacial tension for immiscible and partially miscible liquid systems by Dissipative Particle Dynamics[END_REF][START_REF] Maiti | Bead-bead interaction parameters in dissipative particle dynamics: relation to bead-size, solubility parameter, and surface tension[END_REF][START_REF] Alasiri | Dissipative particle dynamics (DPD) study of the interfacial tension for alkane/water systems by using COSMO-RS to calculate interaction parameters[END_REF] . For example,

Rezaei and Modarress 27 have compared the approach of Travis et al. with that of Groot and

Warren on IFT reproduction for several water/hydrocarbon systems. In the latter case, the Flory-Huggins parameters have been calculated using Hildebrand solubility parameters as shown in equation (17).

,

where is the mean volume of a bead. Rezaei and Modarress [START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF] have shown that parameterization of interactions between unlike particles using Hildebrand solubility parameters with both approaches leads to quantitative values of IFT in good agreement with experimental data. However, Hildebrand solubility parameters are not suitable for interaction between polar molecules or mixtures involving hydrogen bonds. Consequently, only immiscible binary systems such as water/hydrocarbon can be handled. To deal with more complex systems, alternative approaches exist using Hansen solubility parameters [START_REF] Hansen | Solubility Parameters: A User's Handbook[END_REF] . According to Hansen, the solubility parameter can be divided into three components: a term for dispersion , a term for polarity and a term for hydrogen bonding . Thus, the Flory-Huggins parameter can be expressed as:

, , 0.25 , , 0.25 , , , (18)
where 1 as suggested by Hansen. This approach has been used by Shi et al. [START_REF] Shi | Dissipative particle dynamics study of the water/benzene/caprolactam system in the absence or presence of non-ionic surfactants[END_REF] for the water/benzene/caprolactam system. However, variation of the IFT calculated with DPD simulations as a function of the caprolactam concentration differs from experimental data. It can be possible that differences between simulation and experimental data result from an insufficient consideration of partial miscibility between compounds. The study of IFT variation as a function of the solute concentration can be performed only if the phase composition is well reproduced in the DPD simulations.

c. New parameterization approach for liquid-liquid equilibrium of ternary systems

In this work, we propose a new methodology to parameterize interactions between unlike beads in order to reproduce the liquid-liquid equilibrium of ternary systems. Interaction parameters for ternary systems labeled solvent 1/solvent 2/solute are obtained as follows:

-Interactions between like particles. The corresponding interaction parameters are calculated from the isothermal compressibility of fluid as proposed by Groot and Warren.

solvent 1/solvent 2 interactions. This interaction is calculated from the Hildebrand solubility parameters. As shown previously, this approach provides quantitative values of IFT in agreement with experimental data for water/hydrocarbon systems.

-Solute/solvent 1 and solute/solvent 2 interactions. These parameters must be determined precisely because they govern the solubility of the solute in bulk phases. In this work, we

propose a new parameterization method using compositional data to determine solute/solvent interactions.

The parameterization of solute/solvent interactions is based on the assumption that solvents which constitute the bulk phases are totally immiscible. The solute is considered to be miscible in both bulk phases. For a solvent 1/solvent 2/solute system, this means that each bulk phase represents a binary system. The first one consists of solvent 1 and solute molecules while the second phase is a mixture of solvent 2 and solute molecules (see Figure 1). to use compositional data in the Flory-Huggins equations. In the FH theory, the Helmholtz free energy of mixing, ∆ , between two molecules A and B can be written according to the equation (19).

∆ ∅ ∅ ∅ ∅ , ∅ ∅ , (19)
where ∅ and are respectively the volume fraction and the number of beads constituting the molecule A, and ∅ and are the volume fraction and the number of beads constituting the molecule B. The chemical potential of compounds A and B are given by equations ( 20) and ( 21), respectively.

μ μ ° ln 1 ∅ 1 1 ∅ , ∅ , (20)
μ μ ° ln 1 ∅ 1 1 ∅ , ∅ , (21)
where μ ° and μ ° are pure compound chemical potentials of molecules A and B, respectively. By applying the FH theory for a liquid-liquid equilibrium on the solvent 1/solvent 2/solute system shown in Figure 1, equation of the chemical potential of the solute, , can be written as function of solvent 1 and solvent 2. At thermodynamic equilibrium, the chemical potential of the solute is equal in the two liquid phases:

μ μ . (22)
Thus, equation ( 22) can be used to link, , , the parameter between solute and solvent 1 with, , , the parameter between solute and solvent 2 as shown in equation (23).

ln 1 ∅ 1 1 ∅ , ∅ ln 1 ∅ 1 1 ∅ , ∅ , (23)
where , and is the number of beads constituting the solvent 1 molecules, solvent 2 molecules and solute molecules, respectively. Finally, χ parameters are related with DPD interaction parameters using equation (14). Note that the use of Flory-Huggins parameters and equation (14) requires that the interactions between like beads are equal. In our work, we proposed two approaches to determine the / parameters:

1) The first approach (#1) combines Hildebrand solubility parameters and experimental data for only one liquid-liquid composition. Thus, one of the two χ parameters is calculated with solubility parameters while the second is obtained by solving equation ( 23). This approach is valid only if one of the solute/solvent interaction parameters can be calculated from the Hildebrand solubility parameters.

2) The second approach (#2) consists in using different experimental composition data, with 2. A system of equations is established with equations ( 23) and two unknowns: , and , . An optimal solution of this system can be obtained with the least squares method.

The , and , parameters are obtained by minimizing the objective function presented in equation (24).

, ;

with

ln 1 ∅ 1 1 ∅ ln 1 ∅ 1 1 ∅ ∅ . (25)
In this work, it is assumed that the χ parameter can be considered as a constant over the composition range studied. A study of the dependence of χ parameter with compositions is carried out in section 3.a. to discuss this assumption.

The use of the Flory-Huggins theory as proposed in approaches #1 and #2 for the parameterization of solute/solvent interactions requires identical self-repulsion between like beads. This assumption may be a limiting factor for reproducing compositional data or for the prediction of IFT. Two other approaches are proposed to clarify this assumption:

3) In the third approach (#3), interactions between like beads are calculated from isothermal compressibility of each fluid using the equation (10). Because parameterization of solute/solvent interactions using equation ( 14) is no longer valid due to the loss of symmetry (i.e.

), one of the solute/solvent parameter is calculated with Hildebrand solubility parameter following the model developed by Travis et al. [START_REF] Travis | New parametrization method for dissipative particle dynamics[END_REF] (see equations ( 15) and ( 16)), and the second solute/solvent parameter is obtained with a direct fitting performed on available compositional data. (see supporting information for additional details).

4) In the fourth approach (#4), interactions between like beads are also calculated from isothermal compressibility of each fluid but the degree of coarse-graining, , is taken into account using equation (12). Solute/solvent parameters are calculated in the same way as in the #3 approach.

A summary of the specificities and input data used for each parameterization approach is presented in Table 1.

Table 1. Summary of input data used for parameterization approaches developed in this work. a ij is the repulsive parameter,  T is the dimensionless isothermal compressibility,  i the Hildebrand solubility parameter and N m is the number of water molecules in one bead.

Input data to compute interaction parameters ensemble, chemical potential of each species are calculated using Widom insertion test method [START_REF] Widom | Some Topics in the Theory of Fluids[END_REF] .

-Some CG-MC simulations are also performed in the osmotic (µ solute N solvent P zz T) ensemble in order to describe a system with a constant number of solvent particles (N solvent ) and a variable number of solute particles, fixing the chemical potential of the solute (µ solute ). Osmotic ensemble has already been used by Rekvig et al. [START_REF] Rekvig | Investigation of Surfactant Efficiency Using Dissipative Particle Dynamics[END_REF] -Finally, DPD simulations are performed in the NVT ensemble in order to compute the interfacial tension.

In this work, CG-MC in the Gibbs (NVT) ensemble and DPD simulations are performed at constant density ( ̅ 3). Therefore, the total pressure of the system varies depending on the composition. This choice was made to simplify the parameterization procedure. However, it is important to notice that an alternative parameterization of DPD interactions can be done by working at constant pressure. In this case, calculation of interaction parameters between like and unlike beads depends on the total density of the system (see equations ( 10), ( 14) and ( 15)).

Consequently, additional bulk phase density data is required to obtain interaction parameters of pure components. Noting that on this basis, any phase property should be calculated at constant pressure using a reference value.

CG-MC simulations were carried out with the molecular simulation package GIBBS [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry: Monte-Carlo Methods[END_REF] . For simulations in the Gibbs (NVT) ensemble, the two subsystems ("boxes") have each an initial dimension of 10 (in DPD units). The volume of each boxes can vary during the simulation but the total volume remains constant. The total number of beads is 6000. For simulations in the osmotic (µ solute N solvent P zz T) ensemble, box dimensions were set to 60, 10 (in DPD units). Two planar water-organic compound interfaces are created normal to the z-axis. Box length in z-direction is six times larger than in the x and y-directions in order to avoid interactions between the two interfaces. Initial boxes containing a total of 18 000 DPD beads are built for different solute/solvent concentrations using the PACKMOL software package [START_REF] Martínez | PACKMOL: a package for building initial configurations for molecular dynamics simulations[END_REF][START_REF] Martínez | Packing optimization for automated generation of complex system's initial configurations for molecular dynamics and docking[END_REF] . DPD simulations in the NVT ensemble were carried out with the molecular dynamics simulation package NEWTON [START_REF] Nguyen | Viscosity of the 1-ethyl-3methylimidazolium bis(trifluoromethylsulfonyl)imide ionic liquid from equilibrium and nonequilibrium molecular dynamics[END_REF] . Initial configurations are derived from simulations in the osmotic (µ solute N solvent P zz T) ensemble. The area of the interface is kept constant ( 10, in DPD units). A modified version of the velocity-Verlet algorithm [START_REF] Groot | Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation[END_REF] governed the equation of motion, and the time step is fixed at δt = 0.001 in DPD units. Constants in the dissipative force γ and random force σ were set to 4.5 and 3, respectively, in order to keep the temperature fixed at 1, thus satisfying the fluctuation-dissipation theorem (equations ( 7) and ( 8)). In all simulations, periodic boundary conditions were imposed in all directions.

The IFT is evaluated with the Irving-Kirkwood 64 method where the simulation box is divided into a number of slabs parallel to the interface. Then, the value of the IFT is obtained by integrating the difference between normal and tangential stresses across in the zdirection as shown in equation ( 26).

1 2 , (26)
where the factor of 1 2 accounts for two existing interfaces in the simulation box. The coefficient used to convert IFT values from DPD to the international system of units is .

The representation of ternary systems with the coarse-grained model follows the procedure reported by Goel et al. [START_REF] Goel | Estimation of interfacial tension for immiscible and partially miscible liquid systems by Dissipative Particle Dynamics[END_REF] . Beads are chosen to represent the same molecular volume. This procedure leads to more accurate results on the interfacial tension [START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF] For each studied system, solute molecule is able to form a hydrogen bond with water. But the type and strength of these interactions are different. Indeed, 1,4-dioxane and acetone molecules are acceptors of hydrogen bonds. In contrast, 2-propanol is mostly a hydrogen bonding donor, while acetic acid has the particularity of being both an acceptor and a donor of hydrogen bonds with water. These solute/water hydrogen bond interactions must be implicitly represented in the DPD simulation in order to both reproduce the miscibility and the solubility of molecules. On the other side, solutes do not form hydrogen bonds with the molecules of the organic solvents except for the chloroform with acetone. Therefore, when parameterization approaches #1, #3 and #4 are followed, it seems more reasonable to calculate the solute/organic solvent interactions with Hildebrand solubility parameters and compute solute/water using compositional data.

Coarse-grained representation of these systems and their characteristics are shown in Figure 2.

To construct coarse-grained model and to parameterize DPD simulations, molecular volume and Hildebrand solubility parameters for each component are required. Values used in our work are given in the Table 2 and were calculated using data extracted from the DIPPR database 68 .

Molecular volumes were obtained from molar volumes, and Hildebrand solubility parameters were calculated using molar vaporization enthalpies ΔH vap according to equation (27).

∆ , (27)
where is the molar volume of molecule i. The right-hand side of equation ( 27) assumes that vapor phase can be considered as an ideal gas. phase for each pair of compositions for the water/benzene/acetic acid system (obtained using approach #2). It is clear that for these systems the assumption of constant parameters over the whole composition range studied is not valid using approach #1. However, for systems containing acetone solute, variations of the / parameter seem to be more reasonable and less dependent on the composition (see Figure S4 and Figure S6 in supporting information).

For example, the / parameter varies between 1.37 and 1.60 for the water/chloroform/acetone system (i.e.

/ varies from 29.83 to 30.61). data are used for the minimization of equation ( 24) following approach #2. Figure 3 (b) shows that approach #2 seems to provide less compositional dependent / parameters than approach #1 for the water/benzene/acetic acid system. Indeed, the value of / parameters seems to be relatively constant regardless the compositions used in the optimization process. Although, the pair of compositions with the lowest solute concentration provide slightly lower parameters than the others. The same trend can be observed for the water/benzene/1,4dioxane for / parameters (see Figure S3 in supporting information). For systems containing hexane molecules, it can be noted that / parameters can vary according to the pairs of compositions used, but the value of these parameters remain around an average value (the value of parameters when all compositions are used in approach #2, Figure S6 andS7).

On one hand, it appears that approach #2 provide less compositional dependent / parameters than the approach #1 for the water/benzene/1,4-dioxane, water/benzene/acetic acid and water/hexane/2-propanol systems. On the other hand, approach #1 seems better suited to water/chloroform/acetone and water/hexane/acetone systems.

Therefore, assumption of constant parameters over the whole composition range is therefore valid for these systems depending on the approach used. Among all the systems studied in this work, only the water/benzene/2-propanol cannot be parameterized using both approaches in order to obtain parameters that are weakly dependent on the composition (see figure S5).

b. Liquid-liquid equilibrium (LLE)

1 ∑ ∑ ̅
with ̅ ∑ where is the molar fraction of solute in the aqueous phase and in the organic phase for the ith composition. N is the number of compositions and a factor 2 is added because the molar fraction of the solute is taken into account for both bulk phases.

[b]

1 ∑ ∑ with ∑
where is the value of the interfacial tension for the ith composition. To calculate , composition in absence of solute is included.

[c] water/acetone interaction parameter from water/chloroform/acetone system (approach #2) is used for the water/hexane/acetone system. [d] Due to a phase separation between the 2-propanol and the organic phase, IFT for the five most concentrated compositions in solute cannot be computed. Experimental data are plotted in black (solid lines) and the results from the CG-MC simulations are in red (dashed lines). Dimensionless DPD interaction parameters used to compute LLE are located below each ternary plot. Numbers in red represent parameters obtained by fitting to the experimental data Data presented in Figure 4 come from the equilibration of a heterogeneous system that decants, following the conodal lines (red dashed line), into two phases represented on the ternary diagram by red squares. Initial compositions are given in the supporting information. In the approach #1, with a high concentration of 1,4-dioxane, the molar fraction of solute is well reproduced in the organic phase but is overestimated in the aqueous phase. In the parameterization approach #2, conodal lines are better reproduced but an increasingly excess of 1,4-dioxane is observed in the organic phase when the concentration of this solute increases. In approaches #1 and #2, a deviation can be observed on the molar fractions of water in benzene with a high solute concentration. Differences can be explained by the main assumption of our parameterization method: the miscibility between the two solvents is not taken into account for the calculation of the solvent/solvent interaction parameters. However, with the approaches #3 and #4, when the compressibility of each fluid is taken into account, the miscibility between water and benzene appears to be slightly better reproduced. Other systems also seem to benefit from the use of different interaction parameters between like beads. The parameterization approach #3 roughly provides the best predicted miscibility between the aqueous and the organic phase.

Performances of CG-MC simulations for the water/chloroform/acetone system are similar to those obtained for the water/benzene/1,4-dixoane system. LLE diagrams are consistent with experimental data when using the four parameterization approaches but the same errors on the calculated compositional data are observed. Molar fractions of acetone in organic phase are overestimated when the solute concentration increases. In addition, solubility of water in the organic phase is not sufficiently well reproduced in our simulations. It can be noted that with the approach #2, the chloroform/acetone interaction parameter is lower than 25, this value denotes a strong attraction between the two compounds, it is important to emphasize that values lower than 25 cannot be obtained using Hildebrand solubility parameters (see equation ( 17)).

For the water/benzene/acetic acid system, LLE diagram parameterized with the approach #2 is in good agreement with experiments, with some deviations in the organic phase. However, LLE diagrams calculated using interaction parameters derived from approaches #1, #3 and #4 are not in agreement with experimental data. In these approaches, interaction parameter / is determined from Hildebrand solubility parameters and it does not depend on compositional data. And as shown in the previous section 3.a., setting interaction parameter between acetic acid and benzene to a constant value leads to a strong dependence on the composition of the acetic acid/water parameter. Using approaches #1, #3 and #4, we arbitrarily selected the second less concentrated composition for the parameterization.

For the water/hexane/acetone system, the parameterization approaches #1 and #2 provide similar solute/solvent interaction parameters and, therefore, results on phase compositions are also similar (results are grouped and noted [#1-#2] in Table 3). Phase compositions obtained with the four parameterization approaches are in very good agreement with the experimental data. The miscibility between the aqueous and organic phases is particularly well reproduced with the approach #4. For this system, we propose a fifth parameterization approach based on the principle of transferability of interaction parameters and called Transferability in Table 3.

Indeed, water/acetone parameter with a water bead containing four water molecules has already been calculated for the water/chloroform/acetone system. Therefore, water/acetone interaction parameter ( / 28.20 with the approach #2) from the water/chloroform/acetone system is used in the Transferability approach and the second solute/solvent parameter ( / ) is directly obtained with equation (23). Results of our simulations show that this method is valid only for compositions with a low acetone molar fraction in bulk phases (i.e., 0.2 in the aqueous and organic phase).

Systems containing 2-propanol beads are the only cases where LLE diagrams cannot be reproduced correctly. For the water/benzene/2-propanol system, LLE diagrams calculated using the four parameterization approaches are very similar although interaction parameters are different. Some deviations on the calculated compositional data in comparison to experiments are observable on solubility of solute in the aqueous and organic phases and also on the miscibility between solvents. For the water/hexane/2-propanol system, solute concentrations in aqueous and organic phases calculated with CG-MC simulations are in good agreement with experimental data but the miscibility between the aqueous and organic phases is not well reproduced. Moreover, it should be noted that LLE diagrams for the water/hexane/2-propanol system parameterized using approach #1 strongly deviates from experimental data with the most concentrated compositions in 2-propanol (i.e., 0.1 in the aqueous and organic phase). Our DPD simulations carried out with an explicit interface (section 3.c. and 3.d.) have

shown that these deviations are due to a phase separation of the solute and the organic phase.

Overall, the four parameterization approaches allow a good reproduction of LLE diagrams in comparison with experimental data as shown by the analysis of regression functions in Table 3. In most cases, approach #2 allows a better reproduction of the LLE diagram compared to approach #1. Approaches #3 and #4 which take into account the isothermal compressibility of each fluid may, in some cases, slightly improve miscibility between solvents.

Another point to be highlighted is that parameterization approaches #1 and #2 can provide different solute/solvent interaction parameters. For instance, the solute/solvent interaction parameters, / and / , obtained with the approach #1 are lower than those obtained with the approach #2 and are predicted to be very close to 25, the value of the interaction between like beads. But overall, compositional data are well reproduced in both cases for water/benzene/1,4-dioxane system. Therefore, there are several sets of parameters that allow to reproduce compositional data with CG-MC simulations. Prediction of interfacial phenomena/properties with DPD simulations thus represents an additional criterion for identifying the best parameterization approach(es) (see section 3.d).

c. Interface compositions

Interface compositions are difficult to characterize experimentally and are generally left unknown. Indeed, interfaces are very thin surfaces at the macroscopic scale, but at the mesoscopic scale, dimensions of interfaces are large, and consequently, construction of the simulations boxes with an explicit interface requires to know the composition of this region.

Therefore, CG-MC simulations in the (µ solute N solvent P zz T) osmotic ensemble have been conducted in order to impose the precise bulk phase composition and to predict the composition at the interface. This step allows a perfect comparison between the interfacial tension values predicted by DPD calculations with experimental data at exactly the same phase compositions. Chemical potential used in the CG-MC simulations and compositions resulting from simulations are given in the supporting information.

To determine the interface composition from numerical simulations, the thickness of the interfacial region should be defined. Following the method used by Rezaei et al. [START_REF] Rezaei | Self-Accumulation of Uncharged Polyaromatic Surfactants at Crude Oil-Water Interface: A Mesoscopic DPD Study[END_REF] , the thickness of the interface is determined with the criterion "90-90" which defines the distance between two positions where the densities of two phases are 90% of their own bulk densities. Interfacial concentration and thickness of the interfaces were extracted from CG-MC simulations and are plotted as a function of the solute bulk concentration, here, the molar fraction of solute in the aqueous phase. Results for the water/benzene/1,4-dioxane system are reported in Figure 5 and data for other systems are given in the supporting information. predictions and were not used to derive any interaction parameters. For all systems, variation of IFT are presented as a function of the bulk solute concentration in the aqueous and in the organic phases and compared together with experimental data [START_REF] Backes | Interfacial tensions in binary and ternary liquid-liquid systems[END_REF][START_REF] Paul | Interfacial tensions in ternary liquid-liquid systems[END_REF][START_REF] Sada | Interfacial tensions of two-phase ternary systems[END_REF] . For water/benzene/1,4-dioxane, water/chloroform/acetone and water/benzene/acetic acid systems, IFT values are also compared with predictive calculations using the COSMO-RS method implicit solvent model combined with Density Functional Theory (DFT) calculations, extracted from the work by Andersson et al. [START_REF] Andersson | First-Principles Prediction of Liquid/Liquid Interfacial Tension[END_REF] . They showed that their methodology is fast, reliable and requires no experimental input data. Difference between calculated and experimental IFT values is given by the relative deviation, % , following the equation ( 28). % 100.

In addition, criteria, , were employed to compare calculated variation of IFT values with experimental data in order to determine the most suitable parameterization approach (Table 3).

In a first step, results on IFT are analyzed on systems that do not include intramolecular forces (Figure 6 to Figure 9). The IFT values of water/benzene and the water/chloroform interface in absence of solute are predicted in a good agreement with experimental data. For instance, Figure 6 shows IFT value for the water/benzene/1,4-dioxane system and the relative deviation is lower than 1.8 % for approaches #1 and #2. The approach #3 underestimates the IFT ( 10.3% ) while the approach #4 overestimates the IFT ( 18.2%). For the water/chloroform/acetone system presented in Figure 7, water/chloroform IFT in absence of solute is slightly overestimated with approaches #1 and #2 with a relative deviation of 10.4%. The approach #3 underestimates the IFT ( 7.6%) while the approach #4 overestimates the IFT ( 18.8%).

In presence of solute, prediction of variation of IFT presented in Figure 6 for the water/benzene/1,4-dioxane system are in very good agreement with experimental data using the parameterization approach #1 with an 0.997. Other parameterization methods underestimate the variation of IFT in comparison to the experimental data. For example, with the composition dependence of interaction parameters that is not represented by using constant interaction parameters over the range of compositions studied. A better reproduction of hydrogen bonds or electrostatics interactions that are not sufficiently well included in the standard DPD model (including neither attractive nor electrostatic terms) could improve the results. The work presented by Kacar and de With 71 seems to be a promising way to take into account hydrogen bonding within alcohol-water mixtures. They proposed to modify the conventional DPD potential by adding a Morse potential term to represent hydrogen bonding interaction.

In the coarse-grained model, hexane molecules are represented by two beads bonded using an intramolecular force. Results on IFT for these two systems are presented in Figure 10 and Predicted IFT values of water/hexane in absence of solute are not in a good agreement with experimental data and are largely underestimated. For instance, Figure 10 shows IFT values for the water/hexane/acetone system and the relative deviation is about 20.9% for approaches #1 and #2, and for the transferability approach. The approach #3 underestimates the IFT by 31.1% in absence of solute. The best result is obtained with the approach #4 with a relative deviation about of 8.3%. These deviations are commonly found in literature for water/alkane systems [START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF][START_REF] Goel | Estimation of interfacial tension for immiscible and partially miscible liquid systems by Dissipative Particle Dynamics[END_REF] .

In presence of solute, variation of IFT for water/hexane/acetone and water/hexane/2propanol are not well reproduced as shown in Figure 10 and Figure 11, respectively. For water/hexane/2-propanol system parameterized using approach #2, prediction of the variation of IFT is in a relative good agreement with experimental data but on a limited range of compositions. A simulated phase separation occurs when the molar fraction of 2-propanol exceeds 0.1 in the aqueous and in the organic phase. Although, parameterization approaches that needs compositional data for only one LLE and Hildebrand solubility parameter, while the second (#2) is based on several compositional data of LLE. These approaches are based on the main hypothesis that bulk phases in LLE are totally immiscible. Influence of interaction parameters between like DPD particles on IFT and miscibility has also been studied in the third approach (#3) using the isothermal compressibility ( ) of each fluid and in a fourth approach (#4) by adding a term corresponding to the degree of coarse-graining for water beads. It is important to remark that the parameterization methods proposed in this work make use of phase composition (or solubility) as input data, independently of the source type (either experimentally or predicted by atomistic molecular simulations, equations of states (EoS) [START_REF] Peng | A New Two-Constant Equation of State[END_REF][START_REF] Renon | Local compositions in thermodynamic excess functions for liquid mixtures[END_REF] or correlative methods such as partition coefficient from Quantitative Structure-Activity Relationship (QSPR methods) [START_REF] Bodor | A new method for the estimation of partition coefficient[END_REF][START_REF] Abraham | Hydrogen Bonding. 32. An Analysis of Water-Octanol and Water-Alkane Partitioning and the Δlog P Parameter of Seiler[END_REF][START_REF] Tetko | Prediction of n -Octanol/Water Partition Coefficients from PHYSPROP Database Using Artificial Neural Networks and E-State Indices[END_REF] .

Concerning the liquid-liquid equilibrium, GC-MC simulations in the Gibbs (NVT) ensemble were carried out in order to compute the liquid-liquid equilibrium. This method allows a precise determination of the phase diagram of the ternary systems analyzed in this work (water/1,4dioxane/benzene, water/acetone/chloroform and water/benzene/acid acetic, water/benzene/2propanol, water/hexane/acetone and water/hexane/2-propanol). When comparing the different methods of parameterization proposed in this work, we can conclude that in average, approach #2 seem to be the most accurate one to reproduced LLE. In addition, the use of different parameters between like beads (a ii a jj ) slightly improve the miscibility between bulk phases.

One of the important contribution of this work is the use of the CG-MC simulations in the osmotic (µ solute N solvent P zz T) ensemble for the preparation of systems with an explicit interface.

This step allows to work at the same bulk concentrations as the experimental data by computing the precise bulk phase compositions as well as the prediction of the interface composition. It is

  (DPD) has been introduced for the first time in 1992 by

Groot and Warren 43 ,

 43 all fluids have the same isothermal compressibility. Authors obtained quantitative results in agreement with experimental data for the critical micelle concentration and aggregation number for several typical surfactants of different chemical structures. Alternatively,

Figure 1 .

 1 Figure 1. Representation of a ternary system where bulk phases are immiscible

a

  Direct fitting means numerical optimization by iterative reduction of the deviation with respect to experimental data. d. Computational simulations details In this work, three different statistical ensembles are used depending on phenomena under investigation (a workflow summarizing the proposed methodology is available in Figure S1 of supporting information):-DPD model is combined with coarse-grained Monte Carlo (CG-MC) technique in order to simulate systems in the Gibbs (NVT) ensemble. This approach has already been used and validated by Wijmans et al.[START_REF] Wijmans | Phase behavior of monomeric mixtures and polymer solutions with soft interaction potentials[END_REF] with beads and soft potentials from DPD models. In the Gibbs ensemble, two separated simulation boxes that can exchange particles are used with a constant total volume . Thus, it is possible to describe equilibrium between two phases without considering explicitly the interface. Gibbs (NVT) ensemble simulations are used in our work in order to compute phase diagram and, thus, to check the relevance of parameterization methods to reproduce composition in bulk phases. Three different types of Monte Carlo moves are used:(1) translation of beads, (2) transfer of beads between the two boxes and (3) concerted change of volume of each box. These movements are described in greater details in the work of Wijmans et al.[START_REF] Wijmans | Phase behavior of monomeric mixtures and polymer solutions with soft interaction potentials[END_REF] . In addition, for hexane molecules which are represented by two beads, rigid body rotation and configurational regrowth moves are added. During CG-MC simulations in the Gibbs (NVT)

  to compute the number of surfactants necessary to reach an imposed IFT value between water and oil phases. CG-MC simulations in osmotic (µ solute N solvent P zz T) ensemble with an explicit interface allow to predict the solute concentration at the interface from known bulk compositions. The imposed chemical potentials are obtained from previous simulations in the Gibbs (NVT) ensemble. Three different types of Monte Carlo moves are used: (1) translation of beads, (2) change of volume along z-axis (perpendicular to the interface) and (3) insertion or removal of solute beads. For hexane molecules, rigid body rotation and configurational regrowth moves are added.

Figure 2 .

 2 Figure 2. Ternary systems studied with their respective coarse-grained representations. The degree of coarse graining , the mean atomic volume of beads and the characteristic length of each coarse-grain representation are given in the last column.

Figure 3

 3 Figure 3 (a) shows dependence of

Figure 3

 3 (a) shows that using approach #1, large variations of the / parameter can be observed as a function of the composition used and the values obtained vary from -3.84 to -14.65 (i.e. / varies from 11.56 to -26.29, negative values cannot be used in the DPD model used in this work). Large variations of the / parameter are also observed for the water/benzene/2propanol and water/hexane/2-propanol systems (see Figure S5 and Figure S7 in supporting information).

Figure 3

 3 Figure 3 (b) presents the compositional dependence of

Figure 4 .

 4 Figure 4. Liquid-liquid equilibrium ternary diagrams for the water/1,4-dioxane/benzene system at 298.15K for (a) approach #1, (b) approach #2, (3) approach #3 and (4) approach #4.

Figure 5 .

 5 Figure 5. (a) Variation of the 1,4-dioxane interfacial concentration as a function of the molar fraction of 1,4-dioxane in the aqueous phase. (b) Variation of the water/benzene interface thickness as a function of the molar fraction of 1,4-dioxane in the aqueous phase.

Figure 5 (

 5 Figure 5 (a) indicates that the interfacial concentration ( ) increases with solute molar fraction in bulk phase.Figure 5 (b) shows that the thickness of the interface increases by following

Figure 6 .

 6 Figure 6. Experimental and predicted variation of water/benzene IFT as a function of the 1,4dioxane concentration in (a) the aqueous phase and (b) the organic phase. Uncertainty on the value of the interfacial tension in DPD simulations is about 0.20 mN/m.

Figure 7 .

 7 Figure 7. Experimental and predicted variation of water/chloroform IFT as a function of acetone concentration in (a) the aqueous phase and (b) the organic phase. Uncertainty on the value of the interfacial tension in DPD simulations is about 0.20 mN/m.

Figure 8 .

 8 Figure 8. Experimental and predicted variation of water/benzene IFT as a function of acetic acid concentration in (a) the aqueous phase and (b) the organic phase. Uncertainty on the value of the interfacial tension in DPD simulations is about 0.20 mN/m.

Figure 9 .

 9 Figure 9. Experimental and predicted variation of water/benzene IFT as a function of 2-propanol concentration in (a) the aqueous phase and (b) the organic phase. Uncertainty on the value of the interfacial tension in DPD simulations is about 0.20 mN/m.

  Figure 11.

Figure 10 .

 10 Figure 10. Experimental and predicted variation of water/hexane IFT as a function of acetone concentration in (a) the aqueous phase and (b) the organic phase. Uncertainty on the value of the interfacial tension in DPD simulations is about 0.20 mN/m.

Figure 11 .

 11 Figure 11. Experimental and predicted variation of water/hexane IFT as a function of 2-propanol concentration in (a) the aqueous phase and (b) the organic phase. Uncertainty on the value of the interfacial tension in DPD simulations is about 0.20 mN/m.

  ). Following the same methodology, Rezaei and Modarress[START_REF] Rezaei | Dissipative particle dynamics (DPD) study of hydrocarbonwater interfacial tension (IFT)[END_REF] have shown that DPD simulations can be used to study the temperature dependence of IFT of water/hydrocarbon systems in the range from 298.15 to 348.15 K at 0.1 MPa. More complex

interfacial systems such as water/surfactant/hydrocarbon systems have already been modeled using DPD simulations. For example, Rekvig et al.

  . Depending on the coarse grained level, a bead of water represents 4 or 5 water molecules. The volume of a bead is fixed as the arithmetic mean of bead volumes as proposed by Rezaei and Modarress 27 . Following Maiti and McGrother 29 , the cut-off radius is given from the volume of the DPD particles with ̅

	water/benzene/acetic	acid,	water/benzene/2-propanol,	water/hexane/acetone	and
	water/hexane/2-propanol. For each of these systems, experimental data of bulk phase
	compositions and corresponding interfacial tension values are available in the literature 65-67 and
	are reported in the supporting information. Compositional data are molar fractions or mass
	fractions of each molecule in bulk phases of liquid-liquid equilibria and are converted to volume
	fractions according to the coarse-grained representation (see Figure 2 and supporting
	information).

/

. The overall DPD density ̅ is 3, which is a value commonly used in DPD simulations.

Above mentioned methodologies for the parameterization of DPD interactions were applied on six different ternary systems (Figure

2

): water/benzene/1,4-dioxane, water/chloroform/acetone,

Table 2 .

 2 Properties of individual components (DIPPR) at 298.15 K and 1 bar.

	Molecules	(Å 3 )	ΔH vap (kJ/mol)	δ (J/cm 3 ) 1/2
	Water	30.07	43.982	47.9
	Benzene	148.58	33.871	18.7
	Chloroform	133.68	31.393	19.0
	Hexane	218.13	31.549	14.9
	1,4-dioxane	142.25	38.595	20.5
	Acetone	122.61	31.166	19.7
	Acetic acid	95.70	24.313	19.5
	2-propanol	127.64	46.081	23.8

Table 3 .

 3 Values of regression functions for the four parameterization approaches to reproduce the experimental data. Phase compositions are used for parameterization, IFT correspond to predictions.

	Number of compositional	Solvent 1	Solvent 2	Solute	Parameterization approaches	[a]	[b]
	data available,
					approach #1	0.994	0.997
	5	water	benzene	1,4-dioxane	approach #2 approach #3	0.994 0.985	0.847 0.899
					approach #4	0.988	0.909
					approach #1	0.971	0.970
	5	water	chloroform	acetone	approach #2 approach #3	0.983 0.995	0.964 0.871
					approach #4	0.989	0.924
					approach #1	0.875	0.677
	5	water	benzene	acetic acid	approach #2 approach #3	0.967 0.759	0.893 0.908
					approach #4	0.768	0.911
					approach #1	0.900	0.805
	7	water	benzene	2-propanol	approach #2 approach #3	0.888 0.938	0.237 0.814
					approach #4	0.954	0.537
					Transferability [c]	0.758	0.905
	11	water	hexane (2 beads)	acetone	approaches [#1-#2] approach #3	0.998 0.997	0.849 0.812
					approach #4	0.998	0.851
					approach #1	0.514	0.818 [d]
	11	water	hexane (2 beads)	2-propanol	approach #2 approach #3	0.988 0.996	0.416 0.557
					approach #4	0.997	0.377

RESULTS AND DISCUSSION

a. Composition dependence of the Flory-Huggins interaction parameters A study of the compositional dependence of / Flory-Huggins parameters determined from experimental compositions was conducted. In this work, / are determined using one (approach #1) or several (approach #2) compositional data of bulk phases of liquid-liquid equilibria. Dependence of / parameters as a function of composition data obtained using approaches #1 and #2 for the water/benzene/acetic acid system is presented in the Figure 3. For other systems studied in this work, the dependence of the / parameters according to approaches #1 and #2 are given in the supporting information. Compositions of bulk phases for liquid-liquid equilibrium are computed using Monte Carlo coarse-grained simulations in the Gibbs ensemble (NVT). The four parameterization approaches are compared in details using water/benzene/1,4-dioxane system, and results for other systems are given in the supporting information. For the water/benzene/1,4-dioxane system, interaction parameters obtained with the approaches #1, #3 and #4 were determined from the less concentrated solute composition. With the approach #2, the four compositions with the lowest concentrations of 1,4-dioxane were selected ; the fifth composition is not taken into account because of the excessive solubility of water in benzene (see supporting information S2 for details). Phase compositions are visualized using ternary diagrams and corresponding interaction parameters used to compute LLE are indicated below each diagram (Figure 4). In addition, the coefficient of regression , is used to quantitfy the accuracy of our parameterization approaches in the reproduction of the solubility of solutes in the aqueous and organic phases (Table 3).

approach #2, predicted IFT values are about two times lower than the experimental data ( 0.847 . Deviations on IFT values obtained with our DPD simulations are of the same order of magnitude as those obtained with other predictive methods such as COSMO-RS (Figure 6 (a)

and (b)).

For the water/chloroform/acetone system (Figure 7), all the parameterization approaches provide a trend of the IFT variation with acetone molar fraction in agreement with experimental data. The best results are obtained with the approaches #1 ( 0.970) and #2 ( 0.964). In addition, DPD simulations predict more precisely the variation of IFT as a function of the organic phase composition than the COSMO-RS method as shown in Figure 7 (b).

In the case of the water/benzene/acetic acid system (Figure 8), IFT are overestimated by about 30% with approaches #2, #3 and #4 ( varies from 0.893 to 0.911) and by about 40% with the approach #1 ( 0.677). Note that prediction of IFT values can be compared to the experimental data only if compositions are reproduced in bulk phases. Although, approach #2 does not provide the best prediction of IFT compared to approaches #3 and #4, bulk phases compositions are better reproduced. DPD simulations provide better predictions of variation of IFT than the COSMO-RS method. In Figure 8 Our parameterization approaches do not seem adapted to the water/benzene/2-propanol system (Figure 9). Indeed, values and variations of IFT are not in agreement with experimental data using the four parameterization approaches. In addition, we have showed in section 3.b. that LLE diagram is not precisely reproduced. We can assume that these deviations can be attributed to we propose on the basis of compositional data allow to reproduce solubility of acetone or 2propanol in the water/hexane liquid-liquid equilibrium, the variation of the IFT is not reproduced. We can assume that prediction of interfacial tension for systems containing hexane molecules seems to require a better parameterization of the intramolecular forces [START_REF] Lee | Parametrization of Chain Molecules in Dissipative Particle Dynamics[END_REF] 23) that derives from the Flory-Huggins theory cannot be exploited, therefore, parameters must be calibrated directly on the experimental compositional data.

CONCLUSIONS

A consistent methodology for parameterization of interactions taking into account the miscibility of species in a diphasic solution and to predict quantitatively the interfacial tension is proposed in this work. Two parameterization approaches based on the introduction of compositional data into equations of the Flory-Huggins theory were applied. The first one (#1) important to remark that the interface composition in equilibrium with the bulk phases is an important property which is not easily accessible experimentally. This is a crucial preliminary step to compute the interfacial tension (by any molecular simulation method).

IFT for the ternary systems were calculated using DPD simulations. Our results show that the proposed parameterization approaches allow us to reproduce the trend of interfacial tension variation as a function of the solute molar fraction in the aqueous and organic phases for the water/benzene/1,4-dioxane, water/chloroform/acetone and water/benzene/acetic acid systems.

Our results can be considered as quantitative since we obtain good agreements between simulated and experimental data ( 0.893). We notice that results on interfacial tension for systems containing intramolecular forces are not accurate enough when compared with the experimental data. The accuracy of the prediction can be improved by an optimization of the intermolecular forces (which was out of the scope of the present work).

Among the parameterization approaches proposed in this work, the approach #2 seems to be the best compromise for reproducing solute solubility in bulk phases and for quantitatively predicting interfacial tension.

Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENT

Authors would like to thank Dr. R. Lugo for fruitful discussions.