Molecular Dynamics Simulations of Retrograde Condensation in

Apr 15, 2015 - Effect of Pore-Size Distribution on Phase Transition of Hydrocarbon Mixtures in Nanoporous Media. Lei Wang , Xiaolong Yin , Keith B. Ne...
0 downloads 10 Views 6MB Size
Download PDF
Article pubs.acs.org/JPCC

Molecular Dynamics Simulations of Retrograde Condensation in Narrow Oil-Wet Nanopores. William R. W. Welch* and Mohammad Piri University of Wyoming, Department of Chemical and Petroleum Engineering, 1000 University Avenue, Department 3295, Laramie, Wyoming 82071, United States S Supporting Information *

ABSTRACT: Molecular dynamics (MD) simulations were performed for a 70/30 wt % ethane/heptane mixture unconfined and confined to 60 nm oil-wet nanocapillaries with square cross sectional widths of 4 nm. Large pressure ranges along both prograde (310 K) and retrograde (365 K) isotherms were examined. For the unconfined fluid at 310 K, compression resulted in steady increases in density of both ethane and heptane up to the predicted condensation point where a sudden phase transition was observed. At 365 K, pressure increases caused increased density in ethane but reduced density in condensed heptane droplets and retrograde phase behavior was observed as a gradual collapse of the denser phase with increased pressure above 60 bar. In nanopores, surfaces with strictly repulsive walls greatly favored association with ethane at all pressures and promoted exclusion of heptane from the narrow pore at pressures inside the saturation curve. At the oil-wet surface, preferential adsorption of heptane was seen, the extent of which depended on both temperature and pressure. At 310 K, capillary condensation, primarily of heptane, was observed at all pressures spanning the saturation curve, maximizing at low pressures. At 365 K, heptane accumulation in the pore peaked at intermediate pressures. At lower pressures, two distinct phases were observed: an adsorbed phase composed largely of the heavier molecule and a fluid phase composed mainly of the lighter component. While both equilibrium adsorption and accumulation of fluid inside the pore was significantly reduced at pressures below 30 bar along the retrograde isotherm, when pressure gradients were induced across the narrow pore, significant clogging was observed at the low pressure end for pressures as low as 10 bar.

■

INTRODUCTION In recent years, demand for domestic petroleum feedstock in North America and the depletion of conventional reservoirs worldwide has prompted increased global interest in unconventional petroleum reservoirs. Tight oil formations previously considered to be too impractical and expensive for large-scale production are now considered important fuel sources. Some of the most formidable challenges in procuring petroleum from these reservoirs stem from their low permeability, and, in shale rocks, the affinity of hydrocarbons to the organic medium in which it is generated.1 It has been shown that these types of nanochannels, which can approach molecular dimensions, constitute significant proportions of the material’s porosity; accordingly, such pore networks are the key components of hydrocarbon storage and mediation of fluid flow in shale formations.2 As demand for hydrocarbon production from shale formations increases, so does the demand for understanding phase behavior and dynamics of petroleum compounds in carbonaceous nanoconfined spaces.3−8 Kerogen in shale oil formations is a compressed amalgamation of large organic molecular networks which decompose over time rendering free hydrocarbons confined to nanoscale pore systems defined by the remaining organic solid.9,10 While © XXXX American Chemical Society

production from shales is typically associated with low molecular weight gases, shale rocks can also contain considerable amounts of heavier species and, in general, complex mixtures that may include water and brine coexist as vapor, liquid, and surface-adsorbed phases.2,11 Among other things, condensation and size exclusion are known factors that complicate fluid behavior in shale pore networks, nonetheless, current understanding of the behavior of petroleum mixtures in shale nanonetworks is mostly speculative.12−14 Of particular interest for petroleum community are retrograde fluids, which exhibit anomalous phase transitions and solubility properties under high pressures and temperatures that fall inside the saturation curve of the PV phase diagram.15−17 Fluid mixtures exhibiting these properties have been known for over a century, but recently interest has increased as condensation of larger hydrocarbons in shale gas formations has generated optimism with respect to increasing production of heavier hydrocarbons as well as consternation about the proclivity for condensates to seal pores and inhibit Received: November 6, 2014 Revised: April 14, 2015

A

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Scheme 1

gas production.15,16 Retrograde condensation is highly sensitive to the concentration of fluid components and oil-wet networks with high surface area to volume ratios can drastically alter fluid compositions; consequently, adsorption and fluid phase behavior of retrograde mixtures, which depend on the migration history of the fluid, is complex.17 Given the obscure nature of deep subterranean pore networks, especially at the nanoscale, computational modeling is a naturally attractive way to examine physics involved in these types of systems. For nanochannels with significant solid surface area involving fluid particles with high Knudsen numbers, accounting for molecular interactions at interfaces is crucial. Consequently, molecular dynamics (MD) simulations suitable for examining nanoscale systems are increasingly practicable as computational resources become more readily available. Fluid dynamics and phase phenomena of nanoconfined Lennard-Jones (L-J) fluids, united atom molecular models, and to a larger extent, water, have been examined by MD,18−23,24 and more so by Monte Carlo methods.25 Surface interactions of confined methane, CO2, and some other hydrocarbons have seen some treatment by MD,26−29 however, similarly confined multicomponent petroleum fluids in nonmineral organic media still remain relatively unexplored by molecular simulation. In this study, we examine liquid/vapor and adsorption behavior of a simple model system, a mixture of 70/30 wt % ethane/heptane, the retrograde phase behavior of which was first thoroughly demonstrated by Kay in 1938.15 First, large simulations of the mixture (154 494 ethane molecules and 19 683 heptane molecules) were performed for prograde (310 K) and retrograde (365 K) isotherms in order to gain insight into the molecular behavior involved in retrograde phase transitions. Subsequently, the same fluid confined to solid, oilwet nanotubes 60 nm long with a square cross sectional width of 4 nm was examined under variable pressures induced by accelerating pistons on either side of the simulation cell as illustrated in Scheme 1. Finally the same type of system was used to simulate pressure gradients through shorter (40 nm) tubes in order to examine condensation processes relevant to extraction from organic shale.

simulations performed the Andersen thermostat (which produced comparable results). Statistics were taken from the final 9 ns of the production runs. Solid walls for the pore with 4 nm square cross sections were constructed as L-J spheres with radius of 0.5 nm in a face centered cubic (fcc) lattice 4 particles thick. To mimic a strongly oil-wet surface, interaction parameters were derived such that the 2 nm thick layer of the coarse grained surface particles produced a force field comparable to that produced by stacks of associated dodecane molecules in the OPLS-UA force field. This resulted in LJ spheres with σ of 0.47 nm having an ε of 3.2 kJ/mol which corresponds well with parameters used in similar coarse grained hydrocarbon models for dodecane.34−36 Solid particles were held in place using harmonic constraints in all directions. For simulations of the ethane/heptane fluid inside nanochannels, stochastic temperature control was employed through the Langevine dynamics integrator, which is necessary for nonequilibrium simulations, particularly when center-of-mass motion is an important factor.37−39 While this can result in nonconservation of momentum, we believe this is not an issue in our study because we do not examine extremely nonequilibrium simulations, but in the case of pressure gradients, quasi-static systems. Consequently, weak stochastic temperature coupling was used in these cases, with the inverse friction constant set at 50 ps to minimize dampening of dynamics while maintaining the desired temperature. In order to modulate pressure, which could vary from one end of a channel to another, pistons composed of 2 nm thick layers of L-J spheres in an fcc lattice were placed on two ends of the pore/fluid system, moving in the positive and negative z directions. Piston particles were restrained in the x and y directions and assigned high association energies with each other, forcing fcc layers to remain intact and allowing for acceleration in the z direction to create pressure in the system. Association energies between piston and fluid particles were set to approximate hard-sphere interactions as were energies between pistons and solid pore walls. The piston induced pressure system was used for both equilibrium and quasi-equilibrium simulations involving nanochannels. In principle, grand canonical Monte Carlo simulations are best suited to examine adsorption behavior of fluids, nevertheless, for the purposes of this study, relevant information about the effects of pressure on phase phenomena of the confined fluid could be adequately obtained by MD using a large enough fluid volume. Because a highly parallelized MD code is readily available in GROMACS, we were able to examine a large number of simulations over large ranges of pressures for long time periods using MD. Accordingly, data presented here are for confined fluids, technically calculated in the NVT ensemble, but effectively, under mobile piston pressures, in the NPT ensemble. Ensuring that large simulations of adsorption in narrow pores have converged to true equilibrium states proved tedious especially for systems with significant phase separation. For

■

METHODS All MD simulations were performed in parallel using the OPLSUA force field30 as implemented in GROMACS 4.6.3 software.31 Interaction cutoffs were set to be 1.5 nm with the function switched to begin continuously decaying to zero at 1.2 nm. In preliminary investigations, an integration time step of 3 fs, while producing higher bond energies (see Supporting Information, Figure S1), produced vapor/liquid phase behavior nearly identical to what was observed using shorter time steps, and the long time step was used for all simulations. Equilibrium simulations of bulk fluids were performed in the NPT ensemble using the velocity-Verlet integrator, the MTTK barostat,32 and the Nose−Hoover thermostat under strong coupling conditions for 45 ns.33 Initial configurations were obtained from B

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C example, initial rapid expulsion due to starting too far from equilibrium may be followed by condensation of droplets outside the pore, prohibiting the proper amount of heptane from re-entering the pore on any reasonable time scale. We were generally able to avoid said problem using the stepwise equilibration from high to low pressure and most simulations converged, such that average numbers of adsorbed particles did not change for the final 25 ns after a total of 250 ns. For the lower pressures at 310 K, equilibration times of 500 ns were required.

■

In Figure 1a, contours for heptane illustrate how a single component fluid should behave; for a pure fluid with strict phase transitions, there are two main regions, a broad region of lower densities and a narrower region, well separated, at higher density. The probability profiles are unimodal and approximately Gaussian in shape, indicating only one phase at each pressure. The gradual shift in the distribution centers between 2.5 and 10 bar is followed by a sudden increase at 12.5 bar. Accordingly, a clear distinction between vapor and liquid phases is predicted for pure heptane between 10 and 12.5 bar at 505 K. This is below the experimental value of 16 bar,15 which is not surprising considering the OPLS-UA force field was not parametrized for hydrocarbons at high pressures. In Figure 1b, profiles of ln(ρ) vs density probability for the binary ethane/heptane system exhibit some marked differences between the two isotherms. At 310 K, distributions are narrower. The high density tail above 20 bar increases slightly in both area and density as pressure increases, but the two phases are never distinct in the density plot under coexistence conditions. There is a sharp transition to the liquid state at 30 bar, similar to what is seen in the single component system. At 365 K, the range of pressures in which phases coexist is broader, consistent with experimental observation for isotherms closer to the critical temperature. Interestingly, at the higher temperature, the high density portion of the profiles, while increasing in area, decreases in density as the system is compressed. Along the 365 K isotherm, at lower pressures, the two phases cover distinct regions of density which merge with increasing pressure, eventually coalescing into a single phase at 72.5 bar. Again, phase transition pressures predicted for both isotherms are underestimated, in both cases by approximately 10 bar.15 Ethane and heptane density profiles are examined separately for each isotherm in Figure 2. At 310 K (Figure 2a), density distributions for ethane generally follow the same trends as those of entire fluid (Figure 1b), the main difference being that the high density tails do not overlap with the liquid phase. In the heptane distributions at 310 K, a sudden phase separation is evidenced by the leftward shift of the low density portion of the contour above 20 bar accompanied by an increase in the area of the high density region of the curve. For heptane, the high density portion does span the liquid density region. These data indicate that condensation in the saturation envelope is almost exclusively due to the larger molecule. Importantly, at 310 K, in both components, regions of highest density increase upon compression, in both volume and density, until the phase transition point. In contrast, at 365 K (Figure 2b), distributions for the ethane component are bimodal and for pressures above 45 bar, the high density components coincide with the uniform distribution seen at the highest pressures. Both low and high density regions in the ethane component increase in density upon compression, while for heptane, the density of the most compact regions decreases, as observed in the entire fluid (Figure 1b). At 365 K, while the densest regions appear to be composed primarily of heptane at lower pressures, as pressure is increased, condensed regions increase in ethane concentration, with a concomitant reduction in heptane density. The different behavior seen in the density distributions for the system at 310 and 365 K can be easily understood from visualizations of the MD trajectories. In Figure 3, snapshots of the simulations at the two different temperatures over relevant pressure ranges are illustrated. At 310 K, as pressure is

RESULTS

Simulations of Unconfined Fluids. In order to quantitatively examine phase separation in NPT simulations of bulk fluids, each simulation cell was divided into a number of smaller cubes, 3375 (153) for the simulations of the ethane/ heptane mixture and 525 (53) for smaller simulations of pure heptane. The populations and particle densities in each small cube were calculated for 200 frames over 9 ns trajectories at equilibrium for each pressure and corresponding frequency histograms were constructed (see Supporting Information Figure S2). Figure 1 shows selections from density histograms for pure heptane (Figure 1a) and for the 70/30 wt % ethane/ heptane mixture (Figure1b) at 310 K (upper panel) and 365 K (lower panel).

Figure 1. Profiles of ln(ρ) vs probability where ρ is particle density in OPLS united atoms/nm3 at selected pressures for (a) pure heptane at 505 K and (b) the 70/30 wt % ethane/heptane mixture at 310 K (top panel) and at 365 K (bottom panel). C

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 3. Representative snapshots of NPT simulation cells depicting condensation in the 70/30 wt % ethane(green)/heptane(blue) mixture at variable pressures along isotherms at (a) 310 and (b) 365 K. Snapshots are scaled to facilitate visualization and the length of one dimension for each cubic cell is indicated at the bottom of each snapshot.

indicate an obvious phase transition at 310 K, show gradual changes between 55 and 72.5 bar at 365 K, followed by a change in the curvature of the function at 72.5 bar. These data are consistent with the information presented in Figures 2 and 3, indicating the transition to purely supercritical fluid at 72.5 bar, underestimating the experimental value of 83 bar.15 In principle, the simulated volumes of liquid and vapor could be calculated from areas under the corresponding portions of density profiles, however because the two parts are often not distinct, and overlapping extensively in the retrograde region, it is difficult to tell exactly at which pressure these simulations predict the maximum condensation to occur at 365 K. Heptane−heptane radial distribution functions indicate maximum structure of associated heptane at 55 bar (Supporting Information, Figure S4). On the other hand, simple visualization of trajectories strongly indicates the maximum liquid volume to be at 60 bar. Equilibrium Simulations in Oil-Wet Pores. In equilibrium adsorption studies, single capillary tubes 60 nm long having square cross sections of 4 nm were initially given hard sphere potentials (no lipophilicity) and flooded with the retrograde fluid pushed by pistons with a force corresponding to 800 bar. The solid spheres of the packed tubes were then given the strongly oil-wet lipophilicity parameter (ε = 3.2 kJ/ mol), and the fluid was allowed to expand at both 310 and 365 K under pressures encompassing the saturation envelope of the bulk fluid for both temperatures. Expansion was controlled in a stepwise fashion such that simulations for each pressure were started from equilibrated conformations at the pressure immediately above in the series. For quasi-equilibrium simulations of confined fluids under pressure gradients, simulations were started from final configurations of equilibrated simulations under 90 bar of pressure on either side and then pistons at each end were given different accelerations. The initial simulated ethane and heptane distribution in pores having repulsive walls were surprising, especially in

Figure 2. Profiles of ln(ρ) vs probability of the 70/30 wt % ethane/ heptane mixture with partial densities for ethane and heptane plotted separately. Simulations were performed at (a) 310 and (b) 365 K, over the pressures indicated.

increased, total heptane association also increases up to 27.5 bar, as indicated by the increase in the area of the liquid portion of the density distributions in Figure 2, and then all at once, the droplets collapse resulting in a uniformly distributed liquid. In contrast, at 365 K, inclusion of hotter ethane into condensed droplets results in gradual expanding of the heptane droplets. Up to approximately 60 bar, the volume of the condensed phase grows; above 60 bar, the high pressure vapor phase effectively breaks apart the liquid and this happens over a broad pressure range; we take this to be the predicted retrograde region. Here, animated visualizations (see Supporting Information) show increasingly relaxed regions of heptane association that persist throughout the simulations up to 72.5 bar. For the retrograde isotherm, two pressures of interest may be approximated (within 2.5 bar) from these simulations, the dew point pressure and the pressure of maximum condensation although the latter is somewhat more difficult to ascertain. Density data (Supporting Information, Figure S3) which D

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

condensation, exclusion of heptane from the pore surface results in heptane accumulating outside the narrow pore (Supporting Information, Figure S5). For subcritical pressures, association with the surface is diminished for both species, but still more so for heptane. For the strongly oil-wet surface, the distribution of heptane is opposite to what is seen for the nonoil-wet surface. At the oilwet surface, there is a higher concentration of heptane at the surface where the distribution mimics the solid structure. For the lower temperature, outside of the first adsorption layer, heptane distribution is fairly uniform and this is generally true across all pressures. The same is true for ethane and at 310 K, the fluid exists as a single liquid phase in the narrow, oil-wet nanotube even at low pressures. At 365 K, however, a different trend is seen. As the pressure decreases toward the lower region of the saturation envelope, heptane molecules become distributed predominantly in the first and second adsorption layers. In contrast, for ethane, at pressures above 20 bar, molecules are evenly distributed outside the first adsorption layer. Figure 5 illustrates adsorption of each species as a function of pressure in the first and second adsorption layers of the oil-wet pore (Figure 5a) as well as total accumulation of each species inside the pore (Figure 5b). The stoichiometric ratio of ethane to heptane atoms (OPLS united-atoms) is approximately 2:1 in the bulk fluid and populations inside the nanopore corresponding to this composition are observed for both temperatures at high pressures. This is more or less true for the second adsorption layers, but not the case for the first adsorption layer where the ethane/heptane ratio is always lower than the fluid composition, indicating heptane is absorbed preferentially at

comparison to those in oil-wet pores. Figure 4 shows density profiles obtained from simulations of the ethane/heptane

Figure 4. Simulated density profiles for heptane and ethane at selected pressures and temperatures in narrow oil-wet and nonoil-wet nanotubes with 4 nm cross sectional widths.

mixture confined to narrow nonoil-wet pores and oil-wet pores at the pro-grade and retrograde temperatures. Heptane does not associate well with the nonoil-wet surface; in fact heptane molecules prefer to be separated from the wall by a layer of ethane at all pressures. Under high pressures, ethane molecules pack tightly such that their distribution indicates the molecular structure of the pore surface (a lattice in this case) and to some extent this is even true for the second layer of ethane. In general, for the nonoil-wet walls, under conditions promoting

Figure 5. Simulated adsorption isotherms for a fixed amount of the 70/30 wt % ethane/heptane mixture inside an oil-wet nanopore with cross sectional width of 4 nm and length of 60 nm. Average adsorption was measured for the layer directly adsorbed to the solid lattice, within 0.6 nm, (L1) as well 0.6 nm further (L2) as indicated by the schematic in the upper left corner of panel a. Simulations were performed between 5 and 30 bar in increments of 5 bar at 310 K with additional points at 40, 50, and 70 bar. For the 365 K isotherm, simulated pressures ranged between 20 and 150 bar in increments of 10 bar with an additional point at 200 bar. Panel b shows plots of the total number of atoms (OPLS united atoms) in the pore for each species at each pressure. E

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C the oil-wet surface. At 310 K, simulated heptane adsorption in the first adsorption layer exceeds ethane adsorption up to 40 bar, and the same is true below approximately 35 bar for heptane accumulation in the rest of the pore space. In these simulations, as pressure is decreased, the concomitant volume expansion is mostly due to ethane, resulting in capillary condensation primarily of heptane at low pressures. Accordingly, there is significant liquid content in the pore at all pressures for the prograde isotherm. At 365 K, heptane accumulation at the surfaces and in the pore in general, is significantly reduced at low pressures, but reaches a maximum at 60 bar. Directly at the surface, heptane adsorption is favored over a broad range of pressures, from 40 to 110 bar and maximizes between 50 and 60 bar. The same general pattern is seen at the second adsorption layer and throughout the pore, but there is a marked difference between the first and second layers. In contrast to what is observed at 310 K for low pressures, below 30 bar, fewer heptane atoms are adsorbed at the first layer, although accounting for stoichiometry, heptane is still favored. This is not the case outside the first layer, where ethane and heptane populations more closely resemble the bulk concentration at low pressures. It is notable that upon expansion from 60 to 50 bar, at 365 K, adsorption of heptane in the first layer remains constant, at the maximum, while the heptane population in the second layer is significantly reduced. As pressure is reduced further to 30 bar, the heptane concentration in the pore, as well as in the second adsorption layer, decreases dramatically while the change in ethane concentration is much less steep. This illustrates a key difference for the confined fluid at the two different temperatures: at 365 K, ethane is not preferentially expelled from the pore at lower pressures, at least, not to the degree that it is at 310 K. For the retrograde temperature, at the point of maximum heptane accumulation, the total amount of fluid in the pore is significantly less than the maximum observed 310 K and capillary condensation is not observed for the lower pressures at 365 K. Quasi-Equilibrium Simulations in Pores under Pressure Gradients. During extraction in shale reservoirs, induced pressure gradients can cause condensation to occur as pressure decreases, and it is this condensation that is of interest. Unsurprisingly, the data above suggest substantially more liquid formation in the small nanopore at 310 K than what is seen at low pressures for the bulk fluid. However, it is less clear, considering the apparent adsorption and fluid phase separation of the heavy and light molecules, what amount of condensation should be expected at the retrograde temperature. In order to examine condensation due to pressure drops in the retrograde fluid, simulations were performed on the fluid at 365 K with different pressures at the two ends of the simulation cell. The important trends in heptane accumulation could be readily seen after 18 ns and these are illustrated in Figure 6. For the oil-wet pore, at pressures as low as 10 bar, a notable accumulation of heptane occurs at the low pressure end of the tube and this occurs to a significant extent up to 50 bar. Similarly to what is seen in the bulk, for pressures in the retrograde region, the increased disruption of condensation occurs as pressure increases and for the supercritical pressure (shown at 110 bar), there appears to be no heptane condensation at the pore junction. These clogging effects are examined more quantitatively in Figure 7, where average numbers of heptane and ethane atoms in 1 nm cross sections of the simulation cell are plotted for a

Figure 6. Truncated snapshots of simulations of the 70/30 wt % ethane (green)/heptane (blue) mixture at 365 K pushed through nanopores by pressure gradients. Rotated images on the left-hand side are zoomed in slightly so heptane accumulation inside the pores can be more easily viewed.

Figure 7. Plot of the number of atoms in 1 nm slices along the long pore axis of 40 nm pores in a selection of the pressure gradient simulations illustrated in Figure 6. Counts of heptane atoms are indicated by solid lines and those for ethane by dashed lines. The enlarged section of the plot at the top of the figure illustrates the effects of increasing pressure on heptane accumulation at the low pressure end of the nanotube.

section of the cell containing the pore in the center. The tall spike on the low pressure end of the pore indicates the density of heptane atoms at the pore entrance. Peak heptane accumulation in this selection of pressures occurs at 30 bar and it is notable, again, that at 10 bar, the pore entrance is congested to a significant extent.

■

DISCUSSION AND CONCLUSION The parameter set used here (OPLS-UA) underestimates phase transition points at the high pressures examined in this study and in order to more accurately describe the phase diagram using molecular simulations, a more robust force field would be required. Indeed, all atom models including explicit hydrogen atoms and partial charges are known to produce superior results. Additionally, for a more detailed description of the phase diagram of the bulk fluids, a Monte Carlo method may be more appropriate and histogram reweighting might be employed in order to access information about other temperatures. F

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

area to induce condensation at lower pressures than what is seen in the bulk is apparent at 310 K, but the same is not predicted at equilibrium for low pressures at 365 K in this pore. Sealing due to condensation under pressure gradients at 365 K in the narrow oil-wet pore may be alleviated under retrograde conditions, but probably not at low pressures. The models presented here offer some valuable insights into retrograde hydrocarbon fluids in narrow oil-wet pore spaces. For applications to more realistic systems, such as shale, the model system should be improved to more accurately predict phase diagrams and extended to account for more complex fluids as well as diverse pore shapes, sizes and compositions.

Nonetheless, for the highly lipophilic or systems examined in this study, in which Columbic influences are minimal, the simple van der Waals type force field provides valuable insights. MD simulations using the Lennard-Jones force field do appear to capture the phenomenon of retrograde phase behavior for this binary fluid. For compression along the retrograde isotherm, the gradual phase transition in the retrograde region is seen as increased disruption of condensed heptane droplets by the high pressure ethane phase resulting in loosely associated heptane clouds which persist over long periods of time between 62.5 and 72.5 bar. The expansion of the heptane phase under compression at 365 K is in stark contrast to what is observed for compression of the same fluid at 310 K, in which liquid density for both components increases with pressure until the system suddenly transforms to a uniform liquid with a sharp phase transition. The way in which the liquid phase collapses in simulations of the retrograde region presented difficulties in examining liquid content based on local densities. It should be pointed out, however that the type of collapsing behavior observed here might not be observed with any retrograde fluid. For more complex fluids with a broad range of molecular weights, the density of the heaviest components need not decrease and a reduction in droplet size might be observed during compression rather than uniform expansion of the entire droplet. Inside the narrow nanopore, nonoil-wet surfaces exhibit a high preference for the lighter molecule. This was observed at both temperatures and all pressures and the nature of the effect is not immediately clear. If the effect is purely entropic, it might be observed to some extent for adsorption in the oil-wet pore as well. However, no preference for ethane at the oil-wet surface greater than what can be predicted based on stoichiometric composition under any conditions is observed. In any case, there may be an enthalpic component to the effect: heptane molecules, having the stronger L-J potential field, should preferentially associate with ethane rather than the repulsive wall; consequently, ethane which is lighter and abundant surrounds heptane molecules or aggregates resulting in an ethane layer at the solid surface. In the oil-wet pore, heptane is preferentially adsorbed to the pore wall under nearly all conditions studied. At 310 K, accumulation of heptane at the pore wall under low pressures is accompanied by a significant increase in heptane concentration throughout the pore and in general capillary condensation of heptane is observed at pressures well below the bulk condensation point. It is apparent from the equilibrium adsorption simulations that heptane accumulation is maximized at approximately the same pressure as what is seen in the bulk fluid under retrograde conditions. On the basis of information obtained from equilibrium isotherms, we may expect reduction in pressure below 50 bar to open up the interior of the pore at 365 K. On the contrary, simulations of the fluid flowing under pressure gradients seem to indicate otherwise; fluid entering a 40 nm long pore under supercritical conditions at one end condenses heptane at the opposite, low pressure end for pressures as low as 10 bar. In these simulations, the maximum accumulation appears to occur at 30 bar, but above 70 bar, such accumulation is significantly reduced. In summary, molecular level interpretations for the differences between retrograde and prograde phase behavior for petroleum fluids in bulk and in carbonaceous nanopores are accessible using simple Lennard-Jones models of a binary ethane/heptane mixture. The ability for high oil-wet surface

■

ASSOCIATED CONTENT

S Supporting Information *

Figues showing effects of time step on bond energies, predicted particle distributions of unconfined fluids in NPT simulations, predicted densities for the 70/30 wt % ethane/heptane system, heptane/heptane radial distribution functions, and truncated snapshots of 70/30 wt % ethane/heptane fluid pressurized in nanopores with non-oil-wet solid surfaces and a movie showing MD animation. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*(W.R.W.W.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge financial support of Hess Corporation and the School of Energy Resources and the College of Engineering and Applied Science at the University of Wyoming. We also thank the Advanced Research Computing Center at the University of Wyoming for extensive computational resources which made this research possible.

■

REFERENCES

(1) Ross, D. J. K.; Bustin, R. M. Shale gas potential of the Lower Jurassic Gordondale Member, northeastern British Columbia, Canada. Bull. Can. Pet. Geol. 2007, 55, 51−75. (2) Clarkson, C. R.; Solano, N.; Bustin, R. M.; Bustin, A. M. M.; Chalmers, G. R. L.; He, L.; Melnichenko, Y. B.; Radliński, A. P.; Blach, T. P. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel 2013, 103, 606−616. (3) Wu, K.; Li, X. A New Method to Predict Water Breakthrough Time in an Edge Water Condensate Gas Reservoir Considering Retrograde Condensation. Pet. Sci. Technol. 2013, 31, 1738−1743. (4) Kamari, E.; Shadizadeh, S. R. An Experimental Phase Diagram of a Gas Condensate Reservoir. Pet. Sci. Technol. 2012, 30, 2114−2121. (5) Curtis, M. E.; Sondergeld, C. H.; Ambrose, R. J.; Rai, C. S. Microstructural investigation of gas shales in two and three dimensions using nanometer-scale resolution imaging. AAPG Bull. 2012, 96, 665− 677. (6) Sondergeld, C. H.; Ambrose, R. J.; Rai, C. S.; Moncrieff, J. MicroStructural Studies of Gas Shales; Society of Petroleum Engineers: Richardson, TX, 2010. (7) Nelson, P. H. Pore-throat sizes in sandstones, tight sandstones, and shales. AAPG Bull. 2009, 93, 329−340. (8) Wang, F. P.; Reed, R. M. Pore Networks and Fluid Flow in Gas Shales; Society of Petroleum Engineers: Richardson, TX, 2009. G

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (9) Pepper, A. S.; Corvi, P. J. Simple kinetic models of petroleum formation. Part I: oil and gas generation from kerogen. Mar. Pet. Geol. 1995, 12, 291−319. (10) Goth, K.; de Leeuw, J. W.; Püttmann, W.; Tegelaar, E. W. Origin of Messel Oil Shale kerogen. Nature 1988, 336, 759−761. (11) Schieber, J. Common Themes in the Formation and Preservation of Intrinsic Porosity in Shales and Mudstones - Illustrated with Examples Across the Phanerozoic; Society of Petroleum Engineers: Richardson, TX, 2010. (12) Keller, L. M.; Holzer, L.; Wepf, R.; Gasser, P. 3D geometry and topology of pore pathways in Opalinus clay: Implications for mass transport. Appl. Clay Sci. 2011, 52, 85−95. (13) Chalmers, G. R. L.; Bustin, R. M. The organic matter distribution and methane capacity of the Lower Cretaceous strata of Northeastern British Columbia, Canada. Int. J. Coal Geol. 2007, 70, 223−239. (14) Casanova, F.; Chiang, C. E.; Li, C.-P.; Roshchin, I. V.; Ruminski, A. M.; Sailor, M. J.; Schuller, I. K. Gas adsorption and capillary condensation in nanoporous alumina films. Nanotechnology 2008, 19, 315709. (15) Kay, W. B. Liquid-Vapor Phase Equilibrium Relations in the Ethane-n-Heptane Syst. Ind. Eng. Chem. 1938, 30, 459−65. (16) Bennion, D. B.; Thomas, F. B.; Schulmeister, B. Retrograde Condensate Dropout Phenomena in Rich Gas Reservoirs-Impact on Recoverable Reserves, Permeability, Diagnosis, and Stimulation Techniques. J. Can. Pet. Technol. 2001, 40. (17) Park, S. J.; Kwak, T. Y.; Mansoori, G. A. Statistical mechanical description of supercritical fluid extraction and retrograde condensation. Int. J. Thermophys. 1987, 8, 449−471. (18) Ambrose, R. J.; Hartman, R. C.; Diaz Campos, M.; Akkutlu, I. Y.; Sondergeld, C. New Pore-scale Considerations for Shale Gas in Place Calculations; In Society of Petroleum Engineers: Richardson, TX, 2010. (19) Vela, S.; Huarte-Larrañaga, F. A molecular dynamics simulation of methane adsorption in single walled carbon nanotube bundles. Carbon 2011, 49, 4544−4553. (20) Coasne, B.; Alba-Simionesco, C.; Audonnet, F.; Dosseh, G.; Gubbins, K. E. Adsorption and Structure of Benzene on Silica Surfaces and in Nanopores. Langmuir 2009, 25, 10648−10659. (21) Vela, S.; Huarte-Larrañaga, F. A molecular dynamics simulation of methane adsorption in single walled carbon nanotube bundles. Carbon 2011, 49, 4544−4553. (22) Russo, J.; Horbach, J.; Sciortino, F.; Succi, S. Nanoflows through disordered media: A joint lattice Boltzmann and molecular dynamics investigation. EPL 2010, 89, 44001. (23) Chu, Y.; Lu, J. F.; Lu, W. Q. Molecular Dynamics Simulation of Fluid Transport in Nanoscale Pore of Porous Medium. In AIP Conference Proceedings; AIP Publishing: 2010; Vol. 1207, pp 895−00. Cui, S. T.; Cummings, P. T.; Cochran, H. D. Molecular simulation of the transition from liquidlike to solidlike behavior in complex fluids confined to nanoscale gaps. J. Chem. Phys. 2001, 114, 7189−7195. (24) Cummings, P. T.; Docherty, H.; Iacovella, C. R.; Singh, J. K. Phase transitions in nanoconfined fluids: The evidence from simulation and theory. AIChE J. 2010, 56, 842−848. (25) Singh, J. K.; Docherty, H.; Cummings,P. T. Phase Transition under ConfinementComputational Nanoscience; RSC: London, 2011; Chapter 4. (26) Coasne, B.; Long, Y.; Gubbins, K. E. Pressure effects in confined nanophases. Mol. Simul. 2014, 40, 721−730. (27) Cracknell, R. F.; Nicholson, D.; Gubbins, K. E. Molecular dynamics study of the self-diffusion of supercritical methane in slitshaped graphitic micropores. J. Chem. Soc., Faraday Trans. 1995, 91, 1377. (28) Cole, D. R.; Chialvo, A. A.; Rother, G.; Vlcek, L.; Cummings, P. T. Supercritical fluid behavior at nanoscale interfaces: Implications for CO2 sequestration in geologic formations. Philos. Mag. 2010, 90, 2339−2363.

(29) Cole, D. R.; Salim, O.; Phan, A.; Rother, G.; Striolo, A.; Vlcek, L. Carbon-bearing fluids at nanoscale interfaces. Procedia Earth Planetary Sci. 2013, 7, 175−178. (30) Jorgensen, W.; Madura, J.; Swenson, C. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (31) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (32) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit reversible integrators for extended systems dynamics. Mol. Phys. 1996, 87, 1117−1157. (33) Hoover, W. Canonical Dynamics - Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (34) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S.-J. The MARTINI coarse-grained force field: Extension to proteins. J. Chem. Theory Comput. 2008, 4, 819− 834. (35) Sedghi, M.; Piri, M. Molecular Dynamics of Wetting Layer Formation and Forced Water Invasion in Angular Nanopores with Mixed Wettability. J. Chem. Phys. 2014, 141, 194703. (36) Stukan, M. R.; Ligneul, P.; Crawshaw, J. P.; Boek, E. S. Spontaneous Imbibition in Nanopores of Different Roughness and Wettability. Langmuir 2010, 26, 13342−13352. (37) Toton, D.; Lorenz, C. D.; Rompotis, N.; Martsinovich, N.; Kantorovich, L. Temperature control in molecular dynamic simulations of non-equilibrium processes. J. Phys.: Condens. Matter 2010, 22, 074205. (38) Hayashi, K.; Shiraishi, T.; Toyoda, K.; Tanaka, F.; Mori, T.; Hata, T. Temperature-controlled molecular dynamics study on velocity-dependent threshold behavior of dynamic nano-friction. Comput. Phys. Commun. 2011, 182, 2032−2035. (39) Detcheverry, F.; Bocquet, L. Thermal fluctuations of hydrodynamic flows in nanochannels. Phys. Rev. E 2013, 88, 012106.

H

DOI: 10.1021/jp511125e J. Phys. Chem. C XXXX, XXX, XXX−XXX