Water Contact Angle Dependence with Hydroxyl Functional Groups on

Article pubs.acs.org/est

Water Contact Angle Dependence with Hydroxyl Functional Groups on Silica Surfaces under CO2 Sequestration Conditions Cong Chen,*,† Ning Zhang,‡ Weizhong Li,† and Yongchen Song† †

Key Laboratory of Ocean Energy Utilization and Energy Conservation of Ministry of Education, Dalian University of Technology, Dalian 116024, P. R. China ‡ School of Petroleum and Chemical Engineering, Dalian University of Technology, Panjin 124221, P. R. China S Supporting Information *

ABSTRACT: Functional groups on silica surfaces under CO2 sequestration conditions are complex due to reactions among supercritical CO2, brine and silica. Molecular dynamics simulations have been performed to investigate the effects of hydroxyl functional groups on wettability. It has been found that wettability shows a strong dependence on functional groups on silica surfaces: silanol number density, space distribution, and deprotonation/protonation degree. For neutral silica surfaces with crystalline structure (Q3, Q3/Q4, Q4), as silanol number density decreases, contact angle increases from 33.5° to 146.7° at 10.5 MPa and 318 K. When Q3 surface changes to an amorphous structure, water contact angle increases 20°. Water contact angle decreases about 12° when 9% of silanol groups on Q3 surface are deprotonated. When the deprotonation degree increases to 50%, water contact angle decreases to 0. The dependence of wettability on silica surface functional groups was used to analyze contact angle measurement ambiguity in literature. The composition of silica surfaces is complicated under CO2 sequestration conditions, the results found in this study may help to better understand wettability of CO2/brine/silica system. found to be 1.8.33 Neutral silanol terminated surfaces can only be found at the point of zero charge pHPZC (2−3).18 At pH below or above pHPZC, protonation or deprotonation process leads to net positive or negative charge at the surface. The pHdependent equilibrium of silica surface groups in aqueous solutions has stimulated a large amount of research.33−35 The type and number density of silanol groups depend on many factors, such as the cleavage plane, thermal pretreatment, humidity, pH, cleaning protocols, etc., and the deprotonationprotonation equilibrium makes the surface composition more complicated.30−34,36 Reservoir conditions vary widely from site to site (substrate type, pH, type of cations, ionic strength, temperature, pressure, and so forth).35 Under the CO2 sequestration environment, the functional groups on silica surfaces may be a combination of Q2, Q3, Q4, and their protonated/deprotonated forms. The surface composition has been proven to affect the wettability on silica surfaces (Q2 surface is hydrophilic27 and Q4 surface is hydrophobic21,37). Wettability on the silica surface with different functional groups should be very complicated. An investigation of the wettability on silica surfaces with different functional groups is helpful to

1. INTRODUCTION The wettability of CO2/water/mineral systems plays an important role when injecting supercritical CO2 into saline aquifers.1,2 Wettability on silica surfaces has been extensively studied using both experimental3−19 and molecular dynamics simulation methods.4,20−28 However, the measured values by different authors are different.4 Several possible factors that cause the uncertainty of measured contact angles have been proposed and surface contamination was suggested as one of the main factors.6,8 In fact, even without contamination, the functional groups on silica surfaces are complicated. Silica exists in various crystalline polymorphs and amorphous phases which exhibit different bulk and surface properties. The silica surface is usually hydroxylated in a water environment to form five possible compositions: Si(OH)4, Q1, Q2, Q3, and Q4.29,30 Three types of hydroxyl groups on silica surfaces are generally found: isolated, geminal, and vicinal hydroxyl groups.30−32 For amorphous silica surfaces, the surface silanol density was found to be 4.6 OH/nm2 at a typical density of isolated, geminal and vicinal silanols of 1.2, 0.6, and 2.8 OH/ nm2, respectively.30,32 At the surface of crystalline polymorphs, the presence of geminal species can be more substantial. The composition of a quartz surface reacted with various aqueous solutions of pH 0−10 was qualitatively and quantitatively evaluated using X-ray photoelectron spectroscopy and the atomic ratio between the surface oxygen and silicon atoms was © 2015 American Chemical Society

Received: Revised: Accepted: Published: 14680

July 28, 2015 October 27, 2015 October 28, 2015 October 28, 2015 DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology Table 1. Summary of Molecular Dynamics Simulation Studies on CO2/Brine/Silica Wettability in Literature reference

surface type

predicted contact angles

21

Q4

24 28 28 24

Q4 pristine silica plane with silicon atoms in the outmost layer bare Si atoms and a few protonated silanol functional groups (1.6 OH/nm2) Q3/Q4 (1.7−4.5 OH/nm2)

(1) 275 K: 20.7 MPa, 39°; 19.1 MPa, 52.4°;a (2) 300 K: 18.9 MPa, 29.9°; 17.3 MPa, 31.9°a 296 K: 7.8 MPa, 19°; 13.9 MPa, 20°; 28.8 MPa, 17° 318 K, ∼11 MPa: pure water, 22.6°; 3 M CaCl2, 35.5°; 3 M NaCl, 29.5° (1) 300 K: 10−20 MPa, ∼80°; 4.1 MPa, 1 M NaCl, 66.5°;b 4.5 MPa, 4 M NaCl, 70° (2) 350 K: 12 MPa, 66°; 17.5 MPa, 76°b 330 K, 20 MPa: water, 145°(31°); 0.78 M NaCl, 141°(26°)c 0.26 M CaCl2, 145° (34°) 300 K, 12−15 MPa, 1 M CaCl2, ∼ 80°b fluid density: 0.2 g/cm3, 115° fluid density: 0−1 g/cm3, 56−87°b same as those predicted on Q4 surface24

27

Q2

23 39

Q2 Q2

26

Q4

a Pressures were computed from three components of the pressure tensor. bThese contact angle data were extracted from graphics. cContact angles outside brackets were calculated by water droplet method and contact angles inside brackets were calculated by CO2 droplet method.

Figure 1. Snapshots of silica surfaces with different functional groups. On Q3-9% and Q3-50% surfaces, a certain number of Na+ was drawn. Si: yellow, O: red, H: white, Na: blue.

better understand the wetting characteristics of CO2/brine/ silica systems during CO2 sequestration, and it is also helpful to better understand the uncertainty of the measured contact angles4,6,8 and the dewetting process.12 However, wettability on silica surfaces with different functional groups is hard to investigate experimentally. This is because the type and concentrations of functional groups are easy to change with reaction or altered by contaminations. A former study investigated water/vapor contact angles on quartz surfaces with different functional groups.38 Quartz plates were pretreated by heating them in vacuum at temperatures from 200 to 1000 °C to generate surfaces with different functional groups, and then contact angles were measured through the drop 1 min after removal from vacuum at room temperature. Under geologic CO2 sequestration conditions, the experiments must be conducted at high temperature and pressure. After sample surface function groups were altered to the desired condition, the sample has to be placed into chamber where surface functional groups may be changed by following heating and pressurizing process in CO2 or brine environment. So, conducting high pressure, high temperature experiments with varying hydroxyl group structures are very challenging.

Alternatively, molecular dynamics (MD) simulation is a better choice in cases where high quality experiments would be very arduous and accurate force fields exist in the system of interest. MD simulation has been applied to predict contact angles on several model silica surfaces under CO2 atmosphere, as summarized in Table 1.21,23−28,39 Most MD simulations of silica surfaces were made on Q223,27,39 and Q421,24,26 surfaces. Only two former studies used silica surfaces with different silanol number densities.24,28 A surface with bare Si atoms and a few protonated silanol functional groups (1.6 OH/nm2) was constructed to predict contact angles as a function of CO2 pressures.28 However, the surface with bare Si atoms was not a typical polymorph for silica.28 A combination of Q3 and Q4 surfaces with silanol densities from 1.7 to 4.5 OH/nm2 was also used.24 But a hydrophilic characteristic for the Q4 surface was predicted which was different from other studies in literature.21,37 To the best of our knowledge, the study of water contact angles on silica surfaces with different silanol functional groups is missing, and the dependence of water contact angles on silanol functional groups is open to question. In the present study, the dependence of water (brine) contact angles on silanol functional groups was investigated under the atmosphere of supercritical CO2. The effects of the 14681

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology

There are two silica surfaces in each simulation box. When the system is at equilibrium, the bottom silica surface is in equilibrium with the droplet phase and the bulk phase. The top silica surface is in equilibrium with the bulk phase. The distances between two silica surfaces were long enough to make sure that the structures near the two silica surfaces were not affected by each other. The construction of simulation boxes makes it possible to simultaneously investigate contact angles and silica−water or silica−CO2 interactions. The contact angles were analyzed based on the zone in which the z coordination was less than a half of the total length in the z direction. The rest zone was used to analyze silica−water interactions (CO2 droplets). A total of 24 simulation boxes were constructed, and their parameters are summarized in Table 2. The diameters of the droplets were around 3.5 nm. The x-dimensions of the simulation boxes were selected to guarantee enough space for the spreading droplets. The z-dimensions were selected to make sure that the distances between the top silica surfaces with droplets were larger than 6 nm to eliminate the possible interaction of top silica surfaces with water or CO2 droplets. Along the y-direction, there were at least 500 water molecules per nm for the water droplet method and 160 CO2 molecules per nm for the CO2 droplet method. Xions is the molarity of NaCl or CaCl2. NSiOH is the number of silanol groups per nm2. Lx, Ly, and Lz are dimensions of simulation boxes in each direction, respectively. The units for Lx, Ly, and Lz are nm. Nwater and NCO2 are numbers of water and CO2 molecules in each simulation box, respectively. 2.3. Simulation Details. An open source molecular dynamics simulation package NAMD was used to conduct all simulations.40 In three dimensions, periodic boundary conditions were applied. Standard techniques such as neighborhood lists and switching function were used to compute nonbonded interactions. The cutoff distance was 13.5 Å. A smooth switching function was used to truncate the van der Waals potential energy smoothly at the cutoff distance to improve energy conservation. The switching process started from 12.0 Å. Neighborhood lists were updated every 10 time steps. Coulombic interactions were calculated using a Particle mesh Ewald (PME) method.41 A cubic PME interpolation order was applied, and the direct sum tolerance was 10−6. The PME grid sizes in three dimensions were selected to make sure the grid spaces were about 1.0 Å. Atoms in silica were held at

silanol density, space distribution, deprotonation, and protonation ratio on the water contact angles were analyzed. The uncertainty of measured water contact angles in literature4,6,8 and dewetting process12 due to CO2/silica/brine reaction were discussed.

2. SIMULATION METHODS 2.1. Silica Surfaces. Six silica surfaces were selected, namely Q3, Q3-amorph, Q3/Q4, Q4, Q3-9%, and Q3-50% surfaces.29 Q3-amorph is an amorphous Q3 surface. Q3/Q4 is a combination of Q3 and Q4 surfaces. Q3-9% and Q3-50% are two Q3 surfaces with deprotonation degrees of 9% and 50%, respectively. The configurations of these silica surfaces are illustrated in Figure 1. 2.2. System Construction. CO2 and water were saturated at the desired temperature and pressure conditions (318 K and 10.5 MPa). CO2 saturated water and water saturated CO2 were used as two fluid phases. Silica surfaces were equilibrated at the desired temperature and pressure conditions before use. Two silica surfaces were parallelly placed. Then a half-cylindrical droplet made by one phase was placed onto one surface, and the rest was filled with the other phase. Both water droplet and CO2 droplet methods were used. A typical simulation box (water droplet) is illustrated in Figure 2.

Figure 2. A snapshot of simulation box. Silica and water were drawn in VDW format (solid van der Waals spheres for atoms) and CO2 molecules were drawn in Lines format (simple lines for bonds and points for atoms). For clarity, only the outermost water molecules were shown.

Table 2. Parameters for Silica Surfaces with Different Functional Groups water droplet method surface Q3

Q3-amorph

Q3/Q4

Q4 Q3-9% Q3-50%

CO2 droplet method

Xions

NSiOH

Lx

Ly

Lz

Nwater

NCO2

Lx

Ly

Lz

Nwater

NCO2

0 3 M NaCl 1 M CaCl2 0 3 M NaCl 1 M CaCl2 0 3 M NaCl 1 M CaCl2 0 3 M NaCl 0.3 M NaCl 0.3 M NaCl

4.7 4.7 4.7 4.7 4.7 4.7 2.4 2.4 2.4 0 0 4.28 2.35

20.19 20.19 20.19 20.19 20.19 20.19 20.23 20.23 20.23 14.09 14.09 20.19 20.19

3.51 3.51 3.51 4.14 4.14 4.14 3.51 3.51 3.51 3.59 3.59 3.51 3.51

13.39 13.39 13.39 12.83 12.83 12.83 13.32 13.32 13.32 14.74 14.74 13.12 13.08

2085 2085 2085 2085 2085 2085 2085 2085 2085 2295 2295 2085 2085

4406 4406 4406 4406 4406 4406 4419 4419 4419 3252 3252 4406 4406

20.02 20.02 20.02 20.15

3.48 3.48 3.48 4.28

13.39 13.39 13.39 12.83

15826 15826 15826 19494

611 611 611 754

20.02 20.02 20.02 20.28 20.28 20.03 20.04

3.48 3.48 3.48 3.64 3.64 3.45 3.48

13.32 13.32 13.32 12.69 12.69 13.12 13.08

15826 15826 15826 16070 16070 15826 15826

611 611 611 611 611 611 611

14682

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology

Figure 3. Snapshots of final configuration of simulation boxes using water droplet method. Pure water: Q3, Q3-amorph, and Q3/Q4.

Figure 4. Snapshots of final configuration of simulation boxes using CO2 droplet method. Pure water: Q3, Q3-amorph, and Q3/Q4.

used. In NPT and NVT simulations, temperature was fixed to 318 K using a Langevin dynamics method53 with a damping coefficient of 5/ps. NPT ensembles were applied during CO2− water saturation and silica surfaces equilibration processes. All other simulations were conducted in NVT ensembles. To release interfacial stress generated during construction, in all simulation boxes, molecules in the two fluid phases were fixed alternately to conduct NVT simulations at the desired temperature. The stress-removing simulation runs usually took 2−3 ns. Then, additional 15 ns runs were performed. The first 12 ns-run was used to equilibrate, and data were collected in last 3 ns. Contact angles were predicted by tangent lines of 2D density profiles of water. The average positions for outermost silicon atoms were chosen as a baseline. To test the potential impact of baseline selection on predicted water contact angles, we changed the position of the baseline by 0.3−0.5 nm, and no significant effects on contact angles were found. The trajectories for 1 ns run were averaged to calculate density profiles. A 3 ns run produces three density profiles. As a result, for each simulation box, there were a total of six contact angle data points from which the average and error values were calculated. All data were processed in VMD software54 using inhouse codes. 2.4. Hydrogen Bond Analysis. In the MD simulations, the definitions of hydrogen bonds are somewhat arbitrary because of the lack of information on electron density.55 In the present study, the geometrical definition was selected to identify a hydrogen bond. In a hydrogen bond XH···A, the molecule in which atom X is involved in was a donor, and the other

their initial positions using the SHAKE algorithm.42 The time step was 1 fs. Nonbonded interactions were calculated at every time step, while full electrostatic evaluations were computed every two time steps. A multiple time step integration technique r-RESPA43 was adopted. Molecular dynamics simulations are sensitive to selection of force fields, and there are many interatomic potential models for silica, CO2, and water in literature.29,44−48 We selected Heinz’s force field to compute interactions between silica and water because it correctly predicted water contact angles of silica/water/vapor systems.29 Heinz’s force field was also optimized for silica surfaces with different silanol functional groups.29 The objective of the study is to investigate water contact angles dependence on silanol functional groups and we think Heinz’s force field is a good choice for this purpose. For CO2, we selected a flexible CO2 potential model,47 as it gave good predictions of CO2− water interfacial tensions.49 Accurate MD simulation predictions of the fluid−fluid interfacial tension require a longer cutoff distance or the inclusion of van der Waals interactions in the reciprocal-space Ewald sum calculation.50 An additional simulation for the Q3 surface with a longer cutoff distance (17 Å) was performed, but no significant effect on the final water contact angle is found. All simulations were performed in NPT or NVT ensembles. In NPT simulations, the system pressure was fixed to 10.5 MPa using a Langevin piston Nose-Hoover method. The method was a combination of the Nose-Hoover constant pressure method51 with the piston fluctuation control implemented using Langevin dynamics.52 A Langevin piston period of 100 time steps and a Langevin piston decay of 50 time steps were 14683

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology Table 3. Predicted Contact Angles for Different Surfaces θwater surface 3

Q

Q3-amorph

Q3/Q4

Q4 Q3-9% Q3-50%

θCO2

Xions

NSiOH

average

error

average

error

θwater + θCO2

0M 3 M NaCl 1 M CaCl2 0M 3 M NaCl 1 M CaCl2 0M 3 M NaCl 1 M CaCl2 0M 3 M NaCl 0.3 M NaCl 0.3 M NaCl

4.7 4.7 4.7 4.7 4.7 4.7 2.4 2.4 2.4 0 0 4.28 2.35

33.5° 36.4° 35.8° 53.1° 56.1° 53.6° 61.5° 72.3° 69.7° 146.7° 153.8° 21.9° ∼0°

2.36° 8.48° 5.91° 3.45° 4.51° 4.64° 7.16° 4.57° 6.32° 4.84° 6.29° 8.41° -

146.0° 145.2° 140.6° 124.6°

2.57° 6.49° 5.65° 7.48°

179.5° 181.6° 176.4° 177.7°

119.3° 110.0° 109.9° 34.1° 27.2°

9.53° 4.23° 9.48° 2.87° 4.30°

180.8° 182.3° 179.6° 180.8° 181.1°

θwater and θCO2 are predicted water and CO2 contact angles on silica surfaces by water droplet and CO2 droplet methods, respectively.

any surface studied, which is in accordance with results in other studies.21 3.2. Discussion. 3.2.1. CA vs Silanol Density. For neutral silica surfaces with crystalline structure (Q3, Q3/Q4, Q4), the water contact angle increases as silanol number density decreases. For pure water, the water contact angle increases from 33.5° to 146.7° when the silanol number density decreases from 4.7 to 0 OH/nm2. Former studies of the Q2 surface with a higher silanol number density (9.4 OH/nm2) had lower water contact angles (30°,27 22.6°,39 and 20°23). For 3 M NaCl solutions, in response to the decreasing of silanol number density, the water contact angle increases from 36.4° to 153.8°. The same trend follows for 1 M CaCl2 solutions. Similar silanol density dependence for water contact angles was also found by silica/water/vapor29 and glass/water/CO219 systems. The dependence of water contact angles with silanol density for neutral silica surfaces may be explained by hydrogen bonds between silanol groups and water. Hydrogen bonds at the silica and water interface were identified, and the hydrogen bonds number densities (mean numbers of hydrogen bonds per nm2) were analyzed. The results for silica-CO2 saturated water interactions are illustrated in Figure 5. As the silanol number density decreases, the mean number of hydrogen bonds between SiOH and water decreases leading to contact angle increase. The hydrogen bonds number density decreases from 7.7 to 2.8 HB/nm2 when the silanol number density decreases from 4.7 to 0 OH/nm2 (Q3, Q3/Q4, Q4).

molecule was an acceptor. By geometry, two lengths (X···A and H···A) and two angles (XH···A and H-X···A) are usually used to define hydrogen bonds.56,57 The length of H···A as well as the angle XH···A were selected. When the length of H···A is less than a threshold, and the angle XH···A is larger than the threshold, H···A is identified as a hydrogen bond. Thresholds of 0.3 nm and 130° were selected, which has been successfully used to simulate hydrogen bonds between water and silanol groups at quartz-water interface.56

3. RESULTS AND DISCUSSIONS 3.1. Predicted Contact Angles. Molecular dynamics simulations were produced and the final configurations for Q3, Q3-amorph, Q3/Q4 and Q3-9% surfaces are illustrated in Figures 3 and 4. The results for other silica surfaces are illustrated in the SI File (S1−S4). Density profiles were calculated and contact angles were predicted. The results are summarized in Table 3. Two typical density profiles can be found in the SI File (S5). It is clear that contact angles on silica surfaces with different functional groups are different. A hydrophilic to hydrophobic transition occurs when the silica surface changes from Q3 to Q4. Contact angles predicted using CO2 droplet method agree well with those using the water droplet method. The sums of water and CO2 contact angles are around 180° and the largest deviation is 3.6°. It should be noted that the good agreement between results obtained by the water droplet method and the CO2 droplet method is based on the fact that no gravity is considered in MD simulations. The result of contact angles for pure water on Q4 surface agrees well with Tenney and Cygan’s result.21 A contact angle of 80° on Q4 surface was found at 10− 20 MPa and 300 K.26 The reason may be that spherical droplets were used, which have been found to cause large line tension and not to reproduce the hydrophobic characteristic of the Q4 surface.21 The effects of salinity were also studied. When the ionic strength increases from 0 to 3M, the water contact angles increase about 7° on Q4 and 8−11° on Q3/Q4 surfaces. The effect of ionic strength on water contact angles on Q3 and Q3amorph surfaces is very weak. The effects of salinity on the water contact angles on the Q4 surface have also been investigated using 0.78 M NaCl and 0.26 M CaCl2.21 The effects seemed to be very insignificant for such small ionic strengths. Cation types do not affect water contact angles on

Figure 5. Water contact angles as a function of mean numbers of hydrogen bonds per nm2 for silica−CO2 saturated water interaction. 14684

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology 3.2.2. CA vs Silanol Space Distribution. The dependence of water contact angles with silanol density is challenged by the result of Q3-amorph surface. The predicted water contact angle on the Q3-amorph surface is 20° larger than that on Q3 surface which has the same silanol density. This implies that beyond the silanol density, the silanol space distribution also affects water wettabilities on silica surfaces. The mean number of hydrogen bonds decreases from 7.7 to 5.7 HB/nm2 when Q3 surface changes to an amorphous structure (Figure 5). There are no hydrogen bonds between silanol groups on the Q3 surface, however, on the Q3-amorph surface, the number density of hydrogen bonds between silanol groups increases to 0.55 HB/nm2. Hydrogen bonds between silanol groups decrease the chance of silanol groups to form hydrogen bonds with water. Properly spaced hydroxyl groups can provide a match to the structure of ice leading to strong hydrogen bonds between substrate and water molecules.37 When Q3 changes to an amorphous structure, the hydroxyl groups may be not properly spaced to form an “ice like” structure which leads to a decreasing number of hydrogen bonds. 3.2.3. CA vs Silanol Deprotonation/Protonation. For the Q3 surface, when 9% silanol groups are deprotonated, water contact angles decrease from 33.5° to 21.9°. When 50% silanol groups are deprotonated, water spreads on the surface leading to a very low contact angle close to 0. When the pH decreased from 6 to 3, protonation occurs on the silica surface, and the contact angle increased from 19° to 23° in the air−brine−silica system at ambient pressure and temperature.58 The increase (decrease) of water contact angles with protonation (deprotonation) may be explained by destabilization (stabilization) of water films on silica surfaces caused by surface charge changing.18 Dewetting of silica surfaces were found during a micromodel study.12 The micromodels were made out of fused silica plates which were amorphous. When the silica plates were saturated by degassed brine with pH of 5.8−6.2, the deprotonation ratio may be high leading to very low contact angles. When the brine was displaced by supercritical CO2, CO2 gradually dissolved into brine to make pH decrease. However, pH values at different locations were not same as fluid flow can change local pH values.59 As a result, dewetting was only found in some locations, and the final water contact angles changed with location.12 Our predicted contact angle on Q3-amorph surface for 3 M NaCl solutions is 56.1°. Contact angles of the dewetting droplets were 55°−76° for 3 M NaCl.12 3.2.4. CA vs Sequestration Conditions. The predicted contact angles on different silica surfaces show a strong dependence of contact angle on silica surface functional groups: Si(OH)2, SiOH, SiO−, SiOSi). In laboratory experiments, the density, type, and distribution of these functional groups are governed by many factors, such as mineral type, cleavage plane, thermal pretreatment, humidity, pH, cleaning protocols, cation type, ionic strength, and so forth. Different functional groups on silica surfaces may be one of the reasons for ambiguity in the contact angle measurements for silica/CO2/brine systems. The charge density of silica surfaces in water, which is controlled by the acid/base equilibrium between surface groups, vanishes for pH values in the range of 2−3.18 When saturated with CO2, pH in brine solutions is in the range of 3. In laboratory experiments, if the CO2 and brine were fully saturated, then the deprotonated ratio of silica surface silanol groups was low.12 Experimentally measured contact angles on quartz surfaces with similar conditions were in the range of 0−47°(advancing).3,5,7,8,10,11,15,60,61 However, the functional groups on silica

surfaces used in these experiments were unknown. We cannot directly compare the predicted contact angles using MD simulations with those obtained by experiment. Several experimental measured contact angles on amorphous silica were 70°17 and 45−70°(0−5 M NaCl).13 The contact angles on amorphous silica were larger than those measured on quartz, showing the effects of silanol space distribution found in this study. In real supercritical CO2 sequestration conditions, the situation becomes more complicated due to the long time mineral-brine reaction, wide sequestration conditions and fluid flow, and heat transfer characteristics. The true surface should be a combination of possible functional groups (Q2, Q3, Q4, SiO−, et al.). To better understand the wettability of CO2/ brine/silica systems, knowledge about surface functional groups is essential. Once the detailed surface composition is known, molecular dynamics simulation is a good method to predict water contact angles.

4. IMPLICATIONS Wettability of geologic minerals in the presence of CO2 and formation brines highly influence the structural and residual trapping process during CO2 sequestration. For silica minerals, in real CO2 sequestration conditions, the true surface should be a combination of possible functional groups (Q2, Q3, Q4, SiO−, et al.). The dependence of wettability on functional groups found in this study implies that the surface functional group alteration due to brine−CO2−silica reactions may have important consequences, such as CO2 residual saturation, relative permeability, and capillary entry pressure.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.5b03646. Figure S1.1, Snapshots of final configuration: 3 M NaCl on Q 3 surface; Figure S1.2, Snapshots of final configuration: 1 M CaCl2 on Q3 surface; Figure S2, Snapshots of final configuration: 3 M NaCl on Q3/Q4 surface; Figure S3, Snapshots of final configuration: pure water on Q4 surface; Figure S4, Snapshots of final configuration: 0.3 M NaCl on Q3-50% surface; Figure S5.1, A density profile for Q3-amorph surface using CO2 droplet method; Figure S5.2, A density profile for Q3amorph surface using water droplet method. (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 86-411-84708774; e-mail: [email protected] (C.C.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Special thanks to Editor Daniel Giammar and anonymous reviewers. This paper is greatly improved with their fruitful comments. The Fundamental Research Funds for the Central Universities (DUT14LAB13), the Doctoral Startup Funds of Liaoning Province (20121021) and National Natural Science Foundation of China (NSFC, 51206016) provide funding supports on this research. The authors give special thanks to Computing Center in Department of Energy and Power 14685

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology

(17) Sutjiadi-Sia, Y.; Jaeger, P.; Eggers, R. Interfacial phenomena of aqueous systems in dense carbon dioxide. J. Supercrit. Fluids 2008, 46, 272−279. (18) Chiquet, P.; Broseta, D.; Thibeau, S. Wettability alteration of caprock minerals by carbon dioxide. Geofluids 2007, 7 (2), 112−122. (19) Dickson, J. L.; Gupta, G.; Horozov, T. S.; Binks, B. P.; Johnston, K. P. Wetting phenomena at the CO2/water/glass interface. Langmuir 2006, 22 (5), 2161−2170. (20) Javanbakht, G.; Sedghi, M.; Welch, W.; Goual, L. Molecular Dynamics Simulations of CO2/Water/Quartz Interfacial Properties: Impact of CO2 Dissolution in Water. Langmuir 2015, 31 (21), 5812− 5819. (21) 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 (3), 2035−2042. (22) Yan, M. Q.; Yang, X. N.; Lu, Y. J. Wetting behavior of water droplet on solid surfaces in solvent environment: A molecular simulation study. Colloids Surf., A 2013, 429, 142−148. (23) Tsuji, S.; Liang, Y. F.; Kunieda, M.; Takahashi, S.; Matsuoka, T. Molecular Dynamics Simulations of the CO2-Water-Silica Interfacial Systems. Energy Procedia 2013, 37, 5435−5442. (24) McCaughan, J.; Iglauer, S.; Bresme, F. Molecular Dynamics Simulation of Water/CO2-quartz Interfacial Properties: Application to Subsurface Gas Injection. Energy Procedia 2013, 37 (0), 5387−5402. (25) Hamm, L. M.; Bourg, I. C.; Wallace, A. F.; Rotenberg, B. Molecular Simulation of CO2- and CO3-Brine-Mineral Systems. In Geochemistry of Geologic CO2 Sequestration; DePaolo, D. J., Cole, D. R., Navrotsky, A., Bourg, I. C., Eds.; Mineralogical Soc. Amer.: Chantilly, 2013; Vol. 77, pp 189−228. (26) Iglauer, S.; Mathew, M. S.; Bresme, F. Molecular dynamics computations of brine-CO2 interfacial tensions and brine-CO2-quartz contact angles and their effects on structural and residual trapping mechanisms in carbon geo-sequestration. J. Colloid Interface Sci. 2012, 386, 405−414. (27) Bagherzadeh, S. A.; Englezos, P.; Alavi, S.; Ripmeester, J. A. Influence of Hydrated Silica Surfaces on Interfacial Water in the Presence of Clathrate Hydrate Forming Gases. J. Phys. Chem. C 2012, 116 (47), 24907−24915. (28) LIU, S.; YANG, X.; QIN, Y. Molecular dynamics simulation of wetting behavior at CO2/water/solid interfaces. Chin. Sci. Bull. 2010, 55 (21), 2252−2257. (29) Emami, F. S.; Puddu, V.; Berry, R. J.; Varshney, V.; Patwardhan, S. V.; Perry, C. C.; Heinz, H. Force Field and a Surface Model Database for Silica to Simulate Interfacial Properties in Atomic Resolution. Chem. Mater. 2014, 26 (8), 2647−2658. (30) Zhuravlev, L. T. The surface chemistry of amorphous silica. Zhuravlev model. Colloids Surf., A 2000, 173 (1−3), 1−38. (31) Dijkstra, T. W.; Duchateau, R.; van Santen, R. A.; Meetsma, A.; Yap, G. P. A. Silsesquioxane models for geminal silica surface silanol sites. A spectroscopic investigation of different types of silanols. J. Am. Chem. Soc. 2002, 124 (33), 9856−9864. (32) Zhuravlev, L. T. CONCENTRATION OF HYDROXYLGROUPS ON THE SURFACE OF AMORPHOUS SILICAS. Langmuir 1987, 3 (3), 316−318. (33) Duval, Y.; Mielczarski, J. A.; Pokrovsky, O. S.; Mielczarski, E.; Ehrhardt, J. J. Evidence of the existence of three types of species at the quartz-aqueous solution interface at pH 0−10: XPS surface group quantification and surface complexation modeling. J. Phys. Chem. B 2002, 106 (11), 2937−2945. (34) Rimola, A.; Costa, D.; Sodupe, M.; Lambert, J. F.; Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments. Chem. Rev. 2013, 113 (6), 4216−4313. (35) Lagstrom, T.; Gmur, T. A.; Quaroni, L.; Goel, A.; Brown, M. A. Surface Vibrational Structure of Colloidal Silica and Its Direct Correlation with Surface Charge Density. Langmuir 2015, 31 (12), 3621−3626. (36) de Vos, W. M.; Cattoz, B.; Avery, M. P.; Cosgrove, T.; Prescott, S. W. Adsorption and Surfactant-Mediated Desorption of Poly-

Engineering of Dalian University of Technology for providing parallel computing environment.

■

ABBREVIATIONS MD molecular dynamics CA contact angle

■

REFERENCES

(1) Tokunaga, T. K.; Wan, J. M. Capillary Pressure and Mineral Wettability Influences on Reservoir CO2 Capacity. In Geochemistry of Geologic CO2 Sequestration; DePaolo, D. J., Cole, D. R., Navrotsky, A., Bourg, I. C., Eds.; Mineralogical Soc. Amer.: Chantilly, 2013; Vol. 77, pp 481−503. (2) Chaudhary, K.; Cardenas, M. B.; Wolfe, W. W.; Maisano, J. A.; Ketcham, R. A.; Bennett, P. C. Pore-scale trapping of supercritical CO2 and the role of grain wettability and shape. Geophys. Res. Lett. 2013, 40 (15), 3878−3882. (3) Sarmadivaleh, M.; Al-Yaseri, A. Z.; Iglauer, S. Influence of temperature and pressure on quartz-water-CO2 contact angle and CO2-water interfacial tension. J. Colloid Interface Sci. 2015, 441, 59−64. (4) Iglauer, S.; Pentland, C. H.; Busch, A. CO2 wettability of seal and reservoir rocks and the implications for carbon geo-sequestration. Water Resour. Res. 2015, 51 (1), 729−774. (5) Al-Yaseri, A.; Sarmadivaleh, M.; Saeedi, A.; Lebedev, M.; Barifcani, A.; Iglauer, S. N2 + CO2 + NaCl brine interfacial tensions and contact angles on quartz at CO2 storage site conditions in the Gippsland basin, Victoria/Australia. J. Pet. Sci. Eng. 2015, 129, 58−62. (6) Wan, J. M.; Kim, Y.; Tokunaga, T. K. Contact angle measurement ambiguity in supercritical CO2-water-mineral systems: Mica as an example. Int. J. Greenhouse Gas Control 2014, 31, 128−137. (7) Saraji, S.; Piri, M.; Goual, L. The effects of SO2 contamination, brine salinity, pressure, and temperature on dynamic contact angles and interfacial tension of supercritical CO2/brine/quartz systems. Int. J. Greenhouse Gas Control 2014, 28, 147−155. (8) Iglauer, S.; Salamah, A.; Sarmadivaleh, M.; Liu, K. Y.; Phan, C. Contamination of silica surfaces: Impact on water-CO2-quartz and glass contact angle measurements. Int. J. Greenhouse Gas Control 2014, 22, 325−328. (9) Wang, S. B.; Tao, Z. Y.; Persily, S. M.; Clarens, A. F. CO2 Adhesion on Hydrated Mineral Surfaces. Environ. Sci. Technol. 2013, 47 (20), 11858−11865. (10) Wang, S. B.; Edwards, I. M.; Clarens, A. F. Wettability Phenomena at the CO2-Brine-Mineral Interface: Implications for Geologic Carbon Sequestration. Environ. Sci. Technol. 2013, 47 (1), 234−241. (11) Saraji, S.; Goual, L.; Piri, M.; Plancher, H. Wettability of Supercritical Carbon Dioxide/Water/Quartz Systems: Simultaneous Measurement of Contact Angle and Interfacial Tension at Reservoir Conditions. Langmuir 2013, 29 (23), 6856−6866. (12) Kim, Y.; Wan, J. M.; Kneafsey, T. J.; Tokunaga, T. K. Dewetting of Silica Surfaces upon Reactions with Supercritical CO2 and Brine: Pore-Scale Studies in Micromodels. Environ. Sci. Technol. 2012, 46 (7), 4228−4235. (13) Jung, J.-W.; Wan, J. Supercritical CO2 and Ionic Strength Effects on Wettability of Silica Surfaces: Equilibrium Contact Angle Measurements. Energy Fuels 2012, 26 (9), 6053−6059. (14) Bikkina, P. K. Contact angle measurements of CO2-waterquartz/calcite systems in the perspective of carbon sequestration. Int. J. Greenhouse Gas Control 2011, 5 (5), 1259−1271. (15) Espinoza, D. N.; Santamarina, J. C. Water-CO2-mineral systems: Interfacial tension, contact angle, and diffusionImplications to CO2 geological storage. Water Resour. Res. 2010, 46, W07537. (16) Chalbaud, C.; Robin, M.; Lombard, J.-M.; Martin, F.; Egermann, P.; Bertin, H. Interfacial tension measurements and wettability evaluation for geological CO2 storage. Adv. Water Resour. 2009, 32, 98−109. 14686

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687

Article

Environmental Science & Technology

(58) Gribanova, E. V. DYNAMIC CONTACT ANGLES TEMPERATURE-DEPENDENCE AND THE INFLUENCE OF THE STATE OF THE ADSORPTION FILM. Adv. Colloid Interface Sci. 1992, 39, 235−255. (59) Lis, D.; Backus, E. H. G.; Hunger, J.; Parekh, S. H.; Bonn, M. Liquid flow along a solid surface reversibly alters interfacial chemistry. Science 2014, 344 (6188), 1138−1142. (60) Farokhpoor, R.; Bjorkvik, B. J. A.; Lindeberg, E.; Torsaeter, O. Wettability behaviour of CO2 at storage conditions. Int. J. Greenhouse Gas Control 2013, 12, 18−25. (61) Al-Yaseri, A. Z.; Lebedev, M.; Barifcani, A.; Iglauer, S. Receding and advancing (CO2 + brine + quartz) contact angles as a function of pressure, temperature, surface roughness, salt type and salinity. J. Chem. Thermodyn. 2015, DOI: 10.1016/j.jct.2015.07.031.

(vinylpyrrolidone) on Plasma- and Piranha-Cleaned Silica Surfaces. Langmuir 2014, 30 (28), 8425−8431. (37) Isaienko, O.; Borguet, E. Hydrophobicity of Hydroxylated Amorphous Fused Silica Surfaces. Langmuir 2013, 29 (25), 7885− 7895. (38) Lamb, R. N.; Furlong, D. N. CONTROLLED WETTABILITY OF QUARTZ SURFACES. J. Chem. Soc., Faraday Trans. 1 1982, 78, 61−73. (39) Chen, C.; Wan, J.; Li, W.; Song, Y. Water Contact Angles on Quartz Surfaces under Supercritical CO2 Sequestration Conditions: Experimental and Molecular Dynamics Simulation Studies. Int. J. Greenhouse Gas Control 2015, 42, 655. (40) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (41) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98 (12), 10089. (42) Ryckaert, J. P. Special geometrical constraints in the molecular dynamics of chain molecules. Mol. Phys. 1985, 55 (3), 549−556. (43) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 1992, 97 (3), 1990− 2001. (44) Kroutil, O.; Chval, Z.; Skelton, A. A.; Predota, M. Computer Simulations of Quartz (101)-Water Interface over a Range of pH Values. J. Phys. Chem. C 2015, 119 (17), 9274−9286. (45) Bourg, I. C.; Steefel, C. I. Molecular Dynamics Simulations of Water Structure and Diffusion in Silica Nanopores. J. Phys. Chem. C 2012, 116 (21), 11556−11564. (46) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. J. Phys. Chem. B 2004, 108 (4), 1255−1266. (47) Nieto-Draghi, C.; de Bruin, T.; Perez-Pellitero, J.; Avalos, J. B.; Mackie, A. D. Thermodynamic and transport properties of carbon dioxide from molecular simulation. J. Chem. Phys. 2007, 126 (6), 06450910.1063/1.2434960 (48) Vlcek, L.; Chialvo, A. A.; Cole, D. R. Optimized Unlike-Pair Interactions for Water-Carbon Dioxide Mixtures Described by the SPC/E and EPM2Models. J. Phys. Chem. B 2011, 115, 8775−8784. (49) Zhao, L.; Lin, S.; Mendenhall, J. D.; Yuet, P. K.; Blankschtein, D. Molecular Dynamics Investigation of the Various Atomic Force Contributions to the Interfacial Tension at the Supercritical CO2Water Interface. J. Phys. Chem. B 2011, 115, 6076−6087. (50) Chen, F.; Smith, P. E. Simulated surface tensions of common water models. J. Chem. Phys. 2007, 126, 22. (51) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101 (5), 4177− 4189. (52) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant pressure molecular dynamics simulation: The Langevin piston method. J. Chem. Phys. 1995, 103 (11), 4613−4621. (53) Brunger, A. T. X-PLOR, 3.1: The Howard Hugher Medical Institute and Department of Molecular Biophysics and Biochemistry; Yale University, 1992. (54) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14 (1), 33−38. (55) Chelli, R.; Procacci, P.; Cardini, G.; Califano, S. Glycerol condensed phases Part II.A molecular dynamics study of the conformational structure and hydrogen bonding. Phys. Chem. Chem. Phys. 1999, 1, 879−885. (56) Ozkanlar, A.; Kelley, M. P.; Clark, A. E. Water Organization and Dynamics on Mineral Surfaces Interrogated by Graph Theoretical Analyses of Intermolecular Chemical Networks. Minerals 2014, 4 (1), 118−129. (57) Chen, C.; Zhang, K.; Wang, P.; Li, W.; Song, Y.; Zhang, N. Ions Hydration and Hydrogen Bonds Structure in NaCl Solutions At Temperatures and Pressures for Carbon Dioxide Sequestration: the effects of solvated CO2 molecules. Mol. Phys. 2014, 112 (1), 74−84. 14687

DOI: 10.1021/acs.est.5b03646 Environ. Sci. Technol. 2015, 49, 14680−14687