Calculation of the Saturation Properties of a Model ... - ACS Publications

Subscriber access provided by Kaohsiung Medical University

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Calculation of the Saturation Properties of a Model Octane-Water System Using Monte Carlo Simulation Wenjing Guo, Prannay Bali, and Jeffrey R Errington J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b01411 • Publication Date (Web): 15 May 2018 Downloaded from http://pubs.acs.org on May 17, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Calculation of the Saturation Properties of a Model OctaneWater System Using Monte Carlo Simulation

by

Wenjing Guo, Prannay Bali and Jeffrey R. Errington*

Department of Chemical and Biological Engineering, University at Buffalo Buffalo, New York 14260-4200

Submitted to: J. Phys. Chem. B

*

Corresponding author, electronic address: [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We use Monte Carlo simulation to compute the saturation properties of a model octane-water system. The system phase separates into water-rich liquid and vapor phases, octane-rich liquid and vapor phases, and water-rich liquid and octane-rich liquid phases at various conditions. We outline a strategy for determining the saturation properties of the mixture over a wide range of temperatures, pressures, and compositions. The approach begins with direct calculations that enable one to locate a single saturation point. A variety of expanded ensemble schemes are then used to trace saturation curves along paths of interest. We show how the overall strategy provides a means to construct pressure-composition diagrams at fixed temperature and temperature-composition diagrams at fixed pressure. In addition, we demonstrate how the approach is used to trace the liquid-liquid-vapor triple line over a wide range of temperatures. Simulation data are compared with experimental data, when available. Overall, our results show that the approach provides an efficient means to calculate the saturation properties of a binary system.

2 ACS Paragon Plus Environment

Page 2 of 61

Page 3 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Water-hydrocarbon mixtures are of considerable interest due to their prominence in many industrial applications, including those related to oil recovery, biological fuels, and environmental remediation.1-4 These mixtures exhibit interesting phase behavior, with the system phase separating into a water-rich liquid and vapor at relatively low pressure and high water concentration, a hydrocarbon-rich liquid and vapor at relatively low pressure and high hydrocarbon concentration, and a water-rich liquid and hydrocarbon-rich liquid at relatively high pressure. A comprehensive understanding of the phase behavior of these mixtures is important for the efficient design of current and emerging applications. For example, in steam flooding,4 steam is injected into reservoirs to displace oil. In this process, heat and mass transfer cause multiphase separation, with the resulting phase distribution linked to the underlying phase behavior. Apart from industrial interest, studying the phase behavior of water-hydrocarbon mixtures provides insight regarding the manner in which molecular-level interactions between polar and non-polar molecules manifests in bulk phase separation.5 Several experimental studies6-12 have been reported regarding the phase behavior of aqueous systems. conditions.

These reports provide saturation data over a relatively wide range of

However, experimental data are often scarce for aqueous mixtures containing

complex molecules and/or for mixtures at high temperature and/or pressure. Equations of state (EOS) are commonly used in engineering practice to describe the phase behavior of light hydrocarbon mixtures because of their computational efficiency. In recent years, considerable effort has been devoted to the development of EOS13-18 that directly account for hydrogen bonding, such as the cubic-plus-association (CPA)13 and the statistical associating fluid theory



(SAFT).18 For example, Papaioannou et al. have used the group-contribution statistical

associating fluid theory (SAFT- ) approach to predict the vapor-liquid and liquid-liquid equilibria for binary water + n-alkane mixtures.19 Studies of n-hexane, n-octane and n-decane

show good agreement with experimental data, but application to heavier hydrocarbons is less successful. This result is a common feature of most equations of state, and can be attributed to the low vapor pressure exhibited by these compounds.20 Monte Carlo (MC) and Molecular dynamics (MD) simulations have proven successful in understanding the phase behavior of water-hydrocarbon mixtures from a molecular level perspective.17 A number of molecular simulation studies have been conducted to investigate the

mutual solubilities of water + -alkane mixtures.5, 21-23 The description of the solubility of alkane in water is quite challenging due to the nonideal behavior of the mixture. The solubility of hydrocarbon in water is several orders of magnitude lower than the corresponding solubility of water in a hydrocarbon-rich liquid. Ferguson et al.22 used the replica-exchange moleculardynamics (REMD) technique together with an incremental Widom insertion scheme to calculate the solubility of -alkanes (ranging from ethane to -docosane) in water under conditions of vapor-liquid equilibrium. Similarly, Ballal et al5 predicted the solubility of water in different

alkanes (methane to n-dodecane) using Monte Carlo simulation. These molecular simulations mainly focused on liquid-vapor equilibrium. Relatively few studies have focused on liquidliquid phase coexistence. In this work, we use a Monte Carlo simulation based approach to characterize the fluid phase behavior of a model water-octane mixture. We employ a series of methods introduced previously by our group24-25 to trace water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid saturation curves over a wide range of temperatures, pressures, and compositions. 4 ACS Paragon Plus Environment

Page 4 of 61

Page 5 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Within the approach, we first complete direct grand canonical (GC) simulations at relatively high temperature to obtain the density or composition probability distribution for a system over a range of conditions that spans the two coexisting phases. Histogram reweighting26 is then used to locate the activities (chemical potentials) that satisfy coexistence conditions. In the next step, we trace a given saturation curve over a relatively wide range of temperatures, pressures, and/or compositions. To determine coexistence conditions, we first formulate an estimate for how relevant field variables vary along a saturation curve using results generated from direct GC simulations. Expanded ensemble simulations are then used to compute the free energies of the coexisting phases over this range of states. We then again leverage histogram reweighting techniques to locate the true saturation conditions at multiple points along the path. Data generated from this approach are used to construct pressure-composition and temperaturecomposition phase diagrams. Within these simulations, it becomes difficult to sample the liquid state at low temperatures due to the moderate size of octane and relatively high densities of the liquids.

To overcome these difficulties, we employ strategies that include reservoir grand

canonical Monte Carlo,27 a growth expanded ensemble scheme,28 and hybrid MC methods29 to facilitate molecular insertions and deletions and sampling of intramolecular degrees of freedom. This work is part of a broader study in which we are interested in the interfacial properties of the water-octane mixture in the vicinity of a mineral substrate. Within the methods we employ to understand the behavior of the fluids at interfaces, properties are typically evaluated at bulk saturation conditions. Studying the bulk phase behavior of the octane-water mixture is the important first step toward understanding the interfacial behavior of the mixture. This paper is organized as follows. In the next section we describe the simulation methods that we use to compute bulk phase coexistence properties. We then discuss the molecular 5 ACS Paragon Plus Environment


models employed and provide details for the simulations performed in this work. We next present our simulation results and draw comparisons with literature data when possible. Finally, we provide concluding remarks.

Simulation Methods Overview. A multi-step approach is used to obtain saturation properties over a wide range of state conditions. In the first step, we identify a relatively small number of coexistence points along a given saturation line (e.g., water-rich liquid and vapor). The method we employ to directly locate specific phase coexistence points is outlined in the Direct GC Simulations section. We then use a series of expanded ensemble methods to trace saturation lines over a wide range of state conditions. These methods require a known saturation point (e.g., one from a direct calculation). The general approach is outlined in the Expanded Ensemble section, and we then detail specific variants. We first describe how to determine the composition dependence of saturation properties at a given temperature by employing activity fraction expanded ensembles (Activity Fraction Expanded Ensemble section). We then outline strategies for determining how coexistence curves evolve with temperature at either fixed activity fraction or fixed pressure (Temperature Expanded Ensemble section). Finally, we describe an approach one can use to trace the triple line of a binary mixture over a wide range of conditions (Triple Line section). Direct GC Simulations. In this work, we focus on a binary mixture in a grand canonical

ensemble in which the temperature , chemical potentials of the two species , , and volume

are fixed and the particle numbers of the species , and energy fluctuate. In what

follows, we use the inverse temperature = 1/ ( is Boltzmann’s constant) and activities 6 ACS Paragon Plus Environment

Page 6 of 61

Page 7 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


= ( represents the component of the molecular partition function of species stemming from integration over momenta) to describe the temperature and chemical potentials, respectively. For convenience, the activities are often expressed in terms of an activity fraction

= /( + ) and activity sum = + . The activity fraction is used to drive changes in the composition of the system.

As the activity fraction shifts from zero to unity, the

composition of the system changes from pure component 1 to pure component 2. To locate a coexistence point, we collect the probability distribution Π( ) at fixed

, , , over a range of that includes both coexisting phases. The normalized probability

( ) is related to the partition function Υ( , , , ) and grand partition function distribution Π Ξ( , , , ) as follows,

( ; ) = Π

$ Υ( , , , ) Ξ( , , , ) #

( 1)

At near-coexistence conditions, Π( ) is bimodal. The weight (area) associated with a given peak is proportional to the grand partition function, and therefore the weight provides the

pressure of the relevant phase % = ln Ξ to within an additive constant. The saturation point is located by identifying the activity that results in equal peak weights for the two phases. The relevant histogram reweighting expression is

ln Π( ; ) = ln Π( ; ( ) + ln ) ( *

where ( is the activity at which the distribution Π( ; ( ) is obtained via simulation.

( 2)

Saturation properties of interest are readily obtained from the Π( ) probability

distribution. For example, the composition + is given by, 7



+ =

〈 〉 〈 〉 + 〈 〉

Page 8 of 61

( 3)

where,

〈 〉 =

∑#$ Π( ) ∑#$ Π( )

( 4)

The absolute pressure at saturation requires knowledge of the thermodynamic properties of pure component 2. Evaluation of Eq. (1) at = 0 gives,

(0; ) = ln Υ(0, , , ) − ln Ξ( , , , ) ln Π

( 5)

where Υ(0, , , ) is the grand partition function of the pure species 2 system, which is related

to the pressure of the single component system %1234 ( , , ) = ln Υ(0, , , ) .

The

pressure of the mixture %( , , , ) = ln Ξ( , , , ) is then obtained from

(0; ) %( , , , ) = %1234 ( , , ) − ln Π

( 6)

In all cases, properties are evaluated at the saturation conditions for the mixture. Within mixture

GC simulations, we regularly sample the = 0 state. Pure component properties are readily available from simulation data for the single component system. We find that the approach outlined above works well for locating phase coexistence points along the three saturation curves of interest in this work. For water-rich liquid-vapor saturation, we take species 1 to be water, whereas for octane-rich liquid-vapor saturation, we take species 1 to be octane. For liquid-liquid saturation, we find it convenient to take species 1 to be water.


Page 9 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Expanded Ensemble. We now focus on a strategy for determining saturation properties over a wide range of temperatures, pressures, and/or compositions. Direct calculations described above enable us to compute the saturation properties at a specific state point. However, it is computationally expensive to conduct direct simulations at each state point of interest. In previous studies,24-25, 30-34 our group introduced a means for efficiently tracing the evolution of a saturation curve over a path of interest (e.g. temperature, pressure, activity fraction). In what follows, we first describe the general approach and then provide specific examples. We begin by establishing a series of subensembles along the path of interest.

Each

subensemble is defined by a set of field variables 5 = ( , , ). In practice, the subensembles are typically separated by a fixed difference in one of the field variables. The goal is to establish the conditions 5678 that satisfy the conditions for phase equilibrium within each subensemble. To

initiate the calculation, we use coexistence data obtained from direct calculations to formulate an estimate 59 for how 5 varies along the saturation line. We then perform independent expanded

ensemble (EE) simulations to obtain the relative grand potential Ω; for each phase along the path. The absolute grand potential for a given phase is constructed by linking the relevant EE

probability distribution Π ∑C; ( ; 59 ) 9

>

( 11) ( 12) ( 13)

where C; (K; 59 ) provides the subensemble probability distribution for variable K. We also note

that for cases in which 59 and 5@ differ significantly and the reweighted histograms fall outside of the range of states sampled with 59 , we find it convenient to use the approximate expressions,

Δ(Ω; )( , , 9 → @ ) = L M( @ − 9 )

ΔΩ; E → @ , , F = −L MEln @ − ln F 9

9

ΔΩ; E , → @ , F = −L MEln @ − ln F 9

9

( 14) ( 15) ( 16)

When representing ( , ) as (, ), the following free energy difference expression can be useful,

ΔΩ

; (,

∑C; ( ; 59 )( @ / 9 )# → , ) = −ln = > ∑C; ( ; 59 ) 9

@

The approximate version is,


( 17)


ΔΩ; (, 9 → @ , ) = −L M(ln @ − ln 9 )

Page 12 of 61

( 18)

Before closing this section, we address a detail related to the reference point. In Eq. (7) above, we took 53 as the exact conditions that are sampled in one of the subensembles. In some cases,

we do not know the absolute free energy Ω; (53 ) at the specific point 53 , but at a relatively close

point 5N . We then use Eq. (7) with,

Ω; (53 ) = Ω; (5N ) − ΔΩ; (53 → 5N )

( 19)

to obtain the free energy at 53 . The probability distributions collected in the subensemble

characterized by 53 are used in the free energy difference expressions above. The reference

grand potential is typically obtained from a direct simulation described above with Ω; (5N ) = −%; (5N ).

Activity Fraction Expanded Ensemble. Activity fraction expanded ensemble (AFEE) simulations are used to determine how saturation properties evolve at fixed temperature. In

AFEE simulations, 5 is characterized by (, ln ) and we aim to identify the (, ln ) that satisfy

phase equilibrium criteria at a specified with ∈ [0,1]. For the water (1)-octane (2) mixture, there are generally three relevant coexistence curves: water-rich liquid-vapor, octane-rich liquidvapor and liquid-liquid coexisting regions. The three coexistence curves meet at a triple point (8 , ln 8 ). We trace the three curves independently. Before initiating the AFEE simulations, it

is useful to have an estimate of the triple point conditions. Data from direct simulations are used to provide an initial estimate. For the liquid-vapor saturation curves, we construct a series of

subensembles separated by a fixed difference in . The range is split into two parts. For the 12 ACS Paragon Plus Environment

Page 13 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


water-rich liquid-vapor region, we focus on the range ∈ [0, 8 + Δ] and for the octane-rich

liquid-vapor regions, we focus on ∈ [8 − Δ, 1], where Δ represents an activity fraction span

over which we attempt to push into the metastable region. For the liquid-liquid saturation curve,

we construct a series of subensembles separated by a fixed difference in ln . We cover the

range ln ∈ [ln 8 − Δ ln , ln N7R ] that spans from the metastable to the stable liquid-liquid coexistence region. We generate an estimate of the saturated activities by assuming a linear relationship

between ln and along all saturation lines. AFEE simulations are then conducted, and the

general method outlined above is used to compute saturation conditions. In practice, we modify at constant , at constant , or at constant to satisfy Eq. (9). In other words, we solve

one equation in one unknown. The data generated from this approach enable one to construct traditional % − +, S diagrams at fixed temperature.

Temperature Expanded Ensemble. Temperature expanded ensemble (TEE) simulations are used to determine how saturation properties evolve with temperature. We discuss two TEE methods: (1) TEE at a fixed activity fraction and (2) TEE at a fixed pressure. For the first case, 5

is characterized by (, ln ) and we aim to deduce (, ln ) over a wide range of temperatures at

a specified . We construct a series of subensembles separated by a fixed difference in and

assume a linear relationship between and ln to formulate an estimate of the saturation curve. TEE simulations are conducted, and the general method outlined above is used to compute

saturation conditions. To solve Eq. (9), we modify at constant and . This approach is applied to each of the three coexistence curves of interest.



In the second method, we aim to trace saturation conditions at constant pressure. It follows

that both and vary with . We construct a series of subensembles separated by a fixed

difference in and assume linear relationships between and ln and between and to formulate an estimate of the saturation curve. TEE simulations are then conducted to obtain the relative grand potentials for the two phases of interest. The practical aspects of the solution differ from the cases outlined above.

The constant pressure constraint requires that we

simultaneously satisfy Eq. (9) and Ω;$ (5@ ) = −%( , where %( is the constant pressure. In

practice, we modify both and (or and ) at constant temperature and identify the conditions that satisfy the two constraints. The data generated from this approach enable one to construct traditional − +, S diagrams at fixed pressure.

In what follows below, we use this approach to trace saturation curves in the vicinity of the triple point. A strategy analogous to that outlined above for AFEE simulations near a triple point is employed. More specifically, we estimate the location of the triple point and then perform TEE simulations that are expected to continue slightly beyond this point. Triple Line. Within the discussion above, we outlined a means to identify a triple point at a given temperature. We now consider how to trace the triple line over a wide range of conditions. This line is constructed using another form of an expanded ensemble approach. The

goal is to determine how and vary with along the triple line. From AFEE simulations, we collect triple points at different temperatures. These data are used to generate an estimate for (, , ) along the triple line. TEE simulations are then conducted to obtain the relative grand

potential for the three phases: water-rich liquid, octane-rich liquid, and vapor phase. The triple point at each temperature requires that the grand potential of the three phases are equal. In other words, we need to satisfy a modified version of Eq.(9), Ω;$ (5@ ) = Ω;? (5@ ) = Ω;T (5@ ). In 14 ACS Paragon Plus Environment

Page 14 of 61

Page 15 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


practice, we modify both and (or and ) at fixed temperature and identify the conditions that satisfy these two constraints.

Model System and Simulation Details Molecular Models. As noted above, this work is part of a broader study in which we consider water-alkane mixtures at mineral surfaces. As such, we pursue a modeling scheme that is consistent with that employed by others35-37 to describe water-alkane-mineral systems. We work with the united-atom TraPPE model38 for octane and the SPC/E model39 for water. The

interaction parameters are provided in Table 1. For both models, the energy of interaction U(V)

between two interaction sites separated by four or more bonds and by a distance V is given by the potential,

U(V) = 4XY Z[

\Y \Y ^ 1 Y ] − [ ] _ `(V) + V V 4aX( V

( 20)

with,

1 h `(V) = b(Ve − V ) (2V + Ve − 3Ve )(Ve − Ve ) 0

V < Ve Ve ≤ V ≤ Ve V > Ve

( 21)

where \Y and XY are size and energy parameters, respectively, and denotes the value of the

partial charge placed at interaction site . Cross interaction parameters are computed using Lorentz-Berthelot combining rules.38 `(V) is a switching function that brings the Lennard-Jones

potential to zero between Ve = 9.2 Å and Ve = 10 Å.



Page 16 of 61

The water model is rigid. The intramolecular angle bending and torsion parameters for octane are taken from the TraPPE-UA force field.38 A harmonic potential is used to describe bond stretching within octane. The equilibrium bond length is taken from the TraPPE-UA forcefield38 and the spring constant is adopted from the OPLS force field.40 This modification was implemented to accommodate single-molecule hybrid MC moves24 used to sample the intramolecular degrees of freedom of octane. Electrostatic interactions are calculated using the Ewald sum method.41 Table 1. Potential parameters for the model studied

Model

Atom

\ (Å)

X (kJmolG )

()

O

3.166

0.650

-0.8476

H

0

0

0.4238

CHh

3.75

0.8149

0

3.95

0.3825

0

37

SPC/E

TraPPE-UA

36

CH

Simulation Details. In this work, all simulations are completed in a cubic box (32 × 32 ×

32Å) with periodic boundary conditions applied in all directions. We note that finite-size effects

with histogram-based methods and molecular models similar to those studied here were examined in an earlier study.24 For direct GC simulations, the MC move mix includes 35% molecular displacements and rotations, 5% single-molecule hybrid MC moves (octane only),24 and 60% molecular insertions (or growth) and deletions (or reductions). For expanded ensemble simulations, the MC move mix includes 53% molecular displacements and rotations, 5% single16 ACS Paragon Plus Environment

Page 17 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecule hybrid MC moves (octane only),24 40% molecular insertions (or growth) and deletions (or reductions), and 2% subensemble change moves. We first examine saturation properties of pure octane and water systems. For pure water,

direct GC simulations are performed at temperatures of = 450 K and 500 K, and the molecule

number spans from zero to approximately 1120 molecules. A TEE scheme is employed to trace saturation properties over a wide range of temperatures. An expanded ensemble is constructed

with a discretization of Δ = 2.0 × 10v J G . For the octane case, direct GC simulations are

performed at = 350 K and 450 K, and the molecule number spans from zero to approximately 110 molecules. Similar TEE simulations are conducted to trace the octane liquid-vapor saturation line. For the binary mixture system, direct simulations are first conducted at a relatively high

temperature ( = 450 K) and selected activity fractions = 0.0 and 0.2 (water-rich liquid-

vapor), = 0.8 and 1.0 (octane-rich liquid-vapor), and = 0.680 and 0.839 (liquid-liquid).

Expanded ensembles with a spacing of ∆ = 0.01 are used for activity fraction expanded ensemble simulations. Three AFEE simulations are performed independently to trace the waterrich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid saturation line. A variety of TEE strategies are applied to the mixture. In all cases, we employ a

subensemble spacing of Δ = 2.0 × 10v J G .

We trace saturation lines at fixed activity

fractions of = 0.2 (water-rich liquid-vapor), 0.937 (octane-rich liquid-vapor), and 0.89 (liquidliquid). A set of calculations are performed to construct temperature-composition diagrams at

pressures of % =0.5, 1.0, 5.0 bar. Finally, a set of simulations is completed to trace the triple line.



Page 18 of 61

Statistical uncertainties are determined by performing four independent sets of simulations. The standard deviation of the four simulation results is taken as estimate of the uncertainty.

Results and Discussion Pure Components. The saturation properties of the octane-water mixture depend sensitively on the molecular models selected for the two components. We begin by examining the ability of the molecular models to capture the saturation properties of the pure components. Figs. 1 and 2 provide the saturated vapor pressures and densities, respectively, for water and octane. These figures contain experimental data42 and simulation results from Economou43 (water), Martin and Siepmann38 (octane), and Rane et al.24 (water, octane).

The literature

simulation results were generated using the full potential model, and therefore account for longrange dispersion interactions. As noted above, this work is part of a broader project focused on interfacial properties of the water-octane system.

Given the challenges associated with

computing long-range dispersion interactions in inhomogeneous systems, we decided to truncate the Lennard-Jones interaction via a switching function. Our results for water are in close agreement with simulation data and reasonable agreement with experimental data.

The

comparison of simulation results suggests that the long-range contribution to the dispersion interaction has a negligible impact on the thermophysical properties of water.

These

observations are consistent with those reported in an earlier study with a different variant of a truncated water model.32


Page 19 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. The saturated vapor pressure against inverse temperature for water and n-octane. Open blue Circles represent the results from direct GC simulation for water and n-octane. Solid blue curves represent the results from temperature expanded ensemble simulations. experiment data.

42

Solid red lines represent data from Rane et al.

24

Solid black lines represent

for the full (non-truncated) models. 43

Open green circles denote simulation results (non-truncated) for water from Economou et al

and n-

38

octane from Martin and Siepman.

The results for octane point to a different outcome.

The saturated properties of the

truncated model are significantly different from those for the full potential. We first estimated the critical temperature z from the density scaling law {|}~ − {71 ~| − z | , with = 0.326. This approach provides a critical temperature for the truncated octane model of 481 K, which is 15% lower than the experimental value. In contrast, the full TraPPE model provides an estimate that is statistically equivalent to experiment. It follows that the saturated liquid (vapor) densities of the truncated model are consistently under(over)estimated. The full TraPPE model slightly overestimates the vapor pressure. This trend is amplified in the truncated model, with the simulation results approximately a factor of 3 to 10 higher than experimental data.

The

observations reported here are consistent with the findings of Lagüe et al.,44 who found that the 19 ACS Paragon Plus Environment


cutoff for the Lennard-Jones potential is more relevant for hexane than for water. It is clear that electrostatic interactions dominate in water and dispersion interactions dominate in octane.

Figure 2. The saturated density against inverse temperature for pure water and n-octane. Symbols are defined in the same manner as in Fig 1.

Octane-Water Mixture. We now explore the three coexistence regions of interest here: water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid. We begin by considering a


Page 20 of 61

Page 21 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


series of direct calculations and then examine the construction of various saturation lines via the expanded ensemble approaches outlined above. We first consider a direct calculation designed to capture a coexistence point along the water-rich liquid-vapor curve. Fig. 3 provides the water particle number probability distribution at the saturation conditions of = 450 K, ln = −10.8833 , and ln = − 9.465(1) ( = 0.200(6), and ln = −9.27(2)). The minimum and maximum water molecule number is set

such that we fully capture the vapor and liquid phases. The octane molecule number is not constrained. The area under the Π( ) probability distribution is proportional to the relative

pressure, and therefore represents the key metric for identifying the coexistence point. Histogram reweighting at constant ln is utilized to locate this point. To determine the total

pressure, we first use simulation data for pure octane to calculate a pressure of % = 0.119779 (8) MPa for the single-component vapor at ln = −10.8833 .

The total

pressure of the mixture % = 0.656(2) MPa then follows by summing the relative pressure of the mixture and that of the single component.



Page 22 of 61

Figure 3. The water particle number probability distribution for the octane-water system at the liquidvapor coexistence condition defined by = 450 K, η = 0.200(6) , and representative configurations are provided along the order parameter path.


ln = −9.27(2) . Three

Page 23 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 4 shows the water particle number probability distributions at = 450 K for two

points along the liquid-liquid saturation curve. The first point ( = 0.681(2) , ln =

−8.310(3)) lies just beyond the triple point and the second ( = 0.836(3), ln = −7.510(4)) is at a relatively high pressure. The minimum and maximum water molecule numbers are again set to capture both coexisting phases and the octane molecule number is not constrained.

Histogram reweighting at constant ln is utilized to locate the saturation point, and the relative coexistence pressure is extracted from the probability distribution.

In this case, we use

simulation data for pure octane to calculate pressures of % = 4.51(1) MPa and % = 22.75(2) MPa for the single-component liquid at the mixture saturation points.

The resulting total

Figure 4. The water particle number probability distribution for the octane-water system at the liquid-liquid

coexistence conditions defined by = 450 K, = 0.836 (3) , and ln = −7.510(4) (green line) and

= 450 K, = 0.681(2) , and ln = −8.310(3) (black line). Three representative configurations are

provided along the order parameter path.

saturation pressures for the mixture are % = 4.77(1) MPa and % = 22.94(2) MPa. We found 23 ACS Paragon Plus Environment


the water molecule number to be a convenient order parameter in this study. Note that the octane molecule number could also serve as a suitable order parameter. We did not explore the relative merits of these two options Activity Fraction Expanded Ensemble. In this section, we are interested in the composition dependence of saturation properties at constant temperature. The activity fraction expanded ensemble approach provides a convenient means to construct relevant phase diagrams at a specified temperature. To initiate the process, we construct estimates for the relationship

between and ln along the three saturation curves. The direct GC simulations at = 450 K described above provide saturation properties at = 0.0 , 0.195, 0.46 and 0.63 (water-rich

liquid-vapor), = 0.6, 0.8 0.93 and 1.0 (octane-rich liquid-vapor), and = 0.68 and 0.893 (liquid-liquid). Fig. 5 shows the approximate saturation curves that result from a polynomial fit

between and ln . These curves are used to generate an estimate of the triple point at 8 =

0.65 and ln 8 = −8.6. For the liquid-vapor coexistence regions, we split into two ranges:

∈ [0.0, 0.7] and ∈ [0.6, 1.0]. The activity fraction changes relatively little along the liquidliquid saturation curve, and therefore we focus on a total activity range of ln ∈ [−8.55, −7.00].

Fig.5 provides the final converged saturation curves. The inset focuses on the region near the triple point wherein the three saturation curves intersect. The true location of the triple point, 8 = 0.643(3) and ln 8 = −8.44(2), is close to the initial estimate. The most precise estimate of the triple point is obtained from the intersection of the two liquid-vapor saturation curves.


Page 24 of 61

Page 25 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Relationship between activity sum and activity fraction along three saturation lines at = 450 K. Circles represent data obtained from direct GC simulation. Solid green, blue, and red curves provide

the final solution obtained from AFEE simulation for the water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid saturation curves, respectively. Dashed lines represent the initial estimate of the saturation line. The inset contains data in the vicinity of the triple point.

Interestingly, the slopes associated with the water-rich liquid-vapor and liquid-liquid saturation curves are relatively close at the triple point. In general, the slope of a saturation curve in Fig. 5 is expressed as,

ln 1 = + ( − 1)(1 − ) ( − 1)

( 22)

where the derivative = ln ⁄ ln is obtained from a Clausius-Clapeyron-like relationship with,

ln Δ{ = =− ln Δ{


( 23)


− { and Δ{ = { − { , with and representing the two coexisting where Δ{ = {

phases. Fig. S1 in the Supplementary Information shows the relationship between ln and

ln at = 450 K. The dashed lines in the inset represent the linear fit of the simulation data

close to the triple point. We note that estimates for generated via Eq. (23) are consistent with

those slopes calculated in Fig. S1. At the triple point, is approximately 96, 0.02, and 11 for the

water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid coexistence curves, respectively. It follows that the first term in Eq. (22) dominates for the water-rich liquid-vapor and liquid-liquid saturation curves near the triple point, resulting in nearly common slopes of ln / ≈ 1/(1 − ). In contrast, the second term in Eq. (22) dominates for the octane-rich liquid-vapor saturation curve, with ln / ≈ −1/.

Fig. 6 provides the mixture phase diagram in the pressure-composition plane. The triple point is characterized by % = 1.65(5) MPa , + +G3}z = 8(2) × 10G^ .

713

= 0.750(1) , +G3}z = 0.9789(1) , and

Fig. S2 in Supplementary Information provides the pressure-

composition phase diagrams at = 300, 350, 400, and 450 K. The qualitative features of the

phase diagrams remain the same over the temperature range. The triple point pressure and the concentration of water in the vapor phase are found to increase with increasing temperature.


Page 26 of 61

Page 27 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Pressure-composition diagram for the octane-water mixture at = 450 K. Solid green, blue,

and red curves represent water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid saturation curves, respectively.

Fig. 7 shows the mutual solubilities of the binary system as a function of pressure at a temperature of 450 K. For the aqueous phase, the simulation solubility results are in good agreement with experimental data.45 In contrast, for the octane-rich phase, the models employed here underestimate the water solubility by 60% at the triple point. This finding is consistent with the study of Ballal et al,5 who found that using Lorentz-Berthelot mixing rules for the cross parameters produces a water solubility that is an order of magnitude lower than the experimental value at 298 K and 1 atm. According to Ballal et al, an effective water Lennard-Jones energy is required to match the experimental water solubility in TraPPE alkanes. The reason for this deviation is not fully understood. Apart from this contribution, the assumption of neglecting the long-range LJ potential in our work may also contribute to the inconsistencies we find for the water-in-alkane solubility. Finally, we note that inadequacy of the Lorentz-Berthelot combining



rules in describing the liquid-vapor saturation properties of noble gas mixtures was previously reported by Delhommelle and Millie.46

Figure 7. Mutual solubilities for the octane-water system at 450 K. Black circles represent experiment 45

data

at the triple point. The blue line represents the solubility of water in the octane-rich phase.

Uncertainties are provided at select pressures. The green shaded band represents the solubility of octane in the water-rich phase. The band incorporates the uncertainty in the solubility calculation.


Page 28 of 61

Page 29 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The solubility of octane in water is around 10G^ , which is several orders of magnitude lower than the corresponding solubility of water in octane. Given the very low solubility of octane in water, we consider the possibility of treating the water-rich phase as pure water. This approximation allows for enhanced simulation efficiency, while foregoing knowledge of the octane solubility in the aqueous phase. To test this approach, we construct the water-rich liquidvapor saturation curve at a temperature of 450 K treating the aqueous phase as a single component. Fig. 8 shows the activity difference Δ ln between the approximate and rigorous cases. We find that Δ ln is statistically zero at all . Similar results were obtained for other

saturation properties, such as the pressure and vapor phase composition. The conclusion is also expected to remain valid at lower temperatures for which the solubility of octane in water is even lower. This analysis suggests that saturation properties can be obtained in a rigorous and efficient manner by utilizing the pure water approximation. This simplifying approximation is used in the calculations that follow.



Figure 8. The difference in the saturated activity sum when using the pure water assumption versus the full mixture calculation. The data correspond to the water-rich liquid-vapor saturation line. The band incorporates the uncertainty in the activity calculation.

Temperature Expanded Ensemble. We now turn our attention to the temperature dependence of bulk coexistence properties. When applying the Gibbs phase rule to a binary system, we find that two intensive variables can be independently varied when two phases coexist (e.g., liquid-vapor, liquid-liquid).

Therefore, for a subensemble characterized by a 30


Page 30 of 61

Page 31 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


specific , we vary either the activity fraction , total activity , or pressure % to locate phase coexistence. These cases are examined below. We first constrain the system to a fixed activity fraction . Here, our task is reduced to

tracing the saturation line in (, ) space at constant . We show an example of this approach for calculating the water-rich liquid-vapor saturation curve at = 0.2 .

From previous

calculations, we have obtained the liquid-vapor saturation line (, ln) for pure water as well as

saturated activities of the mixture at a temperature of = 450 K and = 0.2. The initial

(, ln) guess for the mixture is obtained by shifting the (, ln) relationship for pure water by

an amount Δ ln = ln ( = 0.0) − ln ( = 0.2). In other words, we start with the direct

mixture result at = 450 K and = 0.2 and use the variation in temperature from the pure

water (, ln) curve to construct an estimate of the analogous curve for the mixture. Temperature expanded ensemble simulations are applied to compute the relative grand potential

of the two phases over a wide range of temperatures ∈ [268.3 K, 452.7 K] ( ∈ [1.62 ×

10 J G , 2.70 × 10 J G ]). The estimate of the saturated activities at each temperature is This process is repeated until the (, ln)

refined using the approach outlined above.

relationship converges. Results are provided in Fig. 9. Similar strategies are used to deduce the temperature dependence of the octane-rich liquid-vapor and liquid-liquid saturation lines. For the octane-rich liquid-vapor at = 0.937, the initial guess is obtained by shifting the pure octane

saturation curve to pass through the reference saturation point at = 450 K. For the liquidliquid saturation line at = 0.89, we shift both pure octane and water saturation lines to the

reference activity and take the average of the two curves as our initial guess. The iteration process is then repeated to obtain the actual saturation lines.


Data obtained at lower


temperatures, for example = 300 K, 350 K, and 400 K, are used to form initial guesses in

Figure 9. Relationship between activity sum and inverse temperature along the water-rich liquid-vapor coexistence curve at η = 0.2 (green), the octane-rich liquid-vapor coexistence curve at η = 0.937 (blue),

and the liquid-liquid coexistence curve at η = 0.89 (red). Open circles represent data generated from

direct GC or AFEE simulations. Dashed lines represent the initial estimate for the saturation line (these

are difficult to see for the two liquid-vapor saturation curves).

activity fraction expanded ensemble simulations. We now consider an isobaric approach. In practice, phase diagrams are frequently reported at constant pressure, so a simulation-based isobaric approach is of significant interest. In what

follows, we make a connection with experimental data47 collected at % = 1 bar. Within the approach pursued here, we are interested in how and vary with at fixed grand potential (pressure). We aim to deduce the (, , ) relationship along the three saturation curves featured

here at % = 1 bar. We begin by collecting % = 1 bar data from the single-component and mixture AFEE and TEE simulations described above. We then use these data to formulate an approximate (, , ) relationship. Similar to the AFEE case, this approximation provides an 32


Page 32 of 61

Page 33 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


estimate of the triple point (8 = 0.8, ln 8 = −10.8, 8 = 340 ). For the water-rich liquidvapor, octane-rich liquid-vapor, and liquid-liquid saturation lines, we work with ∈ [0, 0.9], ∈ [0.7, 1.0], and ln ∈ [-10.7, -14.0], respectively. Fig. 10 provides the converged saturation curves.

In all cases, data from the simulations described above are consistent with those

generated from the isobaric TEE simulations.

The three saturation lines intersect at 8 =

Figure 10. Relationship between activity sum and activity fraction at % = 1 bar. The green, blue, and red

curves represent water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid saturation curves, respectively. The inset contains data in the vicinity of the triple point.

0.867(2), ln 8 = −10.80(2), and 8 = 344(1) K.

Fig. 11 provides the % = 1 bar phase diagram in temperature-composition space. The

phase diagram is compared with experimental results.47 Simulation results are qualitatively consistent with experiment, but deviate significantly from a quantitative perspective. For the octane-rich region, the underestimate of the saturated temperature is consistent with the 33 ACS Paragon Plus Environment


Page 34 of 61

significant overestimation of the octane vapor pressure discussed above.

As noted earlier,

truncation of the dispersion interaction contributes significantly to this deviation. The water-rich liquid-vapor coexistence region shows the opposite trend.

Here, the SPC/E water model

underestimates the pure-component vapor pressure. Fig. S3 in the Supplementary Information provides temperature-composition phase

diagrams at pressures of % = 0.5, 1.0, and 5.0 bar. For different pressures, the shapes of phase

diagrams are similar. The triple point temperature and vapor-phase water concentration increase with increasing pressure. This trend is consistent with what we found above for pressure-

Figure 11. Temperature-composition phase diagram at % = 1 bar. Solid lines represent simulation data from this work. Circles denote experimental data.

47

Green, blue, and red curves/points represent water-

rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid phase separation, respectively.

composition phase diagrams. Triple Line. The triple curve represents a signature component of mixture phase behavior.

We obtained triple point conditions at = 450, 400, 350, 300 K from the AFEE simulations 34 ACS Paragon Plus Environment

Page 35 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


described above. In this section, we consider how the tripe point evolves with temperature. We

aim to deduce the (, , ) relationship along the triple line. Polynomial fits with the data from AFEE simulations provide an initial guess for the relationships between (, ) and (, ln)

within the temperature range ∈ [268.3 K, 452.7 K] ( ∈ [1.6 × 10 J G , 2.7 × 10 J G ]). TEE simulations are employed to compute the relative grand potential of the three phases of interest. At each temperature, we refine our estimate of the triple point by locating the (, ln) combination that satisfies the criteria for liquid-liquid-vapor phase coexistence. Fig 12 provides the converged triple line that result from the process described above. As is noted, our TEE results agrees well with the data from AFEE simulations. Fig. 13 provides the water solubility in the octane-rich phase and vapor phase. The water concentration in the two phases is found to increase with increasing temperature. The models used here underestimate the water solubility by around 70% for the octane-rich phase and 90% for the vapor phase. Fig. 14 provides the temperature dependence of the vapor pressure for the mixture along the triple line. The triple point pressure is found to increase with increasing temperature. overestimate the vapor pressure by 47% at relatively high temperature.


The models used here


Figure 12. Relationships between activity fraction and inverse temperature (top panel) and activity sum and inverse temperature (bottom panel) along the triple line. Red circles represent data obtained from AFEE simulations. Solid red and dashed green curves represent the final solution and initial estimate for the triple line, respectively.


Page 36 of 61

Page 37 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 13. Water solubility in the octane-rich liquid phase and vapor phase. Open circles and squares represent experimental data. Solid lines denote data from this work. The blue and marron lines represent the solubility of water in the octane-rich liquid phase and the solubility of water in the vapor phase, respectively.



Figure 14. Temperature dependence of the vapor pressure for the octane-water mixture along the triple line. Open black circles represent experimental data for the mixture. The solid purple line represents simulation data for the mixture. Dashed blue and green lines represent simulation data for pure octane and pure water, respectively.

Conclusion We described Monte Carlo simulation methods for computing the saturation properties of a model alkane-water mixture over a broad range of temperatures, pressures, and compositions. We first examined the saturation properties of pure water and octane systems. The coexisting densities and pressures were compared with experimental data and previous simulation results. We observed that the scheme used here to truncate the dispersion interactions had a relatively minor effect on the saturation properties of water and a significant impact on the saturation properties of octane.

A combination of direct calculation and various expanded ensemble

techniques were used to construct phase diagrams for the fluid mixture. Activity fraction expanded ensemble simulations were used to trace the composition dependence of saturation properties along the water-rich liquid-vapor, octane-rich liquid-vapor, and liquid-liquid saturation lines at constant temperatures. The data obtained from these calculations were used to construct pressure-composition diagrams at several temperatures. Various temperature expanded ensemble simulations were applied to determine how saturation properties evolve with temperature under a specified constraint. These constraints included fixed activity fraction, constant pressure, liquid-liquid-vapor coexistence. The data obtained from these calculations were used to construct temperature-composition phase diagrams and to trace the triple line over a wide range of temperatures.


Page 38 of 61

Page 39 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Recent work by Mahynski and coworkers48-49 provides interesting possibilities for improving the efficiency of the approach pursued here.

The expanded ensemble strategy

outlined above begins by developing an estimate of the saturation curve. We typically construct these estimates by completing direct calculations at multiple conditions and subsequently fitting a curve to the data. Mahynski et al. have developed a robust means to generate estimates of saturation curves over a reasonably broad range of conditions from data collected within a single direct simulation. Employing this technique to assemble an initial estimate of the saturation curve would likely prove more efficient than the approach pursued in this paper. Indeed, the expanded ensemble strategy outlined here and the estimation technique developed by Mahynski et al. would provide a nice method for determining saturation properties over a wide range of conditions. We have extended this work to interfacial systems. More specifically, we have used the interface potential approach31, 50-54 to study the wetting properties of the water-octane mixture at a silica surface. A similar overall strategy is pursued in which direct calculations are used to determine interfacial properties at a single state point and expanded ensembles are subsequently employed to obtain properties over a range of state conditions. We will provide the results from this study in a separate report. The general approach pursued here can also be extended to ternary and higher-order mixtures. For example, in a related study we consider how a third component partitions between water-rich and alkane-rich coexisting phases. A variant of the expanded ensemble approach outlined above is used to compute water-rich and alkane-rich solubilities of the third component as a function of the activity fraction of this component at fixed temperature and pressure. The details of this approach will be presented in a separate report. Other types of multicomponent 39 ACS Paragon Plus Environment


phase behavior could also be explored via the methods presented here. The specific approach employed depends upon the nature of the path one is interested in exploring.

Supporting Information Figure S1-S3 discussed in the text.

Acknowledgements We gratefully acknowledge the financial support of the National Science Foundation (Grant No. CBET-1264323) and the American Chemical Society Petroleum Research Fund (Grant No. ND6-51898). Computational resources were provided in part by the University at Buffalo Center for Computational Research.

References 1.

Aparicio-Martínez, S.; Hall, K. R., Phase equilibria in water containing binary systems from molecular based equations of state. Fluid Phase Equilib. 2007, 254, 112-125.

2.

Nasrifar, K.; Moshfeghian, M., Liquid–liquid equilibria of water–hydrocarbon systems from cubic equations of state. Fluid Phase Equilib. 2002, 193, 261-275.

3.

Míguez, J. M.; dos Ramos, M. C.; Piñeiro, M. M.; Blas, F. J., An examination of the ternary methane + carbon dioxide + water phase diagram using the SAFT-VR approach. J. Phys. Chem. B 2011, 115, 9604-9617.

4.

Barrufet, M. A.; Liu, K.; Rahman, S.; Wu, C., Simultaneous vapor−liquid−liquid equilibria and phase molar densities of a quaternary system of propane + pentane + octane + water. J. Chem. Eng. Data 1996, 41, 918-922.


Page 40 of 61

Page 41 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5.

Ballal, D.; Venkataraman, P.; Fouad, W. A.; Cox, K. R.; Chapman, W. G., Isolating the non-polar contributions to the intermolecular potential for water-alkane interactions. J. Chem. Phys. 2014, 141, 064905.

6.

Tsonopoulos, C., Thermodynamic analysis of the mutual solubilities of normal alkanes and water. Fluid Phase Equilib. 1999, 156, 21-33.

7.

Chapoy, A.; Mohammadi, A. H.; Richon, D.; Tohidi, B., Gas solubility measurement and modeling for methane–water and methane–ethane– n-butane–water systems at low temperature conditions. Fluid Phase Equilib. 2004, 220, 111-119.

8.

Kobayashi, R.; Katz, D., Vapor-liquid equilibria for binary hydrocarbon-water systems. Ind. Eng. Chem. Res. 1953, 45, 440-446.

9.

Mohammadi, A. H.; Chapoy, A.; Tohidi, B.; Richon, D., Measurements and thermodynamic modeling of vapor−liquid equilibria in ethane−water systems from 274.26 to 343.08 K. Ind. Eng. Chem. Res. 2004, 43, 5418-5424.

10.

Ma̧czyński, A.; Wiśniewska-Gocłowska, B.; Góral, M., Recommended liquid–liquid equilibrium data. Part 1. Binary alkane–water systems. J. Phys. Chem. Ref. Data 2004, 33, 549-577.

11.

Góral, M.; Ma̧czyński, A.; Wiśniewska-Gocłowska, B., Recommended liquid–liquid equilibrium data. Part 2. Unsaturated hydrocarbon–water systems. J. Phys. Chem. Ref. Data 2004, 33, 579-591.

12.

Góral, M.; Wiśniewska-Gocłowska, B.; Ma̧czyński, A., Recommended liquid–liquid equilibrium data. Part 3. Alkylbenzene–water systems. J. Phys. Chem. Ref. Data 2004, 33, 1159-1188.

13.

Folas, G. K.; Kontogeorgis, G. M.; Michelsen, M. L.; Stenby, E. H., Vapor–liquid, liquid– liquid and vapor–liquid–liquid equilibrium of binary and multicomponent systems with MEG. Fluid Phase Equilib. 2006, 249, 67-74.

14.

Medeiros, M., Liquid–liquid interfacial tensions of binary water-hydrocarbons mixtures via gradient theory and CPA equation of state. Fluid Phase Equilib. 2015, 392, 65-73.

15.

Sun, R.; Lai, S.; Dubessy, J., Calculation of vapor–liquid equilibrium and PVTx properties of geological fluid system with SAFT-LJ EOS including multi-polar contribution. Part III. Extension to water–light hydrocarbons systems. Geochim. Cosmochim. Acta 2014, 125, 504-518. 41 ACS Paragon Plus Environment


16. Papaioannou, V.; Calado, F.; Lafitte, T.; Dufal, S.; Sadeqzadeh, M.; Jackson, G.; Adjiman, C. S.; Galindo, A., Application of the SAFT-γ Mie group contribution equation of state to fluids of relevance to the oil and gas industry. Fluid Phase Equilib. 2016, 416, 104-119. 17.

Lobanova, O.; Mejía, A.; Jackson, G.; Müller, E. A., SAFT-γ force field for the simulation of molecular fluids 6: Binary and ternary mixtures comprising water, carbon dioxide, and n-alkanes. J. Chem. Thermodyn. 2016, 93, 320-336.

18.

Boulougouris, G. C.; Economou, I. G.; Theodorou, D. N., Engineering a molecular model for water phase equilibrium over a wide temperature range. J. Phys. Chem. B 1998, 102, 1029-1035.

19.

Papaioannou, V.; Adjiman, C. S.; Jackson, G.; Galindo, A., Simultaneous prediction of vapour–liquid and liquid–liquid equilibria (VLE and LLE) of aqueous mixtures with the SAFT-γ group contribution approach. Fluid Phase Equilib. 2011, 306, 82-96.

20.

Lymperiadis, A.; Adjiman, C. S.; Galindo, A.; Jackson, G., A group contribution method for associating chain molecules based on the statistical associating fluid theory (SAFT-γ). J. Chem. Phys, 2007, 127, 234903.

21.

Boulougouris, G. C.; Errington, J. R.; Economou, I. G.; Panagiotopoulos, A. Z.; Theodorou, D. N., Molecular simulation of phase equilibria for water−n-butane and water−n-hexane mixtures. J. Phys. Chem. B 2000, 104, 4958-4963.

22.

Ferguson, A. L.; Debenedetti, P. G.; Panagiotopoulos, A. Z., Solubility and molecular conformations of n-alkane chains in water. J. Phys. Chem. B 2009, 113, 6405-6414.

23.

Johansson, E.; Bolton, K.; Theodorou, D. N.; Ahlström, P., Monte Carlo simulations of equilibrium solubilities and structure of water in n-alkanes and polyethylene. J. Chem. Phys, 2007, 126, 224902.

24.

Rane, K. S.; Murali, S.; Errington, J. R., Monte Carlo simulation methods for computing liquid–vapor saturation properties of model systems. J. Chem. Theory Comput. 2013, 9, 2552-2566.

25.

Kumar, V.; Errington, J. R., Monte Carlo simulation strategies to compute interfacial and bulk properties of binary fluid mixtures. J. Chem. Phys, 2013, 138, 174112.

26.

Ferrenberg, A. M.; Swendsen, R. H., New Monte Carlo technique for studying phase transitions. Phys. Rev. Lett. 1988, 61, 2635-2638.


Page 42 of 61

Page 43 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


27. Errington, J. R.; Panagiotopoulos, A. Z., New intermolecular potential models for benzene and cyclohexane. J. Chem. Phys, 1999, 111, 9731-9738. 28. Escobedo, F. A.; Pablo, J. J. d., Expanded grand canonical and Gibbs ensemble Monte Carlo simulation of polymers. J. Chem. Phys, 1996, 105, 4391-4394. 29.

Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D., Hybrid Monte Carlo. Phys. Lett. B 1987, 195, 216-222.

30.

Kumar, V.; Rane, K. S.; Wierzchowski, S.; Shaik, M.; Errington, J. R., Evaluation of the performance of GAFF and CGenFF in the prediction of liquid–vapor saturation properties of naphthalene derivatives. Ind. Eng. Chem. Res. 2014, 53, 16072-16081.

31.

Rane, K. S.; Kumar, V.; Errington, J. R., Monte Carlo simulation methods for computing the wetting and drying properties of model systems. J. Chem. Phys, 2011, 135, 234102.

32.

Kumar, V.; Errington, J. R., Application of the interface potential approach to calculate the wetting properties of a water model system. Mol Simul. 2013, 39, 1143-1152.

33.

Rane, K. S.; Errington, J. R., Using Monte Carlo simulation to compute liquid–vapor saturation properties of ionic liquids. J. Phys. Chem. B 2013, 117, 8018-8030.

34.

Rane, K. S.; Kumar, V.; Wierzchowski, S.; Shaik, M.; Errington, J. R., Liquid–vapor phase behavior of asphaltene-like molecules. Ind. Eng. Chem. Res. 2014, 53, 17833-17842.

35.

Tenney, C. M.; Cygan, R. T., Molecular simulation of carbon dioxide, brine, and clay mineral interactions and determination of contact angles. Environ. Sci. Technol. 2014, 48, 2035-2042.

36.

Jimenez-Angeles, F.; Firoozabadi, A., Contact angle, liquid film, and liquid-liquid and liquid-solid interfaces in model oil-brine-substrate systems. J. Phys. Chem. C 2016, 120, 11910-11917.

37.

Giovambattista, N.; Debenedetti, P. G.; Rossky, P. J., Effect of surface polarity on water contact angle and interfacial hydration structure. J. Phys. Chem. B 2007, 111, 9581-9587.

38.

Martin, M. G.; Siepmann, J. I., Predicting multicomponent phase equilibria and free energies of transfer for alkanes by molecular simulation. J. Am. Chem. Soc. 1997, 119, 8921-8924.

39.

Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P., The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269-6271.



40. Damm, W.; Frontera, A.; Tirado–Rives, J.; Jorgensen, W. L., OPLS all‐atom force field for carbohydrates. J. Comput. Chem. 1997, 18, 1955-1970. 41. Karasawa, N.; Goddard, W. A., Acceleration of convergence for lattice sums. J. Phys. Chem. 1989, 93, 7320-7327. 42.

John, B. S.; Juhlin, C.; Escobedo, F. A., Phase behavior of colloidal hard perfect tetragonal parallelepipeds. J. Chem. Phys. 2008, 128.

43. Economou, I. G., Monte Carlo simulation of phase equilibria of aqueous systems. Fluid Phase Equilib. 2001, 183-184, 259-269. 44. Lagüe, P.; Pastor, R. W.; Brooks, B. R., Pressure-based long-range correction for LennardJones interactions in molecular dynamics simulations: Application to alkanes and interfaces. J. Phys. Chem. B 2004, 108, 363-368. 45. Heidman, J. L.; Tsonopoulos, C.; Brady, C. J.; Wilson, G. M., High‐temperature mutual solubilities of hydrocarbons and water. Part II: Ethylbenzene, ethylcyclohexane, and n‐ octane. AIChE J. 1985, 31, 376-384. 46.

Delhommelle, J.; Millie, P., Inadequacy of the Lorentz-Berthelot combining rules for accurate predictions of equilibrium properties by molecular simulation. Mol. Phys. 2001, 99, 619-625.

47.

Tu, M.; Fei, D.; Liu, Y.; Wang, J., Phase equilibrium for partially miscible system of octane-water. Gaoxiao Huaxue Gongcheng Xuebao 1998, 12, 325–330.

48.

Mahynski, N. A.; Errington, J. R.; Shen, V. K., Temperature extrapolation of multicomponent grand canonical free energy landscapes. J. Chem. Phys. 2017, 147, 054105.

49.

Mahynski, N. A.; Blanco, M. A.; Errington, J. R.; Shen, V. K., Predicting low-temperature free energy landscapes with flat-histogram Monte Carlo methods. J. Chem. Phys, 2017, 146, 074101.

50.

Grzelak, E. M.; Errington, J. R., Computation of interfacial properties via grand canonical transition matrix Monte Carlo simulation. J. Chem. Phys, 2008, 128, 014710.

51.

Grzelak, E. M.; Errington, J. R., Nanoscale limit to the applicability of Wenzel’s equation. Langmuir 2010, 26, 13297-13304.

52.

Kumar, V.; Errington, J. R., Wetting behavior of water near nonpolar surfaces. J. Phys. Chem. C 2013, 117, 23017-23026. 44 ACS Paragon Plus Environment

Page 44 of 61

Page 45 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


53. Kumar, V.; Errington, J. R., Impact of small-scale geometric roughness on wetting behavior. Langmuir 2013, 29, 11815-11820. 54.

Grzelak, E. M.; Errington, J. R., Calculation of interfacial properties via free-energy-based molecular simulation: The influence of system size. J. Chem. Phys, 2010, 132, 224702.



TOC Graphic


Page 46 of 61

Page 47 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


82x81mm (300 x 300 DPI)



82x81mm (300 x 300 DPI)


Page 48 of 61

Page 49 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


82x81mm (300 x 300 DPI)



82x81mm (300 x 300 DPI)


Page 50 of 61

Page 51 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


82x81mm (300 x 300 DPI)



82x81mm (300 x 300 DPI)


Page 52 of 61

Page 53 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


402x310mm (96 x 96 DPI)



402x310mm (96 x 96 DPI)


Page 54 of 61

Page 55 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


82x81mm (300 x 300 DPI)



82x81mm (300 x 300 DPI)


Page 56 of 61

Page 57 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


82x81mm (300 x 300 DPI)



82x81mm (300 x 300 DPI)


Page 58 of 61

Page 59 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


402x310mm (96 x 96 DPI)



402x310mm (96 x 96 DPI)


Page 60 of 61

Page 61 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


338x190mm (300 x 300 DPI)


Calculation of the Saturation Properties of a Model ... - ACS Publications

Recommend Documents