Nanoconfined Electrolyte Solutions in Porous Hydrophilic Silica

Article pubs.acs.org/JPCC

Nanoconfined Electrolyte Solutions in Porous Hydrophilic Silica Membranes R. Renou,†,⊥ A. Ghoufi,*,‡ A. Szymczyk,*,† H. Zhu,§ J.-C. Neyt,∥ and P. Malfreyt∥ †

Institut des Sciences Chimiques de Rennes, UMR 6226 CNRS, Université de Rennes 1, 263 avenue du Général Leclerc, 35042 Rennes, France and Université Européenne de Bretagne, France ‡ Institut de Physique de Rennes, UMR 6251 CNRS, Université de Rennes 1, 263 avenue du Général Leclerc, 35042 Rennes, France § College of Environmental Science and Engineering, Tongji University, Mingjing Building, 1239 Siping Road, Shanghai, P.R. China ∥ Institut de Chimie de Clermont-Ferrand, ICCF, UMR CNRS 6296, BP 10448, F-63000 Clermont-Ferrand S Supporting Information *

ABSTRACT: Manipulating of electrolyte solutions is ubiquitous in numerous fields ranging from condensed to soft matter. The fundamental understanding of the molecular mechanism ruling the macroscopic properties is inescapable to well rationalize the engineering processes. In confined situations such as nanofluidic or nanofiltration, interactions between the fluids and the surface drastically modify the usual macroscopic properties. In this work we addressed a structural and dielectric analysis of confined NaCl, NaI, MgCl2, and Na2SO4 in nanoporous hydrophilic silica membranes. We highlight a dielectric anisotropy of confined water and an unusual increase in the radial component of permittivity (ε⊥) of confined solutions with respect to the decrease in ε in the bulk solution phase. Structurally, we investigate the size, the polarizability, and the salt effect on water confinement and ions translocation through a cylindrical silica nanopore.

■

INTRODUCTION Understanding of the behavior of water and electrolyte solutions in confined environments is of the utmost importance for a variety of fields ranging from biophysics1,2 to material science.3 Indeed, the development of nanoporous materials requires a detailed molecular-level understanding of the behavior of both solvent and electrolyte at interfaces and under confinement.4 A number of experimental and theoretical studies have provided molecular-level insights into the behavior of interfacial water, and it has been reported that both structural and dynamical properties of interfacial water are significantly affected by the solid substrate features, resulting in a different behavior compared to that of bulk water.5 However, many fundamental properties of confined water have yet to be fully understood.6 Confined electrolyte solutions are also of interest in many practical applications in the field of chemical engineering, including membrane processes.7,8 For example, nanofiltration (NF) has attracted increasing attention over the recent years and has proved to be an energy efficient and environmentally friendly process for desalination and water treatment.9−12 Thanks to the significant improvements performed in research and development of membrane materials over the past decade, NF has today the power to solve many separation problems in an economically viable way. Nevertheless, its current industrial applications are small in number compared to the potential applications that still await. A major reason is that the physical © XXXX American Chemical Society

phenomena involved in the separation by nanoporous membranes are still poorly understood, and for all current descriptions of NF, the approach is more descriptive than truly predictive due to the use of a variety of empirical fitting parameters in theoretical models used to describe transport.13 Indeed, fluids confined in dimensions that are of no more than an order of magnitude larger than molecular sizes exhibit properties that differ considerably from those of bulk fluids,14 and some peculiar behaviors, like for example the layering structure of solvent molecules near a solid surface,15 are not accounted for in the transport theories. Molecular simulations can explicitly account for the molecular nature of both ions and water, and they can also treat surfaces with atomically detailed acurracy.4 As such, they have become widespread in the literature and proved to be a powerful method to investigate confinement16−23 and interfacial effects.24,25 Thus, Ho et al. have studied the structural aspects of aqueous NaCl and CsCl solutions confined in crystalline slit-shaped silica nanopores of varying degrees of protonation,26 while Beu focused on transport in aqueous ionic solutions through carbon nanotubes.27,28 Almost at the same time Videla et al. have explored the confinement effect within the functionalized surface on structural and dynamic properties Received: December 9, 2012 Revised: April 26, 2013

A

dx.doi.org/10.1021/jp403450x | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

of aqueous electrolytes solutions.29 More recently, Bonnaud et al. have shown the impact of the charged surface of a slit silica pore on the diffusion coefficient of water and ions.30 More recently confinement of organic solvents and their aqueous mixtures has been surveyed. For example, Melnikov et al. investigated the effect of the nanoscale confinement on the properties of a binary acetonitrile−water mixture.31 In this work, we report the structural and dielectric behavior of water and aqueous electrolyte solutions under nanoconfinement. We employed all-atom molecular dynamics simulations at equilibrium to investigate water and various aqueous electrolyte solutions confined within hydrophilic silica nanopores which are the focus of significant scientific interest due to their potential in a wide variety of applications including desalination membranes.32 The effect of the electronic polarizability of ions on the dielectric properties of confined salt solutions was also assessed by considering either nonpolarizable force fields or polarizable force fields for ions.

massless shell, connected by a harmonic spring. The core and shell carry different electric charges, the sum of which equals the charge on the original atom. Details about the core−shell model can be found in ref 38. The induced dipole model42 and fluctuation charge model43 are the explicit models involving a large increase in computational time. Parameters are given in Table 2 and Table 3. We checked that the polarizability effect Table 2. Lennard-Jones Parameters for Nonpolarizable Ionsa

MODEL AND COMPUTATIONAL PROCEDURE Model. In this study, an atomistic hydrophilic silica nanopore has been investigated. The silica framework has been built by carving a cylinder of diameter D = 12 Å in cubic amorphous silica of 71 Å along the z axis.33 Saturation of oxygen atoms has been accounted for by considering the procedure developed by Bródka and Zerda34 to get an irregular inner hydrophilic surface of the porous silicate (7.5 SiOH nm−2 which corresponds to the highly hydrated protonated silica pore34,35). Although the silica matrix was subsequently kept rigid, rotation around the Si−O bond of the hydroxyl groups was allowed from the application of the SHAKE constraints algorithm on O−H bonds of the silanol groups.17,36,37 Though the SHAKE algorithm works well for time steps up to 3 fs, we checked the independence in Δt by using Δt = 1 fs, Δt = 2 fs, and a flexible model (stretching [O−H] and bending potential [Si−O−H]) with Δt = 2 fs. These results are given in Table S3 of the Supporting Information (SI). Table S3 (SI) shows that the time step does not impact the total energy, volume, and overall dielectric permittivity obtained with the SHAKE algorithm. Additionally, we show that the constraint O−H is equivalent to the flexible Si−O−H model in terms of overall dielectric permittivity. Charges and Lennard-Jones parameters of the different sites are given in Table 1. Ions were modeled by

Si Onb Ob Hb a

ε

u.e.

Å

kJ mol−1

1.27265 −0.6350 −0.5340 0.2060

0.000000 1.912131 1.912131 0.000000

0.00 2.70 3.00 0.00

2.21 4.42 4.812 5.40 4.42 1.85 3.55 3.15

a

S and Os are the labels of the sulfur and the oxygen atoms in the sulfonate group, respectively. I−m is the modified iodide where σ has been modified.

Table 3. Lennard-Jones and Charges Parameters for the Polarizable Core−Shell Model of Ions38 +

Na Cl− I−

qcore

qshell

ε (kJ mol−1)

σ (Å)

1.69 3.46 4.73

−0.69 −3.46 −4.73

0.132 0.301 0.871

2.92 4.96 5.52

was well allowed for through the modeling of the electrolyte solutions at the air−water interface. In Figure 1 we established that the core−shell model allows us to recover the polarizability effect close to the interface, i.e., a propensity of Cl− with respect to Na+. Interestingly, we show in Figure 1a that the nonpolarizable model of ions is inefficient to reproduce this interfacial propensity. For the unpolarizable force field we used the recent and robust force field developed by Reif and coworkers.44 Intermolecular parameters of the silica nanopore and ions are provided in Tables 1−3. Labels are described in Figure 3. As the polarizability of water molecules does not impact the ion density in the interfacial region,45 we considered here the nonpolarizable TIP4P/2005 water model.46 This one has been found to provide an impressive performance for a variety of physical properties.47 In Figure 1 and Figure 2 we check that the polarizability of water molecules slightly impacts the interfacial location of ions and the axial and radial density of water in confined medium. This result allowed us to opt for an unpolarizable water model. With the core−shell description we found a small Na−Cl separation distance of 2.1 Å, while 4 Å was found by Jungwirth and Tobias.48 This difference is due to the force field parameters of the drude oscillator’s representation. It is possible to modify the Lennard-Jones parameters to recover this distance, but this is not the scope of this study. The main physics ingredient is the ability to reproduce the interfacial ions propensity. Computational Procedure. Ion transport trough the nanopore was investigated from the molecular dynamics simulations into isothermal−isosurface−isobaric statistical ensemble (III).17 With this approach an explicit planar solid (S)/liquid (L) interface was considered such that the volume of

Table 1. Force Field Parameters of the Silica Matrix17,a σ

σ (Å)

0.407 0.493 0.150 0.409 0.409 3.662 1.046 0.837

Na Cl− Br− I− I−m Mg2+ S(SO42−) O(SO42−)

■

q

ε (kJ mol−1) +

Labels (Onb, Ob, and Hb) are define in Figure 3.

considering a polarizable model38,39 and an unpolarizable model.40 Ion polarization was taken into account by means of a core−shell model.38,41 The reader is redirected to ref 41 for a full description of the polarizable core−shell model. With this model, the polarization has been accounted for from an oscillator mimicking the deformation of the electron cloud without significantly increasing the computational effort. Thus, the polarizable atom is represented by a massive core and B



Article

Figure 3. Description of the protonated silica nanopore. Ob is the oxygen bounded to a hydrogen atom (Hb), while Onb is unbounded.

Figure 1. a) and b) Air−water density profiles of water, Cl−, and Na+ for electrolyte solutions in NaCl at 2 M and 298 K from both polarizable and nonpolarizable models. In a) and b) right and left axes correspond to the water and ion densities, respectively. Water was modeled from the TIP4P/2005 model.46 (a) Corresponds to an enlargement of the interfacial region of b).

Figure 4. Illustration of the initial configuration. A silica nanopore surrounded with two reservoirs of electrolyte solution. Red indicates the oxygen atoms, white the hydrogen, and yellow the silicon atoms. Blue and green colors correspond to the Na+ and Cl− ions, respectively.

the silica material was kept fixed while the fluctuations of the liquid density were controlled by an anisotropic barostat along the z direction. With this method, a finite pore with an inlet and an outlet was modeled. As shown in Figure 4, external S/L interfaces were built by inserting two reservoirs of electrolytic solutions surrounding the silica nanopore. The III statistical ensemble was carried out using the Berendsen barostat49 with isotropic periodic boundary conditions at 298 K. Times of 0.1 and 0.5 ps were used for the thermostat and barostat relaxation time, respectively. The Berendsen algorithm allows us to get the small oscillations of the volume in contrast to the Hoover algorithm. Further, the Berendsen algorithm is well-known to conserve total momentum. Though the thermostat does not

generate a correct canonical ensemble, for large systems, the approximation yields roughly correct results. In the Supporting Information we compared the Hoover and Berendsen algorithm. We show that both methods provide similar energy and volume. Additionally, as shown in Figure S1c of the SI the radial water distributions calculated from Hoover and Berendsen simulations are identical. This allows us to conclude that our results are not dependent on the integration’s algorithm. The initial simulation box was a rectangular parallelepipedic box of dimensions LxLyLz (Lx = Ly = 35.5 Å) and Lz = 210 Å. Systems with explicit interfaces consisted of 5555 water molecules and 100 ion pairs corresponding to a concentration of c = 1 mol L−1 in the reservoirs. The LennardJones interactions were truncated at 12 Å, and the electrostatic

Figure 2. Axial (a) and radial (b) water density confined into the silica with a pore radius of 6 Å. C



Article

interactions were computed from Ewald summation.36 The equations of motion were solved using the velocity Verlet algorithm36 with a time step of 2 fs. Averages were performed for the last 10 ns of simulations after 5 ns of equilibration. All simulations have been carried out from the modified DL-POLY package.50

■

LOCAL DIELECTRIC PERMITTIVITY CALCULATION In linear, homogeneous, isotropic media ε is a scalar. However, in linear anisotropic media it is a tensor of rank 2, and in nonhomogeneous media it is a function of position inside the medium. Therefore, it is difficult to consider an overall dielectric constant under cylindrical confinement.18,19,51,52 Thus ε is defined as a symmetric tensor involving a diagonal matrix.18,19 The cylindrical confinement presents a translational invariance according to z and ϕ (ϕ is the angle in cylindrical coordinates), so all averaged fields and observables only depend on r. The dielectric tensor is diagonal with only two unique components, radial (ε⊥(r)) and axial (ε||(r)) to the curved surface. As shown in ref 18 the local operational expression to compute the local permittivity in Cartesian components is c εαα (r ) = 1 + βε0−1[⟨mα (r )Mα⟩0 − ⟨mα (r )⟩0 ⟨Mα⟩0 ]

(1) Figure 5. Radial (a) and axial (b) densities of center of mass of water molecules in both matrices R6 and R12. (b) The axial water profile close to the interface between z = 35 Å and z = 50 Å. The horizontal straight line corresponds to the water density in the bulk phase, while the vertical dashed line corresponds to the location of the planar interface.

where α = {x,y,z}; r = (x + y ) ; m(r) is the polarization density such as M = ∫ Vm(r)dr; M is the total dipolar moment; ⟨...⟩ denotes a statistical average over the different configurations; β is the inverse temperature; and ε0 is the vacuum permittivity. Local permittivity components were computed from the dipolar moments of water molecules and the dipolar contribution of ion pairs. From the computation of radial distribution function (RDF) between Na+ and Cl− we have determined a separated distance of 3.5 Å. RDFs between Na+ and Cl− in bulk and confined phases are given in Figure S3 of the Supporting Information. While the intensities are different, the location of the first peak is similar. Thus, the distance used to calculate the average dipole moment of ions corresponds to the distance of the first shell of coordination (dotted line in Figure S3 of SI). The contribution of the dipole moment of individual ions stemming from the core−shell relaxation was found negligible as shown in Figure S4 of the SI. Figure S4 clearly highlights a very slight influence of the dipolar relaxation of core−shell. M was then computed by considering the individual dipolar moment of water molecules and ion pairs M(t) = ΣNi waterμi(t) + ΣNj pairsμj(t) where M is the total dipole moment of the system given with respect the individual molecular dipole moment of the ith water molecule and jth ions pair. The radial and axial components correspond to the xx (or yy since in cylindrical symmetry the xx and yy components provide the same contribution 18,19) and zz direction, respectively. Quadrupolar moments were not included in our calculations since it has been shown that the contribution is negligible under cylindrical confinement.18 2

2 1/2

We observe an increase in density close to the cylindrical interface with respect to the bulk water and a layering structure. This is imputed to a confinement effect and a pronounced excluded volume near the interface. For R12, by moving away from the interface the density decreases, and layering organization vanishes because the water molecules no longer feel the solid interface. In both matrices the period of oscillation is about 3−4 Å, which basically corresponds to the diameter of the water molecule. In Figure S5 of the SI we show that whatever the chemical nature of the confined molecular fluid the oscillation period corresponds to the molecular diameter. Additionally, Figure S5 of the SI shows that the chemical nature of the interface does not impact the layering organization of the adsorbed fluid. This clearly suggests that a confinement effect is the main ingredient of the layering structure. For R6 the density of bulk water is not recovered at the center of the pore (see Figure 5a). This results from the drop in pore size. Axially, this layering organization is also observed near the planar interface. As shown in Figure 5b, two layers of 3−4 Å are observed. Dweik et al. have already evidenced this planar anchoring at the external planar surface.53 To study the influence of the interfacial planar silanol rate we constructed a matrix of R = 12 Å with 6.6 SiOH/nm2 (R122) at the planar interface. As shown in Figure 6 we observe a shift in the location of the first water layer for R122 with respect to R12. Indeed, the first water layer of R12 is located at 37 Å, i.e., close to the interface located at 35.5 Å. In R122 the first layer is located at 40 Å with R122. From Figure 7a and Figure 7b a homogeneous density in SiOH is observed with R12, while from R122 two peaks in SiOH are observed beyond the interface. As the water molecules hydrate the SiOH groups from hydrogen bonds the water density is then correlated to the SiOH density. Thus the depletion in water density close to

■

RESULTS AND DISCUSSION Confined Water. First of all, we studied the radial structure and the dielectric response of the confined water in two pores of radius R = 6 Å (R6) and R = 12 Å (R12). Let us note that both pores have exactly the same surface chemical composition, i.e., 7.5 SiOH/nm2, corresponding to highly hydrophilic pores.35,37 We report in Figure 5a the radial (ri = (x2i + y2j )1/2) density of the center of mass of water molecules for both R6 and R12 with xi and yi being the coordinates of particle i. D



Article

the planar interface at 37 Å in Figure 7 is connected to the decrease in density of SiOH. This modification in silanol rate in the external planar surface involves a change in behavior of the density water. This effect has been pointed out by Laria et al. with hydrophobic and hydrophilic pores of R = 15 and R = 10 Å.29 We report in Figure 8a the pair correlation functions (g(r)) between the oxygen atoms of water for both the matrices. Calculation of g(r) was not corrected of the excluded volume. The signature of the geometrical confinement, i.e., the existence of spatial regions where particular atomic species are not allowed, known as the excluded volume effect, is that RDFs inside the pore do not oscillate around 1, as expected for a bulk system. This correction is well detailed in ref 54. To highlight the confinement effect, we do not correct the radial distribution functions. RDFs between oxygen atoms of water with corrections are provided in the Supporting Information. Figure 8a depicts that the confinement and size radius affect the arrangement of water molecules into the pore. Indeed, the first and the second hydration shells are located at the same position, while their integration provides different coordination number, i.e., 3.9 for R12, 3.5 for R6, and 4.0 in the bulk phase. Concerning the hydrogen bonds of water (HB) we report in Figure 8b their radial density. Hydrogen bonds were calculated using a quantum criterion17−19,37 such as the distance between the oxygen and hydrogen atoms of two water molecules is

Nanoconfined Electrolyte Solutions in Porous Hydrophilic Silica

Recommend Documents