Miscibility at the immiscible liquid/liquid interface: A molecular dynamics study of thermodynamics and mechanism

It is well known that water is sparingly soluble in even the most nonpolar liquids. Water will, in very small amounts, penetrate into the adjacent “immiscible” phase when a nonpolar liquid is in contact with water vapor or liquid water. A small concentration of water in a nonpolar environment can have a significant impact by initiating corrosion, inhibiting particle growth, poisoning catalytic sites, occupying active sites, or dramatically increasing the energetic barrier to chemical reactions.^1–5 For these reasons, the quantification and removal of trace water in organic solvents continues to attract attention in several fields.^1,6,7 The equilibrium concentration of water dissolved in nonpolar solvents has been quantified experimentally and, to a limited extent, studied by computer simulation.^8–11 However the thermodynamics and, in particular, mechanism of water’s dissolution into an adjacent, immiscible liquid phase remain open questions.

In this work, we focus on the transfer of water from a bulk liquid water phase to a bulk liquid organic phase at the immiscible water/liquid interface. Our initial intuitive understanding of this transfer event may draw from studies of ion transfer across that interface.^12–19 In ion transfer, the ionic species moves across the interface and drags with it a portion of its hydration shell. The ion remains tethered to the aqueous phase by an extruded protrusion of hydrogen bonded water molecules that may extend over 1 nm beyond the Gibbs Dividing Surface (GDS) into the organic phase. When the final water-water hydrogen bond breaks, the ion and part of its hydration shell (typically 1-4 water molecules) diffuse into the organic phase. However, water-water hydrogen bonding interactions are typically significantly weaker than ion-water electrostatic interactions, and this mechanism may not accurately explain the dissolution of water in an organic solvent.

A second logical basis of comparison is the evaporation of water. Water evaporation is also the transfer of water from a bulk liquid phase into a medium of much lower relative permittivity (a vacuum) and therefore may share thermodynamic or mechanistic aspects with the transfer of water into a nonpolar liquid phase. A series of recent theoretical studies have revealed interesting details regarding the evaporation of water.^20–23 Bonn and co-workers²² found that the evaporation of water involves two main steps: (1) A final collision, where the transferring water is close to a hydrogen-bonded surface water molecule and (2) the formation of a new hydrogen bond between the non-transferring water and another nearby water just prior to this final collision. They suggest that this mechanism acts as a “Newton’s cradle,” where the formation of the new hydrogen bond imparts sufficient energy to the transferring water molecule to break its last remaining water-water hydrogen bond, ejecting it from the liquid/vapor interface.

Sakaguchi and Morita also investigated the evaporation of water through Langmuir-Blodgett films of n-alcohols, reporting that the alcohols slightly reduce the net free energy required for water evaporation or condensation through the monolayer.²⁴ Presence of a short-chain alcohol monolayer (butanol) has a very small effect on the monotonic Free Energy Profile (FEP) of water evaporation but long-chain n-alcohols (decanol and hexadecanol) introduce a free energy barrier of a few kcal/mol due to a significant decrease in the mass accommodation coefficient. The FEPs of evaporation through long-chain films resemble the transfer of an ion across the immiscible liquid/liquid interface when the interface is constrained to a flat plane by an externally applied potential.¹⁵ This relationship is interesting since both the interface constraint and the intra-monolayer Van der Waals forces effectively flatten the interface, suppressing fluctuations. Sakaguchi and Morita report that the transport of water through the monolayer is facilitated by a “water finger” protrusion, further strengthening the link to liquid/liquid ion transfer.

To investigate the water dissolution process, we first use equilibrium molecular dynamics simulations to derive the FEP, also known as the Potential of Mean Force (PMF), for the transfer of water molecule(s) across the interface. These calculations provide a mean field description and thus largely ignore the importance of dynamic fluctuations at the liquid/liquid interface. However, as in the case of ion transfer, the FEP could provide evidence that these fluctuations are important to transport across these adjacent phases. We recently studied the free energy of transfer of a chloride ion across the water/nitrobenzene interface, finding that when interfacial fluctuations are suppressed (a flat interface constraint is implemented but co-transfer of the water’s solvation shell is permitted), the net free energy of transfer increases by approximately 60%.¹⁵ In the case of ion transfer, these and other computer simulations^{13,14,16–18} strongly support and quantify the “water finger” or protrusion mechanism. Also in the present work, these equilibrium FEP calculations are accompanied by a few suggested “reaction coordinates”: collective variables that attempt to describe the transfer event using a few degrees of freedom.

We then complement these equilibrium studies with a set of 487 molecular dynamics trajectories, each having captured an independent dissolution event of water into nitrobenzene. Our analysis of these dissolution events is similar to the water evaporation studies recently performed by Bonn and co-workers,²² which allows for easy mechanistic comparison between water/vapor (evaporation) and liquid/liquid (dissolution) water transfer.

The rest of this paper is organized as follows: In Sec. II, we describe and analyze the thermodynamics calculations performed on the interfacial systems. Section III analyzes the collected dissolution events to investigate the mechanism and dynamics of water transferring from bulk water to nitrobenzene. Section IV summarizes and concludes.

Molecular dynamics simulations of water/nitrobenzene, water/hexane, and water liquid/vapor systems were performed using our in-house code. The liquid/liquid systems consist of two adjacent slabs of immiscible liquids in a 31.28 × 31.28 × 300.0 Å rectangular box with the liquid/liquid interface located in the x-y plane at z ≈ 0 Å. The aqueous phases, consisting of 986 water molecules, are located at z < 0 and the organic phases, either 252 nitrobenzene molecules or 250 hexane molecules, are in the z > 0 region. Each liquid phase is in equilibrium with its respective vapor phase with soft reflecting potential walls located 5 Å from the simulation box boundaries in the z direction to prevent vapor phase molecules from crossing into the adjacent vapor phase. The water liquid/vapor system consists of a slab of 986 water molecules in equilibrium with its own vapor contained within a rectangular box with the same dimensions and soft reflecting walls near the z boundaries as the liquid/liquid systems.

All intermolecular interaction potentials are the pairwise sum of Lennard-Jones and Coulomb terms

where r is the distance between atom centers i and j. Lennard-Jones parameters for mixed interactions are generated using Lorentz-Berthelot combining rules, $εij=εiεj1/2$ and $σij=(σi+σj)/2$⁠. The nitrobenzene model has been described previously²⁵ and reproduces the experimental dipole moment and enthalpy of vaporization. For hexane, we use the Optimized Potentials for Liquid Simulations united atom (OPLS-ua) model.²⁶ Our water model is a version of the flexible simple point charge (SPC) model²⁷ with intermolecular potentials as described by Kuchitsu and Morino.²⁸ These particular combinations of water and organic models produce liquid/liquid and liquid vapor interfacial properties (surface tension, density profiles, and orientational distributions) that agree well with experiment.^13,29–31

Starting configurations were first constructed manually and equilibrated by performing a brief MD simulation with short time steps (dt < 0.1 fs) to relax any high energy, nonphysical configurations. This was followed by at least 1 ns of constant temperature simulation with a time step of 0.5 fs to ensure equilibration (e.g., convergence of average total energy, interaction energy between the two liquids, density profiles, etc.) and verify that all applied constraints conserve energy.

We use umbrella sampling^32–34 with an applied biasing potential to calculate the free energy profile of a water molecule transferring from bulk water into the adjacent phase. Our implementation of this technique has recently been described in detail¹⁵ and here we briefly summarize it.

To calculate the free energy profile of this water transfer process, we first select a single water molecule of interest: the “tagged” water. The interval of interest along the z-axis spans the region from bulk water to the bulk organic liquid, −10 Å > z_w > 17 Å, where z_w is the position of the transferring water molecule’s center of mass along the z-axis, and the liquid/liquid interface is located at approximately z = 0 Å. We divide this interval into a series of N overlapping lamellae that are parallel to the x-y plane (parallel to the liquid/liquid interface). The local free energy in each lamella is

where P_n(z_w) is the probability distribution of z_w within lamella n and C_n is a constant. The A_n segments are combined by selecting C_n to minimize the difference between A_n and A_n+1 in the overlapping regions^35,36 to create a complete free energy profile which spans the entire interval of interest.

Poor statistical sampling of a given lamella may occur if the value of A(z_w) varies by a large amount. We apply a biasing potential, U_bias(z_w), to the tagged water to accelerate the exploration of phase space and improve sampling statistics. This modifies the interaction between the tagged water molecule and the surrounding solvent molecules U_wx as follows:

where U_bias(z_w) is a function of the position of the tagged water molecule’s center of mass along the z-axis and is of the form

where the constants α, γ, and ζ are selected so that the biasing potential approximates A(z_w). The free energy profile within a lamella is determined from the biased probability distribution $Pnb(zw)$ by

The simulations in this work use lamellae that are 3 Å wide and overlap by 1 Å. The tagged water is confined within the specified lamella by a window potential that is zero within the boundaries of the lamella along the z-axis and rises rapidly if the tagged water attempts to exit the lamella. All equilibrium data presented in this work represent 10 ns of sampling within each of the overlapping lamella, and all simulations are performed at 298 K.

Figure 1 shows the density profiles of the two liquid/liquid systems. Obviously water is immiscible with each of these adjacent organic phases. In each system, the Gibbs Dividing Surface (GDS) is located at z_GDS = 0 Å. The exact definition of the Gibbs dividing surface is the plane parallel to the interface where the excess density of one liquid on one side exactly matches the deficit on the other side relative to a step function profile.³⁷ However, at the liquid/liquid interface, one may define the GDS with respect to each liquid and get a slightly different value, which to a good approximation is the location where the density of water is 50% of the bulk value.

The interface region of the water/nitrobenzene system [Fig. 1(a)] is noticeably wider than the water/hexane interface [Fig. 1(b)]. The distance over which the water density varies from 90% to 10% of the bulk value is approximately 5 Å in the water/nitrobenzene system and 3 Å in the water/hexane system. This difference in interface width can be directly related to the surface tension of two interfaces: σ = 51 dyn/cm for water/hexane³⁸ and σ = 25 dyn/cm for water/nitrobenzene.³⁹ The interface width is a result of thermal averaging over capillary fluctuations of a molecularly sharp interface.^37,40 All phases achieve bulk density when only a few molecular diameters away from the liquid/liquid interface.

The miscibility of water in the two organic liquids is determined by energetic and entropic contributions but is roughly related to the relative permittivity of the two liquids: ε_r ≈ 6 for hexane and ε_r = 34.8 for nitrobenzene at 298 K. The exact value can be obtained from the net difference in free energy as determined by the FEP calculations. Figure 2(a) shows the FEP of a water molecule transferring from the bulk aqueous phase to the bulk region of the adjacent phase in each of the three interfacial systems, calculated using the umbrella sampling method described in Sec. II A 2. As the tagged water molecule moves from bulk water toward the organic phase (in the +z direction), each free energy curve increases monotonically and plateaus in the organic region. The water/nitrobenzene FEP (blue curve) begins to rise at z ≈ −3 Å (about one water layer before the GDS), exhibits a slight and broad maximum at z ≈ 7 Å and reaches a plateau at z ≈ 10 Å. The net free energy of transfer is 3.8 kcal/mol. The experimentally determined solubility of nitrobenzene in water, 0.178 M,⁹ corresponds to a net free energy of transfer of 2.36 kcal/mol. This difference in net free energy is not large in magnitude but is a significant percentage of the actual value. The water/hexane FEP (red curve) plateaus at z ≈ 7 Å and had the largest net free energy of transfer, 7.6 kcal/mol. The experimental solubility data, 0.0362 M water in hexane at 298 K,⁸ correspond to a net free energy of transfer of 4.54 kcal/mol, significantly lower than the value determined by our simulations. Finally, the green curve in Fig. 2(a) represents evaporation of water, the transfer of a tagged water molecule from bulk water into an adjacent water vapor phase, also plateau at z ≈ 7 Å. The net free energy of transfer for evaporation is 6.5 kcal/mol, which compares well to the experimentally determined value of 6.32 kcal/mol.¹¹ 

The difference between the water/hexane liquid/liquid net free energy in Fig. 2(a) and the corresponding value determined from solubility data reflects a recent and ongoing discussion regarding the accuracy of the force fields used in water/alkane simulations. Chapman and co-workers¹⁰ identified the poor performance of fixed-charge force fields and Lorentz-Berthelot combining rules when studying the solubility of water in alkanes, noting that the water-alkane interactions are dramatically underestimated. This work was followed by commentary by McDaniel and Yethiraj⁴¹ and a similar report (with similar findings) by Panagiotopoulos and co-workers⁴² on water/CO₂ and water/n-alkane mixtures. Recently, subsequent work by Chapman and co-workers¹¹ confirmed that the electrostatic and induction effects between water and alkanes are critical components to consider when considering this extreme case of water’s solubility in alkanes. Our water/hexane free energy profile in Fig. 2(a) agrees with the underestimation of water’s solubility in an alkane. The solubility of water in nitrobenzene is also underestimated but to a much lesser extent. Since nitrobenzene (unlike hexane) is significantly polar, the electrostatic interaction energy between nitrobenzene and water at the liquid/liquid interface is much larger than the interaction energy resulting from the induced dipoles when water and nitrobenzene are in close contact. As a result, it is possible to account for the total energy using fixed empirical parameters.

The calculated net free energy of transfer of our non-polarizable water across the liquid/vapor interface compares very well with experimental values. This thermodynamic result agrees with the mechanistic simulations performed by Bonn and co-workers, who directly compared the performance of a nonpolarizable SPC water model with an ab initio based polarizable water force field.²² They reported that their conclusions hold equally well regardless of force field.

Figure 2(b) shows the average number of water molecules in the solvation shell of the transferring water. We define surrounding water molecules to be a member of the tagged water molecule’s solvation shell if the oxygen-oxygen distance is less than 3.4 Å, which is the first minimum of g_OO(r) in bulk water. In bulk water, the tagged water has about 5 water molecules in its solvation shell. This population decreases quite rapidly as the tagged water moves to the interface and eventually reaches zero when the final water-water hydrogen bond breaks and the transferring water is “dissolved” in the adjacent phase. The blue curve, representing the water/nitrobenzene phase, maintains a nonzero population deeper into the adjacent phase than in the hexane and water vapor systems, which parallels the behavior of the FEP curves shown in the top panel. The more polar nitrobenzene phase allows the transferring water to drag neighboring waters further into the organic phase, suggesting that the transfer of water into an adjacent, immiscible, weakly polar phase is accompanied by a larger local protrusion of the aqueous phase into the organic phase. This protrusion appears to be larger than in the case of water/nitrobenzene than in water/hexane or water/water vapor but much smaller than the protrusions or “water fingers” that accompany the transfer of a small ion across the immiscible water/organic interface.^13,15,16 The rapid decay to zero suggests that in all systems, the transfer of a single water molecule is energetically favored over the simultaneous transfer of several waters. This is the same behavior as in the case of evaporation^20,22 but is in marked contrast with the transfer of negative or positive ions,^15,18 which keep part of their hydration shell due to strong water-ion electrostatic forces. Nonzero solvation shell values in the region of z ≈ 15 Å exist in the water/nitrobenzene system but are not due to co-transferred water. In these cases, other water molecules have dissolved into the nitrobenzene phase and “found” the transferred water molecule.

To further support the notion of one-at-a-time transfer of water in the water/nitrobenzene system, we computed the FEP for the transfer of a water “dimer” ((H₂O)₂) and a water “trimer” ((H₂O)₃). In each of these co-transfer simulations, we tethered one or two of the surrounding water molecules to the tagged water by implementing a spherical, harmonic restraining potential around the tagged water’s center of mass that kept the center of masses of the tethered and tagged waters within 3.4 Å of each other. (We note that the restraining potential is applied to the centers of mass and not, e.g., the oxygen atom centers so that the orientation of the molecules is not affected by the introduction of a tethering force.)

Figure 3(a) shows the free energy of transfer along the z-axis for the unconstrained system [solid curve, the same as the blue curve in Fig. 2(a)], (H₂O)₂ “dimer” (dashed curve), and (H₂O)₃ “trimer” (dotted curve). The results confirm that the transfer of a single water molecule is preferred since the addition of co-transferred waters increases the net free energy of transfer. As mentioned earlier, the net free energy of transfer for water across the water/nitrobenzene interface is 3.8 kcal/mol. This value increases to 5.9 with the co-transfer of 1 solvation shell water and further increases to 7.4 with the co-transfer of 2 solvation shell water molecules. Figure 3(b) shows the corresponding average number of water molecules in the transferring water’s hydration shell. In bulk water, the restraining potential results in a noticeably larger solvation shell population in the (H₂O)₂ and (H₂O)₃ systems. Aside from this artifact, the solvation shell populations of the (H₂O)₂ and (H₂O)₃ systems gradually reduce toward 1 and 2, respectively, indicating that even the imposition of an artificial solvation shell does not result in the co-transfer of additional water molecules.

At the simplest level, the observation that water dissolves one molecule at a time and does not co-transfer additional water molecules like ions do can be attributed to the much weaker water-water vs. water-ion electrostatic interactions. However, the increase in free energy of transfer from 3.8 kcal/mol to 5.9 kcal/mol represents a more complicated balance of energy and entropic contributions, which a detailed molecular theory could provide.

A major component of the free energy of transfer is the interaction energy between the single water molecule and the rest of the system. We denote by U_W(z) and U_O(z) the average interaction energy between the tagged water molecule at the location z and the rest of the water (W) and organic solvent (O), respectively. The bulk values are denoted by $UWbulk=UW(−∞),UObulk=UO(∞),$ all considered at the same temperature. These are determined by molecular dynamics simulations of a single water molecule in the uniform bulk liquid using the same force field parameters described in Sec. II A 1. We find $UWbulk$ = −24.5 kcal/mol, $UHexbulk$ = −2.0 kcal/mol, and $UNBbulk$ = −9.8 kcal/mol. These values suggest that for transfer to nitrobenzene $ΔUW→NB$ = −9.8 − (−24.5) = 14.7 kcal/mol, $ΔSW→NB$ = (ΔU_t – ΔA_t)/T = 37 cal/mol ⋅ K. For transfer to hexane: $ΔUW→Hex$ = −2.0 − (−24.5) = 22.5 kcal/mol, $ΔSW→Hex$ = (ΔU_t – ΔA_t)/T = 50 cal/mol ⋅ K. The entropy increase corresponds to the much greater translational freedom of the water molecules in the organic solvents than in bulk water. The smaller value in nitrobenzene compared with hexane is due to the ordering imposed by the water-nitrobenzene hydrogen bonding.

For the interaction of (H₂O)₂ with the two bulk phases, we find $UD-Wbulk$ = −42.3 kcal/mol and $UD-NBbulk$ = −17.18 kcal/mol, values that are slightly less negative than twice the single water molecule values. For the trimer, we get $UT-Wbulk$ = −60.0 kcal/mol and $UT-NBbulk$ = −20.4 kcal/mol. These values suggest that for transfer of (H₂O)₂ to nitrobenzene, $ΔUD-W→D-NB$ = −17.18 − (−42.3) = 25.1 kcal/mol, $ΔSD-W→D-NB$ = (ΔU_t – ΔA_t)/T = 64 cal/mol ⋅ K. For the transfer of (H₂O)₃ to nitrobenzene, $ΔUT-W→T-NB$ = −20.4 − (−60.0) = 39.6 kcal/mol, $ΔST-W→T-NB$ = (ΔU_t – ΔA_t)/T = 108 cal/mol ⋅ K. As additional water molecules are co-transferred, the increase in the unfavorable energy change outweighs the favorable increase in entropy leading to the increase in the free energy of transfer.

As the tagged water molecule is transferred across the interface, U_W and U_O vary. While the number of water molecules in the first hydration shell of the tagged water molecule provides useful information about the change in the local environment experienced by this water, a more accurate indicator utilizes the total interaction energy defined above. To this end, we define a solvation coordinate s for the transferred water molecule in each solvent,^14,15

where O is “NB” for nitrobenzene and “Hex” for hexane.

Figure 4(a) shows the value of s_W for all simulated systems where $...$ indicates the ensemble average value collected as a function of the tagged water’s position along the simulation z-axis. Since the first solvation shell constitutes most (about 70%) of the solvent-solute interactions, the curves in Fig. 4(a) closely track the corresponding solvation shell population curves shown in Figs. 2(b) and 3(b). For all systems, the value of $sW$ begins at 1, while the tagged water molecule is in the bulk aqueous region and decreases as it moves into the adjacent phase. For the systems without imposed co-transfer of additional water molecules (solid curves), this decrease is quite dramatic once the tagged molecule crosses the interface since the tagged molecule is effectively stripped of its solvation shell waters. In the cases of (H₂O)₂ and (H₂O)₃, the decrease is more gradual and the value of $sW$ plateaus at a value significantly higher than zero, indicating the large contribution to interaction energy made by solvation shell waters. For example, in the case of (H₂O)₃, two solvation shell waters are present at approximately z = 12 Å and $sW$ ≈ 0.5. This indicates that the strength of the interaction between these two solvation shell waters and the tagged water is about half of the total tagged water-solvent interaction energy in a bulk water environment.

Figure 4(b) shows $sO$ along the simulation z-axis for each of the liquid/liquid systems. For all systems, $sO$ begins at zero in the bulk aqueous phase and increases monotonically, approximately linearly, until the tagged water has transferred into the organic phase. In the water/nitrobenzene systems (blue curves), the early part of this increase (in the region of −4 Å < z < 3 Å) is very similar regardless of the co-transferring water constraint. In the water/hexane system (red curve), the increase of $sO$ begins earlier along the z-axis because the hydration shell is stripped from the transferring water earlier (along the z-axis) than in the water/nitrobenzene system.

To gain geometric insight into the perturbation of the water structure during the transfer, we also describe the transfer process using the “water finger” coordinate w, introduced recently by Morita and co-workers,¹⁶ in the study of ion transfer across the immiscible water/organic liquid/liquid interface. Here we slightly modify the definition of the water finger coordinate to adapt it for the water transfer event. The coordinate w is an instantaneous measure of the transferring species’ connectivity to the aqueous phase. Its value is the largest distance that must be traversed if one were to start at the transferring species and move to the bulk aqueous phase, where the only allowable moves are “hops” to other water molecule oxygen atoms. To calculate this coordinate, an undirected mathematical graph is defined in which each water oxygen atom is a vertex and the distances between the oxygen atom centers are the edges. The coordinate w is the minimum threshold distance required to give a connected path between the transferring water and the aqueous phase. Figures 5(a) and 5(b) show the distance described by w (yellow) in representative system snapshots. When the transferring water is connected to the aqueous phase by a network of water-water hydrogen bonds [Fig. 5(a)], the value of w is approximately the O–O distance of a water-water hydrogen bond. After the last water-water hydrogen bond breaks [Fig. 5(b)], the value of w increases rapidly. We refer the reader to Refs. 15–17 and the supplementary material of this work for more detailed, technical discussions of w.

Figure 5(c) shows the ensemble-averaged values of w as a function of the tagged water’s z position for the interfacial systems. Since the transfer of water occurs one-at-a-time, the value of w in the region of z > 5 Å (tagged water moving into the organic phase) is equal to the length of the final water-water hydrogen bond (when intact) and, after this last bond breaks, w increases rapidly, approximately linearly with the tagged water molecule’s z position. The most obvious feature distinguishing the systems in Fig. 5(c) is the spatially delayed increase of $w$ in the water/nitrobenzene system (blue curve). This corresponds to the tagged water molecule separating from the aqueous phase further into the adjacent phase (recalling that the GDS is located at z ≈ 0 Å). Larger interfacial protrusions (“water fingers”) accompany water transfer in the water/nitrobenzene system due to the lower interfacial surface tension and the final water-water hydrogen bond breaks when the tagged water is at approximately z = 6.5 Å, as opposed to z = 5.9 Å in the case of the water/hexane and water/vapor interfaces. This suggests that water transferring into nitrobenzene does so from a larger protrusion than in the other systems and this will be considered further in Sec. III.

To study the dynamics and mechanism of dissolution of water into nitrobenzene, we use an approach similar to the one used by Mason²⁰ and Bonn and co-workers²² to study water evaporation. We selected 43 independent, equilibrated water/nitrobenzene configurations, ensuring that each starting configuration contains no water molecules dissolved in the nitrobenzene phase. Trajectories are run at constant total energy for 100 ps from each starting configuration by selecting random velocities from a Boltzmann distribution and storing the positions and velocities of the 200 water molecules with the highest z position and 100 nitrobenzene molecules with the lowest z position every 4.0 fs. This selection captures the interfacial region of interest, approximately −6 Å < z < 17 Å in the simulation box (z_GDS ≈ 0). A dissolution event is defined whenever any water molecule reaches z > 10 Å, and the nearest water molecule is over 6 Å away, at which point the trajectory is terminated. These cutoffs ensure that the transferring water has diffused into the organic phase. In the case of no detected dissolution event within the 100 ps interval, new random velocities were selected and the process repeats. In total, 487 dissolution events were collected from 5302 trajectories, which correspond to a total of 507 ns of simulation time.

Our analysis of the collected dissolution events employs a methodology similar to that used by Bonn and co-workers.²² In the final configuration of the recorded trajectories, the “dissolved” water molecule that has transferred into the nitrobenzene phase is labeled “A.” We then travel backwards in time in the trajectory to locate the last water-water interaction prior to the dissolution into nitrobenzene. When the distance between the oxygen atoms of water A and another water is less than 3.3 Å and this water-water distance is the minimum for at least 50 fs while moving backwards through the acquired trajectory, we define this configuration to be t = 0, the moment of the last water-water interaction, and define water B: the last water molecule to “interact” with water A prior to dissolution. This procedure ensures that we do not misinterpret O–O bond “vibrations” as the actual final breakup. The vector pointing from O_B to O_A at time zero, r₀, and the midpoint of the O_B–O_A line segment at time zero, M, are stored and are used as reference points for the dynamics of the dissolution events. Figure 6 is a representative snapshot of a dissolution event, indicating the positions of water molecules A and B. Figure 6(a) was acquired just prior to the breakup of the A–B hydrogen bond (green dotted line). In Fig. 6(b), 40 fs later, the molecules have moved and the hydrogen bond is no longer intact.

Before discussing the dissolution mechanism in detail, a brief comment concerning the dissolution rate is warranted. Since 487 dissolution events were collected in a total of 507 ns trajectories, the rate of dissolution we observed is 0.961 water molecules/ns, which, given the simulation box cross section of (3.128 nm),² corresponds to a dissolution flux of 163 mol∙m⁻²∙s⁻¹. This may be compared with the evaporation flux by using the fact that at equilibrium, the number of evaporating water molecules is equal to the number of collisions of vapor phase molecules with the surface. The later is given by $ρv/4$⁠,⁴³ where $ρ$ is the water vapor density (approximately equal to p_vap/RT with p_vap being the vapor pressure at temperature T) and $v$ = (8RT/πM)^1/2 is the average molecular velocity.⁴³ Using the experimental vapor pressure of water at 298 K, p_vap = 23.8 Torr, we find $ρ$ = 1.28 mol∙m⁻³ and the average velocity $v$ = 592 m/s and thus an evaporation flux of 189 mol∙m⁻²∙s⁻¹ at 298 K. This is surprisingly similar to the dissolution flux given the difference between the free energy barrier to evaporation (6.5 kcal/mol) and dissolution (3.8 kcal/mol). It is possible that our dissolution flux is underestimated because of our conservative definition of the dissolution event, which requires that the transferred water diffuses a few Å into the organic phase. The low dissolution flux can also be due to recrossings of the transition state barrier caused by collisions with nearby solvent molecules, a much less likely event for evaporation.

To examine the dissolution mechanism, we begin by considering several kinematic parameters that allow a direct comparison with the mechanism of water evaporation described in Ref. 22. Figure 7(a) is a cartoon schematic showing some of the relevant geometrical coordinates, and panels b-d show the time dependence of several parameters characterizing the kinematics of the “last collision” before dissolution during the time interval t = 0 ± 100 fs. In all panels, the red curves represent water A (the dissolving water) and the blue curves represent water B (the last water molecule to interact with A). The error bars represent one standard deviation of 487 dissolution events. Figure 7(b) shows the angle θ_A(t) between the instantaneous center of mass velocity vector of A, denoted by v_A(t), and r₀, defined above (red curve) and the same for molecule B (blue curve). The curves in 7b confirm the collision of waters A and B at time zero and show behavior quite similar to the final water-water collision in the evaporation of water.²² About 20 fs prior to the collision, each molecule approaches the r₀ vector about 30° degrees from normal, they collide, moving past the θ = 90° point and exit the collision at nearly the same angle at which the other water entered, resulting in a nearly symmetrical region at t = 0 ± 20 fs.

Figure 7(c) shows the total kinetic energy (⁠$Ek=∑imivi2/2$⁠, where the sum is over the three water atoms) of each water molecule during the final collision. On average water B enters with more kinetic energy and appears to transfer part of this energy to water A during the collision. In this regard the E_k curves are similar in shape to related data obtained for the evaporation of water in Ref. 22. The main difference between the evaporation and dissolution E_k curves is that on average, water B’s kinetic energy remains constant (near the minimum value) after the evaporation collision but water B’s kinetic energy is partially restored after the dissolution collision. This “recovery” of B’s average E_k partially reflects a major mechanistic difference between evaporation and dissolution. In evaporation, water A is “launched” into the vapor phase from the surface by a final collision event occurring near the GDS and water B remains near the GDS. In contrast, as our equilibrium calculations suggested and will further shown below, dissolution into nitrobenzene is accompanied by larger interfacial protrusion of water into the nitrobenzene phase. The pendant water molecule at the “tip” of this protrusion is essentially inserted or deposited into the organic phase and the final water-water interaction exists somewhere between the two possible extremes: (1) Water A is ejected from the tip of the protrusion by the final collision. (2) The protrusion retracts, leaving behind water A. It is between these two extremes that we see the average kinetic energy behavior as shown in Fig. 7(c). In (1), E_k of B would remain constant after the collision while A would gain significantly, and in (2), E_k of A would remain constant after this final collision while E_k of B somewhat increases.

Figure 7(d) shows the average distance of waters A and B from the collision midpoint M during the final collision. Both water A and water B are further from the collision midpoint at +100 fs (after the collision) than they were at −100 fs (water A much more noticeably). We again compare these results with the evaporation data in Ref. 22 where water A, the transferring water, is much further from M at +100 fs, but water B is considerably closer to M at +100 fs than at −100 fs. This reveals further mechanistic insight: In the evaporation collision, water B stays near the collision site, having “launched” A into the vapor phase. In the dissolution event, both waters B and A recoil from the collision at speeds greater than their approach, A having gained momentum from B and B being pulled back into the aqueous phase by its neighbors in the retracting water protrusion (to which B remains hydrogen bonded).

We next examine in detail the hydrogen bonding involved in the dissolution process. As outlined in Sec. I, Bonn and co-workers recently reported on the importance of water-water hydrogen bonding dynamics in water evaporation.²² In their simulations, water B forms a new hydrogen bond with a nearby water molecule at ∼44 fs before the final collision (t = 0). The energy released by the formation of this new hydrogen bond is transferred to the A–B hydrogen bond, the breakup of which results in the evaporation of water A. We use a similar approach in our analysis and also include for the possibility of weak water-nitrobenzene hydrogen bonding.

For water-water hydrogen bonds, we use a common geometrical definition where a hydrogen bond is present if the oxygen-oxygen distance is less than 3.4 Å and the H–O⋯O angle is less than 30°.¹⁵ Note that this definition is more restrictive than the condition for dissolution, which only considers the oxygen-oxygen distance. Figure 8 shows the time variation of average hydrogen bond populations around the time of the final water-water collision. The green curve in Fig. 8 represents the hydrogen bond between waters A and B. Since at time t the A–B hydrogen bond either exists (h(t) = 1) or does not (h(t) = 0), this curve represents the ensemble average $ht$ during the 487 dissolution trajectories. The curve is relatively flat for several hundred fs prior to the detected collision and begins to rapidly decay to zero after about 40 fs. The red and blue curves are the average number of hydrogen bonds between waters A and B and all other water molecules, respectively (not counting A–B hydrogen bonds). The average number of water hydrogen bonds with water A is near zero, indicating that the transferring water only has a strong interaction with water B and suggesting that the predominant final configuration has water A tethered to the aqueous phase by 1 hydrogen bond. Water B participates in an average of 1.8 hydrogen bonds with other water molecules, and this number slightly increases after the dissolution of the A–B hydrogen bond. Figure 8 also considers hydrogen bonding between waters A and B and the surrounding nitrobenzene molecules. For water-nitrobenzene hydrogen bonds, we use the same definition as water-water hydrogen bonds above, where nitrobenzene may participate as a hydrogen bond acceptor. The hydrogen bond population of nitrobenzene-water A (black curve) increases slightly as water A is transferred into the organic phase and able to form the bond more freely. At the same time, the water B-nitrobenzene hydrogen bond population (orange curve) also slightly increases due to the increased exposure of water B to the nitrobenzene phase.

The hydrogen bonding populations in Fig. 8 suggest a rather simple mechanism. Water A is tethered to the aqueous phase by a single hydrogen bond with water B, this hydrogen bond dissociates, and A diffuses into the nitrobenzene phase.

As the final A–B hydrogen bond elongates and breaks, it may adopt two configurations: water A may act as either a hydrogen bond donor or acceptor. The top two panels in Fig. 9 show representative snapshots at t = 0, where water A is acting as a hydrogen bond acceptor [Fig. 9(a)] or donor [Fig. 9(b)]. Hydrogen bonds are indicated by the green dotted lines and nitrobenzene molecules are made invisible so that the difference in the two configurations is clear.

Figure 9(c) shows the average number of A–B hydrogen bonds during the final collision and dissolution event (green curve, the same as the green curve in Fig. 8), and the role of water A in the A–B hydrogen bond is specified by red (acceptor) and blue (donor) curves. We find that the transferring water is the acceptor in about 63% of the final A–B hydrogen bonds. The fact that twice as many dissolving water molecules break an acceptor bond than a donor bond may be attributed to the interaction of water A with nitrobenzene. When A acts as an acceptor, its two positively charged hydrogen atoms face the nitrobenzene phase where they may interact with the negatively charged nitro groups, whereas when A acts as a donor, only one hydrogen atom is interacting with the negative end of the nitro group.

Figure 10 provides further support for the mechanism described above by examining some of the relevant distances between the dissolved water and other nearby molecules. Figure 10(a) shows the shortest distance between water A’s oxygen and the nearest water oxygen (blue curve) and the shortest distance between a water A hydrogen and the nearest water hydrogen atom (red curve). Both the O–O and H–H distances remain rather constant as the collision event is approached, forming a minimum at t = 0 (an artifact of the way the collision is defined, as explained above), and then monotonically increasing, signaling the breakup of the hydrogen bond and subsequent diffusion of the water molecule. The other curves in Fig. 10(a) focus on the interactions with the organic liquid. The green and black curves depict the distance between the transferring water molecule’s oxygen atom and the nearest nitrobenzene oxygen and any nitrobenzene atom, respectively. These curves are flat, indicating that the surrounding nitrobenzene molecules do not appear to engage water A during its final collision and dissociation. The water-nitrobenzene curves also indicate that the negatively charged nitro oxygen atoms are not the nitrobenzene moiety nearest to the transferring water as may be expected when considering the partial charges used in our model.

Finally, in Fig. 10(b), we provide the justification for the assertion that the transferring water emerges from a tip of a water finger. This figure shows the average position of the water A oxygen on the z-axis (red curve) and the distance between waters A and B (blue curve) around the time of the final collision. Again using 3.4 Å as the oxygen-oxygen cutoff distance for both hydration shell membership and maximum water-water hydrogen bond length, we find (average of the 487 transfer events) that the transferring water separates from the bulk at approximately 95 fs after the final collision and at an average z position of 6.0 ± 0.3 Å. This location is only slightly closer to the GDS than the location where the average finger coordinate reaches the cutoff of 3.4 Å: $zw=3.4 Å$ = 6.5 Å in our equilibrium calculations (see Fig. 5). The fact that these two locations are comparable in value is not surprising since in a simple, single-water transfer event w ≈ r_OO = 3.4 Å when the transfer is nearly complete [blue curve, Fig. 8(b)]. Symbolically, $zw=3.4 Å=6.0±0.3 Å≈zw=3.4 Å=6.5 Å$⁠. Thus, the breakup of the water molecule from bulk water during the non-equilibrium dissolution trajectories is taking place from a protrusion’s tip whose location, on average, is very similar to the breakup position as determined by the equilibrium free energy transfer calculations.

The transfer of water from the aqueous phase to an adjacent, immiscible liquid phase has not received much attention, but the thermodynamics and mechanism of this transfer event share features with two more well-studied systems: water evaporation and ion transfer across the liquid/liquid interface. Our molecular dynamics simulations show that, like evaporation, water prefers to transfer to the adjacent phase one-at-a-time: The net free energy of transfer for co-transfer with a single member of its hydration shell [an (H₂O)₂ dimer] is larger than single water transfer by about 3.5 k_BT. In addition to the transferring water’s position relative to the GDS, the normalized solvation energies may be used as reasonable reaction coordinates that correlate strongly with position-based free energy profiles.

Dynamic studies that focused on the interaction between the dissolving water molecules and other interfacial water molecules prior to the transfer event show that the local structure surrounding the transferring water is more similar to a transferring ion than evaporating water. The transfer from water to nitrobenzene occurs via a protrusion into the organic phase, with the last water-water hydrogen bond breaking at about 6 Å beyond the GDS. This breakup occurs approximately 95 fs after the last strong water-water interaction and (unlike in the evaporation of water) involves the dissolution of a pendant water molecule from a water protrusion into the organic phase and retraction of the protrusion.

See supplementary material for (S1) representative snapshots of the water/nitrobenzene and water/hexane liquid/liquid simulation boxes and (S2) water finger coordinate calculation methodology.

This work is supported by the National Science Foundation through Grant No. CHE-136076. J.J.K. also acknowledges R. J. Vigna for helpful discussions and the University of California, Santa Cruz Chancellor’s Dissertation Year Fellowship for financial support.