Molecular dynamics for ion-tuned wettability in oil / brine / rock systems

The dependence of wettability on brine ionic composition in oil/brine/rock systems, which is denoted as ion-tuned wettability, has important applications in geoscience. Due to the involvement of non-continuum effects, molecular dynamics (MD) simulation is necessary to improve understanding of its mechanism. This work establishes a reliable molecular dynamics (MD) framework to study ion-tuned wettability. We prove that our model system can well represent the wettability of the real oil/brine/rock system, and that the wettability alteration from the MD results is qualitatively consistent with the EDL repulsion theory when ion-binding does not exist. In the process to establish the MD framework, our work suggests that adding counter ions to balance interface charges is good for studying concentration effect on wettability, and the contact angle defined from mass center coordinates is effective to measure the wettability of the nano-scale MD system. This work provides the basis for on-going studies using MD s...

Up to now, the mechanism of ion-tuned wettability is still not clear.The complexity comes from the variety of ion-related interactions in OBR systems, including electrical double layer (EDL) repulsion, [16][17][18] ion-binding, [19][20][21] and acid/base interaction. 19,20EDL repulsion essentially comes from the osmotic pressure produced by EDL overlap between similarly charged oil and rock.Ion-binding is the connection by multivalent ions between similarly charged groups of oil and rock, while acid/base interaction is the Coulombic attraction between oppositely charged groups.The influence of ionic composition on ion-binding and acid/base interaction is usually called multicomponent ion exchange (MIE).It is still under debate whether EDL repulsion [22][23][24] or MIE 11,21,25 dominates ion-tuned wettability.
Ion-tuned wettability is hard to be fully understood by the classical theories and current experiments.For example, the continuum theory can obtain the contribution of EDL repulsion to wettability alteration quantitatively, but that of MIE is not considered.Actually the EDL repulsion theory itself is still far from perfect to predict wettability, as it artificially assumes the water film thickness down

Interaction functions
The interactions between atoms in the MD simulation include bonds, angles, dihedrals, Lenard-Jones and Coulombic interactions.
For bonds and angles, the harmonic form of interaction energy is used, where K b and K a are the harmonic constants, r 0 and ϕ 0 the equilibrium bond length and angle in gas phase, r the distance between bonded atoms, and ϕ the angle formed by three adjacent bonded atoms.For dihedrals, the OPLS form of interaction energy 27 is used, where are the energy constants and φ is the dihedral formed by four adjacent bonded atoms in a chain.Lenard-Jones (LJ) potential is a typical description of van der Waals force, with the form of where ε ij and σ ij are the LJ potential parameters, r ij the distance between two non-bonded atoms, and r LJ ,c the cutoff distance to calculate E LJ ,ij .For a homogeneous system, tail correction is usually used to include the contribution of long-range LJ potential (r ij > r LJ ,c ) to pressure and energy. 28owever, since the system in this work is highly heterogeneous with multiphase interfaces and EDL, we use a relatively large cutoff r LJ ,c = 1.5 nm instead of the tail correction.To calculate the mixing LJ parameters, we assume The short-range Coulombic force between atoms with charges q i and q j has the following interaction function, where 0 is the vacuum permittivity.The PPPM (Particle-particle-particle mesh) scheme [29][30][31] is used to calculate the long-range Coulomb force with a precision of 10 4 , and the slab boundary 32 is used for the non-periodic direction (z direction in Figure 1).We set r COUL ,c = r LJ ,c for a best calculation performance.

Force field parameters
Force field parameters determine the properties of MD system, and must be carefully selected to reproduce the properties that we care.We will introduce the selection of force field parameters of brine, oil, and rock respectively in this section.
We choose a flexible three-point model, SPC/Fw, 33 to model water, mainly for two reasons.First, SPC/Fw is a flexible model which can well describe the polarization of water in a highly heterogeneous system. 33In addition, SPC/Fw may match several properties of water well, as shown in Table I, among which permittivity is a very important electrical property that influences EDL structure.At 1 atm and 300 K, the relative permittivity predicted by SPC/Fw ( r = 82.3) is very close to the experimental data ( r = 78.5),similar to SPC/E model 34 ( r = 72.0)and much better than SPC model 35 ( r = 63.4),TIP3P model 36 ( r = 99.4) and flexible or non-flexible modified TIP3P model in CHARMM force field 37 ( r = 160.7 and r = 104.3).The force field for ions comes from Ref. 38.The calculation of water properties for different water models is described in detail in the supplementary material.
We apply all-atom OPLS force field 27 (the version released with BOSS 4.8) for oil, modeling decane (CH 3 (CH 2 ) 8 CH 3 ) and decanoic acid (CH 3 (CH 2 ) 8 COOH).The surface tension of decane/water interface from MD results with the above force field is close to that from experiments, as shown in Table I. 33,39,40 Simulation details to get the surface tension are in the supplementary material.
We use α-quartz to represent the rock surface, adopting the CHARMM water contact angle model. 41,42This model charges α-quartz by deprotonation of some surface hydroxyls, consistent with electrokinetics theories such as Basic Stern model 43 and triple layer model. 44,45It also fully defines harmonic bonds and angles, so flexible surface can be produced.

B. Model system
With the above basic stones, we construct the oil/brine/α-quartz system for MD simulation, as shown in Figure 1.We will introduce the geometry, interface charging and ion assignment of this system in this section.

Geometry
In our simulations, we construct five layers of α-quartz lattice, with crystal orientation [1 0 0] parallel with x axis and [0 0 1] parallel with y axis.Hydrogens or hydroxyls are added to the bare oxygens or silicons on the crystal surface.In this way, all the surface hydroxyls are existed in germinal silanol groups.The upper most oxygens of α-quartz lattice are set at z = 0.2003 nm, whose absolute value is close to the sum of van der Waals radius of H of a hydroxyl (0.11 nm) 46 and the bond length of a hydroxyl (0.096 nm), 42 so that the brine/α-quartz interface is approximately located at z = 0.
Atoms of α-quartz satisfying z > 0.6 nm are constricted by flexible bonds and angles, while the other atoms of α-quartz are fixed still to give a stable base.
The y-direction thickness of the system W is fixed as 3.2412 nm.The oil column has totally 263 CH 3 (CH 2 ) 8 CH 3 and 20 CH 3 (CH 2 ) 8 COOH.These oil molecules can form a cylinder with radius R 0 ≈ 3 nm and length W, assuming the homogeneous oil density at 0.728 g cm 3 .Initially the water molecules are located in the volume V brine = LH − πR 2 0 W , and the top vacuum with height H = 1 nm is reserved to release the internal pressure.
Here we give two sizes of the system.For size S, L = 14.7300 nm, H = 12.0000 nm, with totally 16249 water molecules.For size L, L = 20.0622nm, H = 15.0000nm, with totally 30756 water molecules.The choosing of system size for different cases will be described in Sec. 3.

Interface charging
To charge the α-quartz surface, we artificially remove hydrogens from some germinal silanol groups and then redistribute the site charges according to Ref. 41, 42.We consider two surface charge densities of α-quartz: σ 0,qz = 0 C m 2 and σ 0,qz = 0.067 C m 2 .Only the region satisfying x < (L 4.91)/2 nm or x > (L + 4.91)/2 nm can be charged, as shown in Figure 2. Since in the results θ is in the range 35-60 degrees, approximately only region of the α-quartz contacting with brine can be charged.This charge assignment mimics the condition when oil is first attached onto the solid and brine enters later.It may differ from the real conditions since no reaction equilibrium is considered, but it is still reasonable to use it to get some preliminary results.Oil/brine interface is not charged in this work, to avoid possible ion-bindings or acid/base interaction, since we want to compare results with EDL repulsion theory.

Ion assignment
In this work, we consider the KCl solution, with similar size and diffusion coefficient for cation and anion, as the brine phase with expected bulk concentrations 0.05 M, 0.2 M and 1 M at 300 K, and 0.85 M at 373.15 K.We denote c set as the expected bulk ion concentration.
We assign ions using the rule where N K and N Cl are the number of ions in the system, and Q is total charge number of the brine.The reason to choose this ion assignment method and whether the real bulk concentration can be controlled near c set will be fully discussed in Sec.III.System size and ion assignment should be considered together.The system should be large enough to (1) avoid EDL overlap between brine/oil interface and its periodic images or brine/vapor interface, and (2) accommodate enough ions to control bulk concentration.Size L could satisfy the first condition for c set > 0.05 M, while size S for c set > 0.2 M. When the first condition is satisfied, the second is also satisfied using our ion assignment method, as discussed in Sec.III.We do not use size smaller than S, considering the finite thickness of interface.System size is set the same for charged and uncharged cases with the same c set .The size and ion numbers for all the simulated cases are summarized in Table II.

C. Running procedures
We use the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) package developed by Sandia laboratory for the classical MD simulation.We set different running procedures for charged and uncharged cases.For Cases 1-4 (σ 0,qz = 0 C m 2 ), initially oil molecules are arrayed in a √ πR 0 × √ πR 0 × W box attached onto z = 0 plane, and ions are randomly distributed in water molecules which are arrayed in the volume of L × W × H not occupied by oil.Energy minimization is performed first, then 0.1 ns nvt for oil to fully remove initial strain of oil, and finally 50 ns nvt for the whole system.We output the topology of the system when it has reached equilibrium (here we take t = 30.1 ns), and then charge α-quartz and add cations to the top vacuum to balance system charge, to give the initial state of Cases 5-7 (σ 0,qz = 0.067 C m 2 ).For Cases 5-7, we simply simulate 50 ns nvt.Temperature is kept 373.15K for Case 4 and 300 K for all the other cases.For all the cases, time step is 1 fs, the lower part of α-quartz (z < 0.6 nm) is fixed and does not contribute to the temperature of the system, and data of the last 40 ns is used for contact angle calculation.The MD results used by Sec.III-V all belong to the cases listed in Table II.

A. Difficulties and challenges
Bulk ion concentration is the concentration where ion distribution is not influenced by interface charge in an equilibrium system.According to the expected bulk concentration c set , we determine ion numbers of the system by an ion assignment method, while the real bulk concentration c real is unknown before simulation.c real may differ from c set , and the reasons are discussed in the following.We expect the difference is as small as possible to study concentration effect on wettability.Some symbols used here have been proposed in Sec.II.B.
When Q = 0, a natural way to assign ions is N K = N Cl = c set V brine .Because of the finite thickness of oil/brine interface and brine/vapor interface, ions are not totally homogeneously distributed in V brine .Therefore, c real may slightly differ from c set .
When Q 0, charged interface attracts co-ions and repel counter-ions, making the relation between ion numbers and bulk concentration vaguer.Thus the choosing of ion assignment method should be careful.Here we assume Q > 0. Three possible strategies are: (1) add cations to balance interface charge, N K = c set V brine + Q, N Cl = c set V brine ; (2) add cations and delete anions to respectively balance half of interface charge, It can be roughly judged that the first strategy makes c real > c set while the third one make c real < c set .The concrete difference between c real and c set for these three strategies can be evaluated by EDL theory.
In summary, the bulk concentration control is an important step to study ion-tuned wettability by MD simulation.In the following, we will first use the classical EDL theory to prove that the first ion assignment method is relatively good for bulk concentration control; then we will introduce how to estimate c real and evaluate the ion assignment method from MD results.

B. Theoretical analysis
In this part, we use the classical EDL theory to evaluate the three ion assignment methods as mentioned in Sec. A. Since N K N Cl = Q, N K + N Cl alone can determine the ion assignment method, When the system reaches equilibrium, the ion distribution in the diffuse layer should follow the Poisson-Boltzmann (PB) equation: where ψ is the electrical potential, e the electron charge, N A the Avogadro constant, and k b the Boltzmann constant.To simplify the calculation, instead of what we simulate as shown in Figure 1, here we analyze the brine film confined between two very large parallel plates with area S p ; one plate is charged with slip plane at x = 0 and equivalent surface charge density σ d (the opposite of surface charge of diffuse layer), while the other is uncharged and locates at x = H p .Thus V brine = S p H p .No oil is considered.The boundary conditions of Eq.( 7) for this system are given according to the Gauss Law, Thus the PB equation can be solved to get the profile of ψ.Next, we obtain Combining Eq.( 6) and Eq.( 9), we get c real = f 2 −1 ( f 1 (c set )), whose numerical solution is shown in Figure 3.We can see that ion assignment method I exhibits the best performance to control the bulk concentration, because (1) for low concentrations, only method I could control c real at the same magnitude of c set ; (2) for high concentrations, the ability of method I to control bulk concentration could be effectively improved by enlarging H p .

C. MD results: c real estimation and ion-assignment method validation
We need to estimate c real from the MD results, in order to validate the ion assignment method and to evaluate concentration effect more accurately.At the equilibrium state, local ion concentrations satisfy the Boltzmann distribution, Thus we have 47 To estimate the bulk concentration from MD results, we do the average: where ∆t means time average in time period ∆t, Ω means space average in space Ω.Ω is the brine volume satisfying 1 for different ion assignment methods when (1) H p = 6 nm and (3) H p = 12 nm.The dash line, dot line, and dash-dot line respectively represent ion assignment method I, II, and III.The black solid line represents c real = c set as a reference.Essential parameters for calculation are: T = 300 K, r = 78.5.

FIG. 4. MD results
showing the effectiveness of ion assignment method I to control bulk concentration.Parameters of the simulated cases are summarized in Table II.
(unit: nm), where R c , x c and z c are curvature radius and coordinates of curvature center of oil drop, and x c could range from x c, min to x c, max since oil drop may gently move in x direction.The smallest element of Ω is 0.2 × W × 0.2(unit: nm).
We use the MD results to see the effectiveness of ion assignment method I to control the bulk concentration, as shown in Figure 4.It is clear that the bulk ion concentration can be controlled very well.

IV. CONTACT ANGLE MEASUREMENTS
For nanodrops, accurate measurement of contact angle is a big challenge.For macroscopic drops, contact angle is usually measured by analyzing drop-surrounding interface shape, such as fitting the whole interface with Young-Laplace drop shape equation or fitting part of the interface near contact line with polynomial.However, at nano-scale, the interface shows finite thickness and irregular near the contact line, leading to the difficulties to determine the contact angle from interface shape.In Sec.A, we derive the contact angle from mass center coordinates, and in Sec.B and Sec.C we evaluate its fluctuations and representativeness of wettability from the MD results.

A. Method
The height of the mass center z m of the oil drop is a possible evaluation of the wettability of the system shown in Figure 1, and its larger value means more water-wet conditions.z m can be easily calculated in MD simulations.To present a more intuitional view of wettability, we derive the waterphase contact angle θ from z m .Here we assume a homogeneous density of oil, 0.728 g cm 3 at 300 K and 0.624 g cm 3 at 373.13 K, so that the cross section of the drop in xz plane is πR 2 0 with R 0 = 3 nm at 300 K and R 0 =3.141 nm at 373.15 K.The location of the brine/α-quartz interface is z = 0 as defined in Sec.II.B.1.Then we can get the two following integrations: where R c is the curvature radius of oil/brine interface.θ and R c can be simultaneously obtained from known R 0 and z m by combination of Eq.( 14) and Eq.( 15).

B. MD results: Fluctuations
Using a time-averaged θ in finite time to represent the real wettability of the MD system in equilibrium, we need to evaluate its fluctuations by analyzing standard error of the mean(SEM). 28y averaging θ every time period τ for n τ times, a sample of contact angle θ τ ,i (i = 1, 2, . . ., n τ ) is obtained.We denote θ averaged on the whole sample time ∆t = τn τ as θ ∆t , and the expectation of θ ∆t as θ * .Then standard error of θ ∆t is calculated by For a constant ∆t, it can be proved by statistics that the expectation of s ( θ ∆t ) is independent on τ, under the condition that τ is larger than correlation time of θ to make θ τ ,i (i = 1, 2, . . ., n τ ) independent.Assuming θ respects normal distribution, a 90% confidence interval of θ * can be estimated as θ ∆t − 1.645s ( θ ∆t ) , θ ∆t + 1.645s ( θ ∆t ) .
With the above basis, we can analyze the fluctuations of θ ∆t obtained from MD results, taking Case 3 and Case 7 as two examples.We first estimate the correlation time of θ as 2 ns, since s( θ ∆t ) converges when τ > 2 ns, as shown in Figure 5(a).In the following, we will always calculate s( θ ∆t ) using τ = 2 ns.The error development with time is shown in Figure 5(b) and (c).It is clear that near t = 50 ns, s( θ ∆t ) is already near 1 degree and θ ∆t almost does not change with elongated t.The 90% confidence intervals of θ * for Case 3 and Case 7 are respectively [49.7, 52.9] and [43.5, 46.8].
To simplify the symbols, in the following we will replace θ * by θ.

C. MD results: Representativeness of wettability
Due to the assumptions of homogeneous density and constant curvature, the defined contact angle may not well represent the drop shape, or wettability of the system.Here we evaluate this problem using the density contour plot obtained from MD simulation, as shown in Figure 6.The black lines FIG. 6. MD results of the density distribution (unit: g cm 3 ) for the liquid in Case 2. The black lines are the oil/brine interface calculated from θ and R c and the liquid/α-quartz interface defined by z = 0. represent the oil/brine interface calculated from θ and R c , and the liquid/α-quartz interface defined by z = 0. Clearly the black arc well coincides with the interface indicated by density distribution.Therefore, we ensure that the measured contact angle can well represent the wettability of the simulated system.

V. WETTABILITY VALIDATION
In this section, we will validate the wettability of the model system by comparing with experiments and EDL repulsion theory.

A. Comparison with experiments
Firstly, in order to check the representativeness of this simulated system for the real oil/brine/rock system, we compare the contact angle of the MD results with available experimental data.We carefully select the experiments whose materials and conditions are similar to our MD simulation.To compare Case 1-4, we select experiments with a low-pH brine for negligible surface charge; to compare Case 5-7, we select experiments with no charged polar groups in oil and a surface charge density of solid close to 0.067 C m 2 .As shown in Figure 7, the contact angle of MD simulation and the selected experiments are in the same level for an uncharged system, and in the same trend with concentration for charged solid surfaces.Actually the experimentally measured contact angle is largely dependent on measurement procedures.Readers may refer to Ref. 5,14,48,49 for more information.It is confident to say that the wettability of our modeled system is close to that of the real oil/brine/rock system.

B. Comparison with EDL repulsion theory
Next, we will briefly introduce the continuum theory of the contribution of EDL repulsion to ion-tuned wettability, and compare it with the MD results when negligible MIE effects exist.
The EDL repulsion between oil and rock is the only mechanism of ion-tuned wettability with a whole theory to quantitatively link brine salinity and pH with contact angle.Theory of EDL repulsion was first developed at the beginning of colloid science, applied in the famous DLVO theory to predict colloid cohesion. 46,50,51Then researchers specialized in oil industry combine a simplified hydration force model and wetting theory with DLVO theory, to quantitatively predict wettability alteration of OBR system by brine. 16,17,52Though this paper is only concerned about OBR systems, the theory can be extended to other fluid/fluid/solid systems, such as brine/vapor/quartz system. 53he following theoretical framework is based on previous classical theory, 16,17,52 with some assumptions clarified below.EDL repulsion emerges when EDLs overlap in a confined space.Assume constant σ 0,qz , the EDL repulsion pressure generated between plate oil/brine interface and plate brine/α-quartz interface can be calculated by the contact value theorem: 46 where c s i is the concentration of i-th ion on α-quartz surface, which can be easily calculated assuming both Stern layer and diffuse layer obey PB equation.
Assume h eq as the equilibrium water film thickness between oil and rock, and then the EDL repulsion energy between oil/brine interface and brine/rock interface can be calculated by Note that this equation can be seen as an integration of disjoining pressure in the transition zone near contact line, where the curvature is strongly influenced by surface forces, as indicated by the augmented Young-Laplace equation. 16,54To predict contact angle, the artificial assumption of h eq down to angstroms 52 can be an intrinsic defect of the theory.Assume the water-phase contact angle is θ 0 when w EDL = 0, and apply the Young-Dupre's equation, 16,54 we get  5) is estimated based on triple layer model of silica, 44,45 according to the concentration and the pH.
where γ ob is the surface tension of the oil/brine interface.γ ob can be calculated by the standard MD simulations, as described in the supplementary material.For the system shown in Figure 1, the γ ob is only weakly related with concentration, with values 55.1, 54.9, 56.5 mN m 1 for 0.05, 0.2, 1 M KCl solution.To predict θ of Case 5-7, the MD results of Case 1-3 are used to give the θ 0 .With above equations, the contribution of EDL repulsion to ion-tuned wettability can be calculated.Note that EDL repulsion can also be produced between one charged interface and one uncharged interfaces (Case 5-7), as indicated by this theory.Before comparing the predicted contact angle with MD results, we first prove the formation of EDL in the MD simulation by the ionic distribution.Figure 8 shows the distribution of K + and Cl -for case 6.It is clear that K + ions are attracted to the charged brine/α-quartz interface, while Cl -ions are repelled, and thus the existence of EDL is proved.In the supplementary material, we also present a quantitative comparison between the PB equation and MD results for a simple brine/α-quartz system without the oil drop.Next, we can compare the predicted contact angle with MD results, as shown in Figure 7 causes the contact angle decrease, and this decrease is more significant for lower concentration at the same surface charge density.Besides, the comparison is also quantitatively consistent with the assumed h eq in spite of the many assumptions in the theory.

VI. CONCLUSIONS
In this work, we build a reliable MD framework to study the ion-tuned wettability of oil/brine/rock systems.To start, we carefully choose the force field to build the MD system.Then we emphasize that how to assign ion numbers to balance interface charge can largely influence the bulk ion concentration of the system, which is important when studying the concentration effect on wettability.From both theoretical analysis and validation from MD results, we suggest to add counter ions to equilibrium interface charge.With the built system, MD simulations are run, and the contact angle is calculated by converting center mass coordinates.The fluctuation analysis indicates that the correlation time of the contact angle is around 2 ns, and the error can be controlled well within 40 ns sample time.In addition, the measured contact angle is well consistent with the interface shape indicated by the density contour plot.Finally, we validate the wettability of the MD system in two aspects.We prove that the wettability of our model system is close to that of the real oil/brine/rock system; and that the wettability alteration is qualitatively consistent with EDL repulsion theory when the quartz is charged but the oil is uncharged, in which case ion-binding does not exist.To conclude, this MD framework can be reasonably used to reveal the mechanism of ion-tuned wettability in the future.Our

FIG. 1 .
FIG. 1.The model system for MD simulation.Besides the geometry parameters defined here, the y-direction thickness is denoted as W.
FIG. 8. MD results of the ionic distribution (unit: M) for (a) K + , and (b) Cl -in Case 6.The white lines are the oil/brine interface calculated from θ and R c and the liquid/α-quartz interface defined by z = 0.

TABLE I .
Properties of water, decane, and water/decane interface calculated by MD simulation with SPC/Fw water model and opls-aa decane model, at T = 300K and p = 1 atm.Experiment data for water properties comes from Ref. 33, for decane comes from Ref. 39, and for water/decane interface tension (T = 300.65K) comes from Ref. 40.

TABLE II .
A summary of all the simulated cases for the model system.N K and N Cl are respectively the total number of K + and Cl in the simulated system.Case 4 is kept at T = 373.15K, while other cases are all kept at T = 300 K.