The structure of a concentrated solution of NaCl in D2O was investigated by in situ high-pressure neutron diffraction with chlorine isotope substitution to give site-specific information on the coordination environment of the chloride ion. A broad range of densities was explored by first increasing the temperature from 323 to 423 K at 0.1 kbar and then increasing the pressure from 0.1 to 33.8 kbar at 423 K, thus mapping a cyclic variation in the static dielectric constant of the pure solvent. The experimental work was complemented by molecular dynamics simulations using the TIP4P/2005 model for water, which were validated against the measured equation of state and diffraction results. Pressure-induced anion ordering is observed, which is accompanied by a dramatic increase in the Cl–O and O–O coordination numbers. With the aid of bond-distance resolved bond-angle maps, it is found that the increased coordination numbers do not originate from a sizable alteration to the number of either Cl⋯D–O or O⋯D–O hydrogen bonds but from the appearance of non-hydrogen-bonded configurations. Increased pressure leads to a marked decrease in the self-diffusion coefficients but has only a moderate effect on the ion–water residence times. Contact ion pairs are observed under all conditions, mostly in the form of charge-neutral NaCl0 units, and coexist with solvent-separated Na+–Na+ and Cl–Cl ion pairs. The exchange of water molecules with Na+ adopts a concerted mechanism under ambient conditions but becomes non-concerted as the state conditions are changed. Our findings are important for understanding the role of extreme conditions in geochemical processes.

Sodium chloride in water forms the basis of a wide variety of the Earth’s geofluids, such as brine, for which the state conditions can vary over a broad range of temperatures and pressures.1 On the ocean floor at continental margins, for instance, methane-containing clathrate-hydrate structures can form at a temperature T below ∼288 K and pressure p above ∼120 bars.2 In comparison, on the ocean floor at hydrothermal vents, supercritical conditions can exist where the temperature and pressure exceed values corresponding to the critical point of seawater at 680 K and 298 bars.3 At subduction zones, temperatures of ∼523 K at 10 kbar and ∼773 K at 20 kbar are typical along subduction pathways,4 and the salinity of aqueous fluids can reach the equivalent of 20 wt. % NaCl,5 corresponding to an aqueous solution of 4.28m NaCl (1m = 1 molal). The concentration of halite-bearing fluid inclusions in hydrothermal ore-forming systems can be more than 50 wt. % NaCl, where these inclusions were trapped within growing crystals and therefore provide information on the nature of the fluid phase from which the crystals were formed.6–8 Further afield, a subsurface NaCl-rich ocean may be a prominent geological feature on Jupiter’s moon Europa.9 

It is a challenge to create realistic models for the geochemical processes that take place under extreme conditions, given their scope and the paucity of experimental information on the structure of geofluids. These fluids include aqueous solutions of NaCl and the adaptations that originate from, e.g., the dissolution of CO2 and/or silicates.1,10 Here, it is desirable to have a molecular-scale understanding of fundamental phenomena such as ion speciation, e.g., the fractions of fully dissociated single ions vs diatomic ion pairs (NaCl0) and triple clusters (Na2Cl+ or NaCl2).11 The nature of this speciation is sensitive to temperature and pressure largely because of changes to the dielectric constant of water, which decreases with increasing temperature at fixed pressure but increases with pressure at fixed temperature.12,13

Knowledge of the ion speciation is necessary to constrain the parameters used in models for predicting the standard-state thermodynamic properties of ions and complexes. Here, a widely used model is the Helgeson–Kirkham–Flowers equation of state,14 which delivers results that depend on an accurate description of the dielectric constant of water and its variation with temperature and pressure.15 Until recently, this dependence restricted its use to pressures of up to 5 kbar.13 In addition, the model uses the dielectric continuum theory of Born,16 which assumes that the dielectric constant of the solvent near to an ion is the same as that in bulk solution, an over-simplification that can affect the calculation of equilibrium constants under extreme conditions.17 These constants express the relative concentrations of ion complexes and are necessary to understand phenomena that include the dissolution of rock, the transport and deposition of ore minerals, and the treatment of corrosion.10,12,15 From a topological point of view, it is important to understand how the hydrogen-bonded network of the aqueous solution reorganizes in order to accommodate the pressure-induced packing of the chemical species.18 

The structure of aqueous NaCl solutions under extreme conditions is difficult to access by experiment because of the corrosive nature of the samples and the need to contain them within apparatus that is suitable for the chosen structural probe. Here, we address this challenge by applying the method of neutron diffraction with chlorine isotope substitution to measure the structure of an aqueous solution of 5m NaCl in D2O with (i) temperature increasing from 323 to 423 K at 0.1 kbar and (ii) pressure increasing from 0.1 to 33.8 kbar at 423 K. Under these conditions, the static dielectric constant of water ɛ is cycled, decreasing to 56% of its ambient conditions value with increasing temperature before returning to this value with increasing pressure (see Sec. IV B 4). The structure of water also varies dramatically, with the nearest-neighbor O–O coordination number increasing from 4.419 to ≃11–12.20–22 The neutron diffraction results are interpreted with the aid of molecular dynamics simulations using the TIP4P/2005 model for the water molecules.23 The Na+ and Cl ions were given full formal charges of +e and −e, respectively, where e is the elementary charge, and the Lennard-Jones parameters were taken from the work of Yagasaki et al.24 

The neutron diffraction method was chosen for the experiments because it provides site-specific information on the coordination environment of the Cl ion,25 and neutrons are readily transmitted through high-pressure apparatus. In contrast, there are challenges in the application of x-ray absorption spectroscopy to solutions containing ions of light elements because of the small energy of the K-edge, which limits the transmission of x rays through the solution and its container.26,27 The structure of aqueous NaCl solutions has recently been investigated by x-ray diffraction at room temperature and pressures in the gigapascal regime (1 GPa = 10 kbar),28 although site-specific structural information was not accessible and the measured diffraction patterns were restricted to a limited range of scattering vectors.

This paper is organized as follows. The essential neutron diffraction theory is given in Sec. II, and the molecular dynamics simulations are outlined in Sec. III. The neutron diffraction experiments are described in Sec. IV. In the data analysis, the density was taken from experiment or, at the highest pressures, from the molecular dynamics simulations. In Sec. V, the experimental results are compared to those obtained from the simulations in order to test the validity of the force field used in the molecular dynamics work. The simulations are then used to give a more comprehensive account of the effect of temperature and pressure on the structure and dynamics. The results are discussed in Sec. VI, where issues include the nature of the hydrogen-bonded network, the effect of pressure on the ordering of the anionic (Cl and O) species, and the Na+–water exchange mechanism. Conclusions are drawn in Sec. VII.

In a neutron diffraction experiment on a solution of NaCl in D2O for which the incident neutron wavelength λ is fixed, the diffracted intensity is measured as a function of the scattering angle 2θ.29 Let Isc(θ) and Ic(θ) represent the background-corrected intensities measured for the sample (s) in its container (c) and for the empty container, respectively. The differential scattering cross section for the sample can be written as30 

dσdΩ=1NsAs,sc(θ)Isc(θ)a(θ)Msc(θ)Ac,sc(θ)Ac,c(θ)Ic(θ)a(θ)Mc(θ),
(1)

where Ns is the number of illuminated sample atoms; a(θ) is a normalization factor measured by reference to a vanadium standard;31,Ai,j(θ) denotes the attenuation factor for neutrons that are scattered in medium i and attenuated, through absorption and scattering, in medium(s) j; and Mi(θ) is the cross section for multiple scattering in medium(s) i. For cylindrical samples, Ai,j(θ) can be calculated by using the method of Paalman and Pings32 and Mi(θ) can be evaluated within the quasi-isotropic approximation by using the method of Soper and Egelstaff.33 The use of Eq. (1) ensures that the attenuation corrections are applied to once-scattered neutrons.

The differential scattering cross section of Eq. (1) can also be written as

dσ/dΩ=F(Q)+P(Q),
(2)

where θ is related to the magnitude of the scattering vector via Q = (4π/λ)sin θ. The total structure factor

F(Q)=αβcαcβbαbβSαβ(Q)1
(3)

describes the structure of the system, where cα and bα are the atomic fraction and bound coherent neutron scattering length of chemical species α and Sαβ(Q) is a Faber–Ziman partial structure factor.34 The term P(Q) has contributions from (i) the self-scattering from individual nuclei and (ii) the inelastic scattering associated with both the self and distinct parts of dσ/dΩ, for which there is no exact theory.29,35,36

The total pair-distribution function is obtained from the Fourier transform

G(r)=12π2ρr0dQQF(Q)M(Q)sin(Qr),
(4)

where r is a distance in real space and ρ is the atomic number density of the sample. The modification function M(Q) = 1 for QQmax and M(Q) = 0 for Q > Qmax originates from the maximum value Qmax that can be accessed by using a diffractometer. If Qmax is sufficiently large that M(Q) does not truncate any of the large-Q oscillations in F(Q), it follows that

G(r)=αβcαcβbαbβgαβ(r)1,
(5)

where gαβ (r) is a partial pair-distribution function. The low-r limit of Eq. (5) is given by G(r → 0) = −∑αβcαcβbαbβ. The mean coordination number of atoms of type β contained within a volume defined by two concentric spheres of radii rmin and rmax centered on an atom of type α is given by

n̄αβ=4πρcβrminrmaxdrr2gαβ(r).
(6)

The running coordination number is found by setting rmin = 0 and by finding the change in n̄αβ as rmax is increased from zero.

Let neutron diffraction experiments be performed on two solutions that are identical in every respect, except that chlorine with isotopic abundance Cl is replaced by chlorine with isotopic abundance Cl, where b†Cl > b‡Cl. The samples differ only in the isotopic enrichment of chlorine, so their scattering and absorption cross sections are very similar. In consequence, there is little change to either the magnitude or θ-dependence of As,sc(θ), Ac,sc(θ), and Msc(θ). If the measured samples in container intensities are denoted by Isc(θ) and Isc(θ), it follows from Eq. (1) that

ΔFCl(Q)dσ/dΩ|Cldσ/dΩ|Cl1NsAs,sc(θ)Isc(θ)Isc(θ)a(θ),
(7)

where As,sc(θ) can be calculated for either of the samples. The difference function method therefore leads to a cancellation of systematic errors that arise from the container, background, and multiple scattering corrections, a feature that will be exploited in the following work on aqueous solutions under extreme conditions.

From Eq. (2), it follows that Eq. (7) can be re-written as

ΔFCl(Q)=ΔFCl(Q)+ΔPCl(Q),
(8)

where the first-difference function is given by

ΔFCl(Q)=FCl(Q)FCl(Q)=ASClD(Q)1+BSClO(Q)1+CSClNa(Q)1+DSClCl(Q)1,
(9)

with A = 2cClcDbDΔbCl, B = 2cClcObOΔbCl, C = 2cClcNabNaΔbCl, D = cCl2bCl2bCl2, and ΔbCl = b†Clb‡Cl. Hence, the main contribution from inelastic scattering is eliminated on forming Eq. (8) to leave a small residual correction ΔPCl(Q)cClbCl2̄bCl2̄, where bα2̄=bα2+binc,α2 and binc,α is the incoherent scattering length of chemical species α.37 

As shown in Eq. (9), the first-difference function simplifies the complexity of correlations associated with F(Q) and provides site-specific information on the coordination environment of the chloride ion. The corresponding real-space information is obtained from the Fourier transform

ΔGCl(r)=12π2ρr0dQQΔFCl(Q)sin(Qr)=AgClD(r)1+BgClO(r)1+CgClNa(r)1+DgClCl(r)1.
(10)

The low-r limit is given by ΔGCl(r → 0) = −(A + B + C + D). In the implementation of Eq. (10), a Lorch modification function M(Q) = sin(πQ/Qmax)/(πQ/Qmax) for QQmax and M(Q) = 0 for Q > Qmax can be used to smooth Fourier transform artifacts.38,39

Molecular dynamics simulations were performed using the MetalWalls code.40 The 5m solution was represented by a system of 82 NaCl ion pairs and 818 D2O molecules. Pure water was represented by a system of 1000 D2O molecules. The rigid TIP4P/2005 model23 was used for the water molecules, in which the intra-molecular O–D bond length is fixed at 0.9572 Å and the intra-molecular D–O–D bond angle is fixed at 104.52°. This model was selected because it yields a phase diagram in good agreement with experiment.41 For the Na+ and Cl ions, the Lennard-Jones parameters were selected from the work of Yagasaki et al.24 because they were specifically designed to reproduce the solubility of the salt, and Lorentz–Berthelot mixing rules were used.

The systems were simulated at the following state points: (298 K, 0.001 kbar); (323 K, 0.001 kbar); (373 K, 0.001 kbar); (423 K, 0.001 kbar); (423 K, 1 kbar); (423 K, 2.4 kbar); (423 K, 10 kbar); (423 K, 22.4 kbar); (423 K, 27.9 kbar); and (423 K, 33.8 kbar). In all cases, a first NpT simulation of 200 ps was performed using an isotropic barostat (made of a Nosé–Hoover chain of length 5 with a time constant of 500 fs) coupled to a thermostat (also made of a Nosé–Hoover chain of length 5 with a time constant of 100 fs),42 where N denotes the number of particles. For each state point, the equilibrium volume V from the NpT simulation was selected and the system was simulated for 2 ns in the NVT ensemble using a thermostat43 made of a Nosé–Hoover chain of length 5 with a time constant of 1000 fs. The time step was 2 fs for all the simulations. The Coulomb interactions were calculated using an Ewald summation, and a 12 Å cutoff was used for the Lennard-Jones interactions (a tail correction was applied for the calculation of the pressure tensor). The rigidity of the water molecules was enforced using the rattle algorithm.44 

Note that in the NpT simulations, the fluctuation in pressure is typically ±0.5 kbar, so the calculated parameters at ambient pressure (0.001 kbar) and 0.1 kbar cannot be discriminated. Accordingly, the lowest pressure point will be labeled 0 kbar.

Each G(r) function was initially obtained by combining the calculated gαβ(r) functions using Eq. (5). The corresponding F(Q) function was then obtained by (i) Fourier transforming G(r) and by (ii) calculating the Fourier components that contribute toward the Sαβ(Q) functions in the low-Q region45 and combining them according to Eq. (3). The datasets were then joined smoothly at Q ∼ 2.5 Å−1 using the direct results at smaller Q and the Fourier transformed results at larger Q. In this way, it is possible to avoid the Fourier transform artifacts that originate from the finite size of the simulation cell. To provide a direct comparison with the diffraction results, the final F(Q) function was Fourier transformed to G(r) using Eq. (4) where the modification function M(Q) was chosen to match that used in the experiments.

The bond-angle distribution B(ϕ) for bond angle ϕ was calculated using a cutoff distance between atoms of type α and β set at the minimum after the first peak in gαβ(r) found for the solution at 423 K and 0 kbar. Here, the calculated values are proportional to the number of bonds between ϕ and ϕ + Δϕ, which is dependent on the solid angle ΔΩ ∝ sin(ϕ) subtended at that value of ϕ. Each bond-angle distribution is therefore plotted as B(ϕ)/sin(ϕ) to compensate for the effect of ΔΩ. In this way, a finite bond-angle distribution at ϕ ≃ 180° will not be artificially suppressed.46 

Gaseous H*Cl was prepared by the action of concentrated H2SO4 on dry salt in a flask connected to a vacuum line (∼10−5 Torr), where *Cl represents the isotope 35Cl (99% enrichment), chlorine of natural isotopic abundance natCl, a 1:1 mixture of 35Cl and 37Cl isotopes mixCl, or the isotope 37Cl (99% enrichment). The vacuum line was connected via a Young’s tap to glassware comprising two joined flasks. One flask was cooled with liquid nitrogen to solidify the gaseous H*Cl. The glassware was then sealed and removed from the vacuum line, distilled water was introduced into the second flask, and solidified H*Cl was warmed gently to ensure slow dissolution in the water, which was stirred vigorously. The H*Cl solution was neutralized using a 0.5 N standardized solution of NaOH in a titration setup, where the pH of the product was monitored continuously and the titration was terminated when the pH ≃ 6. Finally, the neutralized Na*Cl solution was dried thoroughly at 353 K, and the product was dissolved in D2O (99.9% enrichment, Sigma-Aldrich) to give a solution of 5m Na*Cl in D2O. The use of D2O avoids the large incoherent scattering cross section associated with light hydrogen and reduces inelastic scattering.

The neutron diffraction experiments were performed using two different high-pressure setups mounted on the instrument D4c at the Institut Laue-Langevin.47 This instrument was chosen on account of its (i) high count-rate across a wide Q-range, (ii) ability to access this Q-range with high-pressure setups,30 and (iii) high stability, which leads to reproducibility of the measured diffraction patterns.48 

1. Ti–Zr pressure cell experiment

The first setup employed a cylindrical Ti0.676Zr0.324 cell of inner diameter 5 mm and wall thickness 5.5 mm49 to access pressures of up to 1 kbar and temperatures of up to 423 K. The Ti–Zr alloy was chosen because it has a mean coherent neutron scattering length of zero, and its oxide layer offers good corrosion resistance to aqueous solutions under high pressure and temperature conditions.50 The cell was warmed using cartridge heaters set into clamps at the top and bottom of the cell. A separator of novel design was used to prevent mixing between the sample and the pressurizing fluid (Fluorinert FC-770).50 The cell, separator, and associated components were filled using a procedure that was designed to avoid trapping gas bubbles within the sample in the high-pressure system.51 A locking-pin placed in the base of the cell enabled it to be mounted with a reproducible orientation, which is important because the Ti–Zr alloy is not mixed ideally, so there are small Bragg peaks that show preferred orientation.52 The incident neutron wavelength was 0.4987(1) Å.

Diffraction patterns were measured for solutions of 5m NanatCl or Na37Cl in D2O held within the Ti–Zr cell either (i) at a pressure of 0.1 kbar and a temperature of 323, 373, or 423 K or (ii) at a pressure of 0.5 or 1.0 kbar and a temperature of 423 K; the empty Ti–Zr cell at a temperature of 323, 373, or 423 K; the empty instrument; and a cylindrical vanadium rod of diameter 6.37 mm for normalization purposes. The sample in cell and empty cell run times were optimized by using the procedure described in Ref. 53. The data analysis followed the procedure described in Sec. II.

2. Paris–Edinburgh press experiment

The second setup employed a VX5 Paris–Edinburgh (PE) press with single toroid anvils to access pressures of up to 33.8 kbar, where the anvils were warmed using cartridge heaters.54 Encapsulated Ti0.676Zr0.324 gaskets were used to contain the sample.55 Here, a pipette was used to load 30 µl of sample into each half of the gasket. The high surface tension of the solution enabled one of the gaskets to be inverted and placed on top of the other, thus ensuring that no gas bubbles were trapped. Excess sample was wiped away, and the gaskets and their outer Ti–Zr toroid were sealed in the anvils of the press using a small applied load. The incident neutron wavelength was 0.4971(1) Å.

Diffraction patterns were measured for solutions of 5m Na35Cl or NamixCl in D2O held within matched Ti–Zr gaskets at a pressure of 2.4, 10.0, 22.4, 27.9, or 33.8 kbar and a temperature of 423 K; an empty uncompressed Ti–Zr gasket; several empty Ti–Zr gaskets that had been recovered from different high-pressure conditions; the empty anvils with different spacings; and different sized vanadium pellets held within Ti–Zr gaskets for normalization purposes. The sample in gasket and empty gasket run times were optimized using the procedure described in Ref. 53. The data analysis followed the procedure described in Sec. II. The pressure–volume equation of state for the Ti–Zr alloy was taken from Ref. 52. At pressures greater than 33.8 kbar, the appearance of Bragg peaks indicated crystallization of the solutions.

3. Pressure calibration

For the Ti–Zr cell experiment, the pressure at the sample position was measured using a transducer. For the PE press experiment, the pressure was determined from the load applied to the anvils by using a calibration based on the F(Q) functions measured for D2O under the same experimental conditions as for the aqueous solutions.56 First, F(Q) was measured for D2O at 423 K for different applied loads, and the position of the principal peak at QPP ∼ 2.1–2.5 Å−1 was noted. At the largest applied load, crystallization was observed and a pressure of 33.8 kbar was identified from the equation of state for D2O.57 Second, the pressure dependence of QPP was measured for D2O at 423 K and 0.1, 0.5, 1.0, 1.5, or 2.0 kbar using the Ti–Zr cell. A comparison of the measured QPP values enabled a pressure of 2.4 kbar to be identified for the smallest load applied in the PE press experiment. A linear dependence between pressure and the applied load was assumed.

4. Solution density

A solution of 5m NaCl in D2O has the same ratio of NaCl to water molecules as a solution of 5.5585m NaCl in H2O, corresponding to a 24.5 wt. % solution of NaCl in light water. At pressures of up to 1 kbar, the number density of the solution in heavy water was taken to be the same as the equivalent solution in light water58 (Table I). The ambient condition value compares to ρ = 0.0950(3) Å−3 as measured using an Anton-Paar DMA 4500 M density meter. At higher pressures, ρ was taken from the equation of state found from the molecular dynamics simulations of this work (Sec. III). As shown in Fig. 1, the results are in good accord with the equations of state reported by Potter and Brown58 and Mantegazzi et al.5 The latter was parameterized for state conditions that cover this work using data for water and for aqueous solutions of 1 and 3m NaCl at temperatures of up to 673 K for pressures in the range 5–45 kbar.

TABLE I.

Temperature T and pressure p dependence of the atomic number density ρ for a solution of 5m NaCl in D2O and its ratio with the ambient conditions number density ρ0. The static dielectric constant of water ɛ is also given, as taken from Ref. 59 for pressures <10 kbar and from Ref. 13 for pressures ≥10 kbar.

SolutionWater
SetupT (K)p (kbar)ρ−3)ρ/ρ0ɛ
298 0.001 0.0954 1.000 78.46 
298 0.10 0.0957 1.003 78.85 
Ti–Zr cell 323(1) 0.10(1) 0.0947 0.993 70.27 
373(2) 0.10(1) 0.0923 0.968 55.76 
423(4) 0.10(1) 0.0895 0.939 44.30 
423(4) 0.50(1) 0.0918 0.962 45.67 
423(4) 1.00(2) 0.0938 0.983 47.14 
PE press 423(4) 2.4(5) 0.0949 0.995 50.4 
423(4) 10.0(5) 0.1061 1.112 61.0 
423(4) 22.4(5) 0.1166 1.223 71.6 
423(4) 27.9(5) 0.1209 1.267 75.3 
423(4) 33.8(5) 0.1237 1.297 79.0 
SolutionWater
SetupT (K)p (kbar)ρ−3)ρ/ρ0ɛ
298 0.001 0.0954 1.000 78.46 
298 0.10 0.0957 1.003 78.85 
Ti–Zr cell 323(1) 0.10(1) 0.0947 0.993 70.27 
373(2) 0.10(1) 0.0923 0.968 55.76 
423(4) 0.10(1) 0.0895 0.939 44.30 
423(4) 0.50(1) 0.0918 0.962 45.67 
423(4) 1.00(2) 0.0938 0.983 47.14 
PE press 423(4) 2.4(5) 0.0949 0.995 50.4 
423(4) 10.0(5) 0.1061 1.112 61.0 
423(4) 22.4(5) 0.1166 1.223 71.6 
423(4) 27.9(5) 0.1209 1.267 75.3 
423(4) 33.8(5) 0.1237 1.297 79.0 
FIG. 1.

Equation of state for a solution of 5m NaCl in D2O. The main panel shows the variation of the number density with the pressure at 423 K, and the inset shows the variation of the number density with the temperature at 0.1 kbar. The datasets are from the molecular dynamics simulations of the present work, the work of Potter and Brown,58 and the work of Mantegazzi et al.5 

FIG. 1.

Equation of state for a solution of 5m NaCl in D2O. The main panel shows the variation of the number density with the pressure at 423 K, and the inset shows the variation of the number density with the temperature at 0.1 kbar. The datasets are from the molecular dynamics simulations of the present work, the work of Potter and Brown,58 and the work of Mantegazzi et al.5 

Close modal

As indicated in Table I, the state conditions enabled access to a broad reduced-density ρ/ρ0 range of 0.94–1.30, where ρ0 is the density at 298 K and 0.001 kbar. Table I also lists ɛ for pure water, which decreases to 56% of its ambient pressure value with temperature increasing to 423 K at 0.1 kbar, before returning to this value with pressure increasing to 33.8 kbar at 423 K.

5. Scattering lengths and weighting factors

The coherent neutron scattering lengths used in the data analysis were b35Cl = 11.65(2) fm, bnatCl = 9.5770(8) fm, bmixCl = 7.365(32) fm, b37Cl = 3.08(6) fm, bNa = 3.63(2) fm, bD = 6.661(4) fm, and bO = 5.803(4) fm.60 The weighting factors for the partial structure factors in Eqs. (9) and (10) are listed in Table II (1b = 10−28 m2).

TABLE II.

Weighting factors in Eqs. (9) and (10) for the high-pressure experiments on solutions of 5m NaCl in D2O. For the Ti–Zr cell experiment, ΔFCl(Q) = FnatCl(Q) − F37Cl(Q). For the PE press experiment, ΔFCl(Q) = F35Cl(Q) − FmixCl(Q).

Coefficient
SetupA (mb)B (mb)C (mb)D (mb)
Ti–Zr cell 16.92(16) 7.37(7) 0.462(5) 0.805(4) 
PE press 11.16(10) 4.86(4) 0.305(3) 0.798(6) 
Coefficient
SetupA (mb)B (mb)C (mb)D (mb)
Ti–Zr cell 16.92(16) 7.37(7) 0.462(5) 0.805(4) 
PE press 11.16(10) 4.86(4) 0.305(3) 0.798(6) 

The measured F(Q) functions for the experiments using the Ti–Zr cell and PE press are shown in Figs. 2 and 3, respectively. The inelasticity correction P(Q) of Eq. (2) was estimated by fitting a fifth order polynomial to Qdσ/dΩ. The results show a notable shift of the first peak position in F(Q) to a larger Q-value with increasing pressure.

FIG. 2.

Total structure factors F(Q) for solutions of 5m (a) NanatCl in D2O and (b) Na37Cl in D2O. The points with the vertical error bars give the functions measured using the Ti–Zr cell, and the black curves give spline fits. The error bars are smaller than the curve thickness at most Q-values. The red curves show the simulated functions. The high-temperature functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in F(Q) at 0.1 kbar and 323 K.

FIG. 2.

Total structure factors F(Q) for solutions of 5m (a) NanatCl in D2O and (b) Na37Cl in D2O. The points with the vertical error bars give the functions measured using the Ti–Zr cell, and the black curves give spline fits. The error bars are smaller than the curve thickness at most Q-values. The red curves show the simulated functions. The high-temperature functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in F(Q) at 0.1 kbar and 323 K.

Close modal
FIG. 3.

Total structure factors F(Q) for solutions of 5m (a) Na35Cl in D2O and (b) NamixCl in D2O. The points with the vertical error bars give the functions measured using the PE press, and the black curves give spline fits. The error bars are smaller than the curve thickness at most Q-values. The red curves show the simulated functions. The high-pressure functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in F(Q) at 2.4 kbar and 423 K.

FIG. 3.

Total structure factors F(Q) for solutions of 5m (a) Na35Cl in D2O and (b) NamixCl in D2O. The points with the vertical error bars give the functions measured using the PE press, and the black curves give spline fits. The error bars are smaller than the curve thickness at most Q-values. The red curves show the simulated functions. The high-pressure functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in F(Q) at 2.4 kbar and 423 K.

Close modal

The G(r) functions for the solutions of 5m NanatCl or Na35Cl in D2O are shown in Fig. 4. The first peak at 0.95(1) Å is identified with the intra-molecular O–D distances within water molecules, and integration of this peak gives an O–D coordination number n̄OD = 2.02(2). The O–D bond distance is shorter than the value of 0.985(5) Å measured for pure D2O using the method of neutron diffraction with oxygen isotope substitution, i.e., in experiments designed to eliminate the main effects of inelastic scattering on the bond length.48,61 A foreshortened O–D bond length is typical of reactor-based neutron diffraction experiments on water or aqueous solutions, an artifact that originates from the difficulty in accurately correcting for the effects of recoil in neutron scattering events.35,62 Such bond-length issues are absent for ΔGCl(r) because the main effects of inelastic scattering are eliminated when forming the first-difference function ΔFCl(Q) (Sec. II).

FIG. 4.

Total pair-distribution functions G(r) for solutions of 5m (a) Na35Cl in D2O and (b) NanatCl in D2O. The black curves were obtained by Fourier transforming the measured F(Q) functions given by the black curves in Figs. 3(a) and 2(a), respectively, and setting the low-r oscillations (chained blue curves) to the calculated limit G(r → 0). The red curves in (a) and (b) were obtained by Fourier transforming the simulated F(Q) functions given by the red curves in Figs. 3(a) and 2(a), respectively, after truncation at a node in F(Q) (Qmax = 20 Å−1). The green broken vertical line marks the O–D bond distance in the TIP4P/2005 model.

FIG. 4.

Total pair-distribution functions G(r) for solutions of 5m (a) Na35Cl in D2O and (b) NanatCl in D2O. The black curves were obtained by Fourier transforming the measured F(Q) functions given by the black curves in Figs. 3(a) and 2(a), respectively, and setting the low-r oscillations (chained blue curves) to the calculated limit G(r → 0). The red curves in (a) and (b) were obtained by Fourier transforming the simulated F(Q) functions given by the red curves in Figs. 3(a) and 2(a), respectively, after truncation at a node in F(Q) (Qmax = 20 Å−1). The green broken vertical line marks the O–D bond distance in the TIP4P/2005 model.

Close modal

The molecular dynamics simulations reproduce the main features in the measured F(Q) functions (Figs. 2 and 3). In accordance with the use of a rigid model for the water molecules, the large-Q oscillations are, however, more pronounced than found by the experiment, which leads [via Eq. (4)] to a sharper first peak in G(r) and to oscillations at smaller r-values about the G(r → 0) limit (Fig. 4).

The ΔFCl(Q) functions obtained directly from the measured differential scattering cross sections are shown in Fig. 5. The associated ΔFCl(Q) functions [Eq. (8)] did not show a Q-dependent slope, indicating that the (small) light water content of the solutions was balanced.63,64 The functions show a measurable contrast between the total structure factors, even for the experiment using the PE press where the volume of sample illuminated by the incident neutron beam was <7% of the volume illuminated in the Ti–Zr cell experiment. With increasing pressure at 423 K, the first peak in ΔFCl(Q) at ≃2.1 Å−1 sharpens and moves toward a larger Q-value.

FIG. 5.

Difference functions ΔFCl(Q) for solutions of 5m NaCl in D2O as measured using (a) the PE press or (b) the Ti–Zr cell. The points with the vertical error bars show the measured data points, and the black curves show spline fits to which a Lorch function38 has been applied with Qmax set at either (a) 9.75 Å−1 or (b) 21 Å−1. The red curves show the results from the simulations. In (a) and (b), the high-pressure and high-temperature datasets are offset vertically for clarity of presentation, respectively. The green broken vertical line gives the position of the first peak in ΔFCl(Q) at 0.1 kbar and 323 K.

FIG. 5.

Difference functions ΔFCl(Q) for solutions of 5m NaCl in D2O as measured using (a) the PE press or (b) the Ti–Zr cell. The points with the vertical error bars show the measured data points, and the black curves show spline fits to which a Lorch function38 has been applied with Qmax set at either (a) 9.75 Å−1 or (b) 21 Å−1. The red curves show the results from the simulations. In (a) and (b), the high-pressure and high-temperature datasets are offset vertically for clarity of presentation, respectively. The green broken vertical line gives the position of the first peak in ΔFCl(Q) at 0.1 kbar and 323 K.

Close modal

The difference functions ΔGCl(r) measured in the experiment using the Ti–Zr cell are shown in Fig. 6. The large cutoff value Qmax = 21 Å−1 of the Lorch modification function leads to a distinct first peak in ΔGCl(r). This peak is attributed to Cl–D correlations37,65,66 and gives the bond lengths and coordination numbers listed in Table III. The second peak in ΔGCl(r) has contributions from overlapping Cl–O and Cl–D correlations (Fig. S1). If the peak position is assigned to the nearest-neighbor Cl–O distance, the calculated D–Cl–O bond angle is consistent with approximately linear Cl⋯D–O hydrogen bonds, within the experimental uncertainty in the Cl–D and Cl–O distances.

FIG. 6.

Difference functions ΔGCl(r) for solutions of 5m NaCl in D2O. The black curves were obtained by Fourier transforming the measured ΔFCl(Q) functions given by the black curves in Fig. 5(b) and setting the low-r oscillations (blue chained curves) to the calculated limit ΔGCl(r → 0). The red curves were obtained by Fourier transforming the simulated ΔFCl(Q) functions given by the red curves in Fig. 5(b) after the application of a Lorch modification function with Qmax = 21 Å−1 as per the measured datasets. The high temperature functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in ΔGCl(r) at 0.1 kbar and 323 K.

FIG. 6.

Difference functions ΔGCl(r) for solutions of 5m NaCl in D2O. The black curves were obtained by Fourier transforming the measured ΔFCl(Q) functions given by the black curves in Fig. 5(b) and setting the low-r oscillations (blue chained curves) to the calculated limit ΔGCl(r → 0). The red curves were obtained by Fourier transforming the simulated ΔFCl(Q) functions given by the red curves in Fig. 5(b) after the application of a Lorch modification function with Qmax = 21 Å−1 as per the measured datasets. The high temperature functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in ΔGCl(r) at 0.1 kbar and 323 K.

Close modal
TABLE III.

Pressure and temperature dependence of the measured Cl–D bond length rClD, position of the second peak r2 in ΔGCl(r), and Cl–D coordination number n̄ClD for a solution of 5m NaCl in D2O. rmin and rmax refer the limits used to obtain n̄ClD [Eq. (6)]. The values are compared to those found for a solution of 5.32m NaCl in D2O under ambient conditions (0.001 kbar and 298 K).37,65,66

p (kbar)T (K)rClD (Å)r2 (Å)n̄ClDrmin (Å)rmax (Å)
0.001 298 2.26(3) 3.20(5) 5.5(4) 1.79 2.87 
0.10(1) 323(1) 2.24(1) 3.52(9) 5.7(3) 1.66 2.82 
373(2) 2.27(1) 3.33(9) 5.5(3) 1.72 2.82 
423(4) 2.28(1) 3.39(9) 4.9(3) 1.72 2.76 
0.50(1) 423(4) 2.24(1) 3.25(9) 5.4(3) 1.66 2.82 
1.00(2)  2.24(1) 3.27(9) 5.2(5) 1.72 2.76 
2.4(5)  2.20(2) 3.18(9) 6.1(7) 1.72 2.76 
10.0(5)  2.22(2) 3.13(9) 8.2(7) 1.66 2.82 
22.4(5)  2.21(2) 3.23(9) 5.2(7) 1.47 2.58 
27.9(5)  2.21(4) 3.15(9) 3.4(7) 1.47 2.52 
33.8(5)  2.22(2) 3.19(9) 4.4(7) 1.47 2.52 
p (kbar)T (K)rClD (Å)r2 (Å)n̄ClDrmin (Å)rmax (Å)
0.001 298 2.26(3) 3.20(5) 5.5(4) 1.79 2.87 
0.10(1) 323(1) 2.24(1) 3.52(9) 5.7(3) 1.66 2.82 
373(2) 2.27(1) 3.33(9) 5.5(3) 1.72 2.82 
423(4) 2.28(1) 3.39(9) 4.9(3) 1.72 2.76 
0.50(1) 423(4) 2.24(1) 3.25(9) 5.4(3) 1.66 2.82 
1.00(2)  2.24(1) 3.27(9) 5.2(5) 1.72 2.76 
2.4(5)  2.20(2) 3.18(9) 6.1(7) 1.72 2.76 
10.0(5)  2.22(2) 3.13(9) 8.2(7) 1.66 2.82 
22.4(5)  2.21(2) 3.23(9) 5.2(7) 1.47 2.58 
27.9(5)  2.21(4) 3.15(9) 3.4(7) 1.47 2.52 
33.8(5)  2.22(2) 3.19(9) 4.4(7) 1.47 2.52 

The difference functions ΔGCl(r) measured in the experiment using the PE press are shown in Fig. 7. Here, the small cutoff value Qmax = 9.75 Å−1 of the Lorch modification function, which was used to avoid the Fourier transform artifacts that originate from high-Q noise in ΔFCl(Q), leads to considerable overlap of the first and second peaks in ΔGCl(r). For this reason, the n̄ClD coordination numbers (Table III) are less reliable than found in the Ti–Zr cell experiment. Nevertheless, the first and second peak positions in ΔGCl(r) are consistent with approximately linear Cl⋯D–O hydrogen bonds, within the experimental uncertainty in the Cl–D and Cl–O distances. The nature of these hydrogen bonds will be discussed in Sec. VI A.

FIG. 7.

Difference functions ΔGCl(r) for solutions of 5m NaCl in D2O. The black curves were obtained by Fourier transforming the measured ΔFCl(Q) functions given by the black curves in Fig. 5(a) and setting the low-r oscillations (blue chained curves) to the calculated limit ΔGCl(r → 0). The red broken and solid curves were obtained by Fourier transforming the simulated ΔFCl(Q) functions given by the red curves in Fig. 5(a) after the application of a Lorch modification function with either (i) Qmax = 9.75 Å−1 as per the measured datasets or (ii) Qmax = 21 Å−1 to provide a comparison with the datasets shown in Fig. 6, respectively. The high-pressure functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in ΔGCl(r) at 2.4 kbar and 423 K.

FIG. 7.

Difference functions ΔGCl(r) for solutions of 5m NaCl in D2O. The black curves were obtained by Fourier transforming the measured ΔFCl(Q) functions given by the black curves in Fig. 5(a) and setting the low-r oscillations (blue chained curves) to the calculated limit ΔGCl(r → 0). The red broken and solid curves were obtained by Fourier transforming the simulated ΔFCl(Q) functions given by the red curves in Fig. 5(a) after the application of a Lorch modification function with either (i) Qmax = 9.75 Å−1 as per the measured datasets or (ii) Qmax = 21 Å−1 to provide a comparison with the datasets shown in Fig. 6, respectively. The high-pressure functions are offset vertically for clarity of presentation. The green broken vertical line gives the position of the first peak in ΔGCl(r) at 2.4 kbar and 423 K.

Close modal

At 0.1 kbar, the first peak in ΔGCl(r) broadens with increasing temperature as the Cl–D bond length increases, the Cl–D coordination number decreases (Table III), and the height of the function at the first minimum (at ≃2.8 Å) increases [Fig. 8(a)]. This increase in height at the first minimum originates primarily from broadening of the first peak in gClD(r) (see Sec. V C).

FIG. 8.

Measured dependence of the difference function ΔGCl(r) − ΔGCl(0) on the (a) temperature at 0.1 kbar vs (b) pressure at 423 K for a solution of 5m NaCl in D2O. The results were obtained using the Ti–Zr cell.

FIG. 8.

Measured dependence of the difference function ΔGCl(r) − ΔGCl(0) on the (a) temperature at 0.1 kbar vs (b) pressure at 423 K for a solution of 5m NaCl in D2O. The results were obtained using the Ti–Zr cell.

Close modal

At 423 K, the relative height of the first and second peaks in ΔGCl(r) changes with increasing pressure, and the first peak broadens [Figs. 7 and 8(b)]. There is little change to the Cl–D bond length as taken from the first peak or shoulder position in ΔGCl(r). It is difficult to assess the effect of pressure on the Cl–D coordination number because of overlap between the first and second peaks in ΔGCl(r).

The molecular dynamics simulations reproduce the main features in the measured ΔFCl(Q) functions (Fig. 5), including the sharpening of the first peak and its shift in position to a larger Q-value with increasing pressure. The simulations also reproduce the main features in the corresponding real-space functions ΔGCl(r) (Figs. 6 and 7), although the first and second peaks are, in general, sharper than found by the experiment.

The dependence of the simulated gClD(r) and gClO(r) functions on temperature and pressure is shown in Fig. 9. As in the experiment, the first peak in gClD(r) broadens with increasing temperature. Large-r oscillations develop with increasing pressure, which are particularly notable for gClO(r).

FIG. 9.

Dependence of the simulated partial pair-distribution functions gClD(r) and gClO(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K.

FIG. 9.

Dependence of the simulated partial pair-distribution functions gClD(r) and gClO(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K.

Close modal

The n̄ClD values obtained by integrating to the first minimum in gClD(r) are summarized in Fig. 10. As in the experiment, n̄ClD decreases with increasing temperature, although the simulated coordination numbers are larger than the measured values. The Cl–D coordination number increases with pressure, where the peak overlap makes it difficult to extract reliable values from the experiment. The n̄ClD and n̄ClO values obtained by integrating to the first minimum in gClD(r) and gClO(r), respectively, are similar for the different temperatures but diverge with increasing pressure, where n̄ClO becomes substantially larger than n̄ClD (Fig. S2).

FIG. 10.

Dependence of the coordination numbers n̄NaO (black open squares), n̄ClD (red open triangles), and n̄NaCl=n̄ClNa (blue open circles) obtained from the simulations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. Each n̄αβ value was obtained by integrating to the first minimum in the corresponding gαβ(r) function. The total coordination numbers of the Na+ ion n̄Na=n̄NaO+n̄NaCl (black solid squares) and Cl ion n̄Cl=n̄ClD+n̄ClNa (red solid triangles) are also shown. The n̄ClD values are compared to those measured in the Ti–Zr cell experiment (magenta open triangles with the vertical error bars) (Table III).

FIG. 10.

Dependence of the coordination numbers n̄NaO (black open squares), n̄ClD (red open triangles), and n̄NaCl=n̄ClNa (blue open circles) obtained from the simulations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. Each n̄αβ value was obtained by integrating to the first minimum in the corresponding gαβ(r) function. The total coordination numbers of the Na+ ion n̄Na=n̄NaO+n̄NaCl (black solid squares) and Cl ion n̄Cl=n̄ClD+n̄ClNa (red solid triangles) are also shown. The n̄ClD values are compared to those measured in the Ti–Zr cell experiment (magenta open triangles with the vertical error bars) (Table III).

Close modal

A similar escalation in n̄ClO with pressure is reported from x-ray diffraction experiments on a 3m solution of NaCl in (light) water at 300 K, where an analysis using empirical potential structure refinement showed n̄ClO = 5.9(1.1) at ambient pressure vs n̄ClO = 16.2(1.6) at 1.7 GPa28 and from other molecular dynamics simulations on aqueous NaCl solutions.67 The discrepancy between the high-pressure n̄ClD and n̄ClO values found in the present work indicates that many of these oxygen atoms are in water molecules to which the chloride ion does not form hydrogen bonds.

Figure 11 shows the Cl–D–D bond-angle distributions where both deuterium atoms belong to the same water molecule [Fig. 12(a)]. There is a peak at ∼145°, indicating that water molecules in the primary solvation shell of the chloride ion are orientated in order to facilitate hydrogen bonding with species in the second solvation shell. Figure 11 also shows the Cl–O–D bond-angle distributions where O and D belong to the same water molecule and D is the nearest-neighbor to the chloride ion [Fig. 12(b)]. At ambient conditions, there is a maximum around zero with most of the angles less than 35°, which indicates that water molecules in the primary solvation shell of the chloride ion are orientated to form single hydrogen bonds. The shoulder that grows around 45°–50° with increasing pressure originates from non-hydrogen-bonded configurations and will be discussed in Sec. VI A.

FIG. 11.

Dependence of the simulated Cl–D–D, Cl–O–D, and Na–O–M bond-angle distributions on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. The insets zoom into the Cl–O–D bond-angle distributions, where the arrow points to the growth of a shoulder with increasing pressure. The distributions were calculated using Cl–D, Cl–O, and Na–O cutoff distances of 2.98, 3.92, and 3.2 Å, respectively.

FIG. 11.

Dependence of the simulated Cl–D–D, Cl–O–D, and Na–O–M bond-angle distributions on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. The insets zoom into the Cl–O–D bond-angle distributions, where the arrow points to the growth of a shoulder with increasing pressure. The distributions were calculated using Cl–D, Cl–O, and Na–O cutoff distances of 2.98, 3.92, and 3.2 Å, respectively.

Close modal
FIG. 12.

Schematic showing the (a) Cl–D–D, (b) Cl–O–D, and (c) Na–O–M bond angles.

FIG. 12.

Schematic showing the (a) Cl–D–D, (b) Cl–O–D, and (c) Na–O–M bond angles.

Close modal

The dependence of the simulated gNaO(r) and gNaD(r) functions on temperature and pressure is shown in Fig. 13, and the corresponding running coordination numbers are shown in Fig. S3. The n̄NaO and n̄NaD/2 values obtained by integrating to the first minimum in gNaO(r) and gNaD(r), respectively, are in accord at all state conditions, i.e., most of the deuterium atoms belong to the same water molecules as the nearest-neighbor oxygen atoms. The coordination number n̄NaO decreases with increasing temperature but increases with pressure (Fig. 10).

FIG. 13.

Dependence of the simulated partial pair-distribution functions gNaO(r) and gNaD(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K.

FIG. 13.

Dependence of the simulated partial pair-distribution functions gNaO(r) and gNaD(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K.

Close modal

Figure 11 shows the Na–O–M bond-angle distributions where O–M represents a line along the dipole moment of the water molecule [Fig. 12(c)]. These distributions show a maximum around 180° for all temperatures at 0 kbar, which indicates that the dipole moment of a solvated water molecule points directly toward the Na+ ion. The maximum shifts to a smaller angle with increasing pressure, indicating the development of a finite tilt angle.

The temperature and pressure dependence of the simulated gNaCl(r), gNaNa(r), and gClCl(r) functions is shown in Fig. 14, and the Na–Cl running coordination number is shown in Fig. S4. Contact Na–Cl pairs occur at a distance of ≃2.8 Å. The Na–Cl coordination number, obtained by integrating to the first minimum in gNaCl(r), does not change markedly with the state conditions (Fig. 10). The fraction of Na+ ions with zero, one, two, or three nearest-neighbor Cl ions is given in Fig. 15. The majority of sodium ions do not form contact ion pairs. Of those that do, most are in charge-neutral NaCl0 units, the fraction of which is more sensitive to changing temperature than to changing pressure (the fraction increases by ∼120% as the temperature increases from 298 to 423 K and decreases by ∼21% as the pressure increases from 0 to 33.8 kbar). For the higher coordinated ion pairs, the Na–Cl–Na bond-angle distribution has broad peaks at ∼85° and 110° under ambient conditions, which become less pronounced with increasing temperature and show little variation with pressure (Fig. S5). The Cl–Na–Cl bond-angle distribution has broad peaks at ∼90° and 180° under ambient conditions, which show some variation with the state conditions (Fig. S5).

FIG. 14.

Dependence of the simulated partial pair-distribution functions gNaCl(r), gNaNa(r), and gClCl(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. At 0 kbar, the arrows point to the first two peaks in gNaNa(r) at about 3.7 and 4.4 Å.

FIG. 14.

Dependence of the simulated partial pair-distribution functions gNaCl(r), gNaNa(r), and gClCl(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. At 0 kbar, the arrows point to the first two peaks in gNaNa(r) at about 3.7 and 4.4 Å.

Close modal
FIG. 15.

Dependence of the fraction of Na+ ions with zero, one, two, or three nearest-neighbor Cl ions, as obtained from the simulations, on the (a) temperature at 0 kbar vs (b) pressure at 423 K. A cutoff distance of 3.5 Å was used in the calculations.

FIG. 15.

Dependence of the fraction of Na+ ions with zero, one, two, or three nearest-neighbor Cl ions, as obtained from the simulations, on the (a) temperature at 0 kbar vs (b) pressure at 423 K. A cutoff distance of 3.5 Å was used in the calculations.

Close modal

The Na–Na pair-distribution function has a split nearest-neighbor peak under ambient conditions (Fig. 14). To investigate the origin of this feature, the percentage of Na+ ion pairs that have common D2O molecules and/or Cl ions was calculated for Na–Na cutoff distances of 4.06 and 4.90 Å, which correspond to the position of the minimum after the first and second peaks in gNaNa(r), respectively, for the solution at 298 K and 0 kbar. Figure 16 shows that, under ambient conditions, the first peak originates primarily from Na+ pairs that share two water molecules in edge-sharing conformations [Fig. 17(a)] with an Na–O–Na bond angle of ∼98° (Fig. S6). The second peak originates primarily from Na+ pairs that share a single water molecule in corner-sharing conformations with an Na–O–Na bond angle of ∼125° (Fig. S6). With increasing temperature, the fraction of edge-sharing motifs decreases as the first and second peaks in gNaNa(r) merge. With increasing pressure, there is a small increase in the fraction of edge-sharing motifs as the first peak in gNaNa(r) moves to a smaller distance (Fig. 18).

FIG. 16.

Temperature dependence of the fraction of Na+ ion pairs obtained from the simulations that share the listed numbers of water molecules and/or chloride ions. The sodium ion separation is either (a) ≤4.06 Å or (b) ≤4.9 Å.

FIG. 16.

Temperature dependence of the fraction of Na+ ion pairs obtained from the simulations that share the listed numbers of water molecules and/or chloride ions. The sodium ion separation is either (a) ≤4.06 Å or (b) ≤4.9 Å.

Close modal
FIG. 17.

Schematics showing (a) the bridging of two Na+ ions by two water molecules to give an edge-sharing configuration and (b) the type of configuration that contributes toward the peak in the Cl–O–Cl bond-angle distribution under ambient conditions (Fig. 19). In (a), the D atoms lie above and below the plane made by the Na+ ions and O atoms.

FIG. 17.

Schematics showing (a) the bridging of two Na+ ions by two water molecules to give an edge-sharing configuration and (b) the type of configuration that contributes toward the peak in the Cl–O–Cl bond-angle distribution under ambient conditions (Fig. 19). In (a), the D atoms lie above and below the plane made by the Na+ ions and O atoms.

Close modal
FIG. 18.

Pressure dependence of the fraction of Na+ ion pairs obtained from the simulations that share the listed numbers of water molecules and/or chloride ions. The sodium ion separation is either (a) ≤4.06 Å or (b) ≤4.9 Å.

FIG. 18.

Pressure dependence of the fraction of Na+ ion pairs obtained from the simulations that share the listed numbers of water molecules and/or chloride ions. The sodium ion separation is either (a) ≤4.06 Å or (b) ≤4.9 Å.

Close modal

Figure 19 shows a sharp peak in the Cl–O–Cl bond-angle distribution at ∼104° under ambient conditions. The peak broadens with increasing temperature, and broadens and shifts position toward a smaller angle with increasing pressure. Figure 20 gives the associated bond distance resolved bond-angle map, i.e., it shows the dependence of the Cl–O–Cl bond angle on the Cl–O distance. The map shows that the peak at ∼104° originates from Cl–O distances around 3.1 Å, which corresponds to the first peak position in gClO(r) (Fig. 9). These observations indicate the presence of solvent-separated chloride ion pairs, in which two Cl ions share a single water molecule: The Cl⋯O–D bonds are almost linear (Sec. V C) and the intra-molecular D–O–D bond angle is 104.52° [Fig. 17(b)]. Such configurations are the most abundant at low pressure but become less prevalent as the pressure increases and other types of solvent-separated Cl–Cl pairs become more numerous (Fig. 21).

FIG. 19.

Dependence of the simulated O–O–O, O–Cl–O, and Cl–O–Cl bond-angle distributions on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. The distributions were calculated using O–O and Cl–O cutoff distances of 4.3 and 3.92 Å, respectively. The arrows mark a Cl–O–Cl bond angle of 104.52°.

FIG. 19.

Dependence of the simulated O–O–O, O–Cl–O, and Cl–O–Cl bond-angle distributions on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. The distributions were calculated using O–O and Cl–O cutoff distances of 4.3 and 3.92 Å, respectively. The arrows mark a Cl–O–Cl bond angle of 104.52°.

Close modal
FIG. 20.

Map showing the simulated dependence of the Cl–O–Cl bond angle on the Cl–O distance. (a) T = 298 K, p = 0 kbar; (b) T = 323 K, p = 0 kbar; (c) T = 373 K, p = 0 kbar; (d) T = 423 K, p = 0 kbar; (e) T = 423 K, p = 10 kbar; and (f) T = 423 K, p = 33.8 kbar. The multiplicities are normalized to give a maximum value of unity.

FIG. 20.

Map showing the simulated dependence of the Cl–O–Cl bond angle on the Cl–O distance. (a) T = 298 K, p = 0 kbar; (b) T = 323 K, p = 0 kbar; (c) T = 373 K, p = 0 kbar; (d) T = 423 K, p = 0 kbar; (e) T = 423 K, p = 10 kbar; and (f) T = 423 K, p = 33.8 kbar. The multiplicities are normalized to give a maximum value of unity.

Close modal
FIG. 21.

Dependence of the fraction of chloride ion pairs obtained from the simulations that share the listed numbers of water molecules on the (a) temperature at 0 kbar vs (b) pressure at 423 K.

FIG. 21.

Dependence of the fraction of chloride ion pairs obtained from the simulations that share the listed numbers of water molecules on the (a) temperature at 0 kbar vs (b) pressure at 423 K.

Close modal

The temperature and pressure dependence of the simulated gOD(r), gDD(r), and gOO(r) functions is shown in Fig. 22, and the O–O running coordination number is shown in Fig. S7. The first peak in gOD(r) broadens with increasing temperature but changes little with increasing pressure. The first peak in gOO(r) also broadens with increasing temperature but sharpens when p ≥ 10 kbar as the nearest-neighbor O–O coordination number increases and pronounced large-r oscillations develop. Large-r oscillations also develop in gClO(r) as the pressure is increased (Fig. 9).

FIG. 22.

Dependence of the simulated partial pair-distribution functions gOD(r), gDD(r), and gOO(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. The sharp peak in gDD(r) at 1.514 Å originates from the intra-molecular correlations within rigid water molecules.

FIG. 22.

Dependence of the simulated partial pair-distribution functions gOD(r), gDD(r), and gOO(r) on the (left-hand column) temperature at 0 kbar vs (right-hand column) pressure at 423 K. The sharp peak in gDD(r) at 1.514 Å originates from the intra-molecular correlations within rigid water molecules.

Close modal

The O–O–D bond-angle distributions are shown in Fig. 23 where D and the central O atom reside in the same water molecule (Fig. 24). At ambient conditions, these distributions show a maximum around zero with most of the angles less than 30°, a typical cutoff value for the acceptor hydrogen bonds in water.68 The shoulder at ∼52° corresponds to non-hydrogen-bonded configurations. The origin of these features will be discussed in Sec. VI A.

FIG. 23.

Dependence of the simulated O–O–D bond-angle distribution on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The insets zoom into the region of a shoulder at ∼52°. The distributions were calculated using an O–O cutoff distance of 4.3 Å.

FIG. 23.

Dependence of the simulated O–O–D bond-angle distribution on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The insets zoom into the region of a shoulder at ∼52°. The distributions were calculated using an O–O cutoff distance of 4.3 Å.

Close modal
FIG. 24.

Schematic showing the O–O–D bond angle. The angle is either large or small depending on whether the water molecule is accepting or donating a hydrogen bond, respectively.

FIG. 24.

Schematic showing the O–O–D bond angle. The angle is either large or small depending on whether the water molecule is accepting or donating a hydrogen bond, respectively.

Close modal

Figure 19 shows that the O–O–O and O–Cl–O bond-angle distributions take similar forms at each of the investigated state points, which suggests that the chloride ions attempt to fit into the hydrogen bonded network of the water molecules.

The Na+ and Cl ions have filled electron shells but differ in their electric charge and ion radius. The ion field strength Z/rion2, where Z = 1 and rion is the ion radius, is larger for Na+ compared to Cl, i.e., 0.961 Å−2 compared to 0.305 Å−2, respectively, using the Shannon69 radius for sixfold coordinated ions. The charge within a water molecule is non-uniformly distributed such that each hydrogen atom has a partial charge of +δe and the oxygen atom is associated with a partial charge of −2δe, and the ion–water orientation is different for cations [Fig. 12(c)] compared to anions [Fig. 12(a)]. Therefore, in dilute solution where the ions are well separated, the residence time of a water molecule in the first solvation shell of an ion is expected to be longer for Na+ vs Cl. In concentrated solution, however, greater structural complexity originates from contact and solvent-separated ion pairing, the detail of which is dependent on the state conditions. It is therefore informative to investigate the effect of this complexity on the ion–water dynamics.

The residence time τ of a water molecule in the first solvation shell of Na+ or Cl was assessed from the simulations by using the cage-correlation function70,71

Cout(t)=Θ1niout0,t.
(11)

Here, a water molecule is deemed to be within the cage of Na+ if the Na–O distance is within a cutoff distance defined by the first minimum in gNaO(r) or the cage of Cl if the Cl–D distance is within a cutoff distance defined by the first minimum in gClD(r). In Eq. (11), niout(0,t) is the number of water molecules that exit the cage of ion i during time t; the Heaviside step function Θ takes the value unity if its argument is greater than zero or the value zero otherwise; and the angular brackets denote an average over all ions and all origins in time. As defined, there is a finite contribution to Eq. (11) if niout(0,t)<1, so Cout(t) measures the average rate at which a single water molecule exits a cage. The converse process is represented by the function Cin(t), which is also defined by Eq. (11) but with niout(0,t) replaced by niin(0,t), the number of water molecules that enter the cage of ion i during time t. By going backward in time from the end of the simulation, all the exits become entrances and vice versa. It follows that Cin(t) = Cout(t).70 

The temperature and pressure dependence of the Cout(t) functions for the Na+ and Cl ions is shown in Figs. S8 and S9, respectively. The contribution from short-time vibrational motion of the water molecules is small such that each function is dominated by the longer-time diffusive motion of the water molecules as represented by Cout(t)=Aexpt/τout, where A is an amplitude. The residence time τout is defined as the time required for the function to decay to 1/e of its original value and was extracted from the integral τout=0dttCout(t)/0dtCout(t).

Figure 25 shows the temperature and pressure dependence of the residence time of water molecules in the cages of the Na+ and Cl ions, denoted by τNaout and τClout, respectively. At ambient conditions, the residence time for Na+ is larger than for Cl, and the value for the chloride ion is consistent with the limit of ≤100 ps found for concentrated aqueous solutions from experiments using incoherent quasi-elastic neutron scattering.72,73 The residence time for Na+ shows a significant decrease with increasing temperature, approaching the value for the chloride ion. In comparison, there is little change to either residence time with increasing pressure.

FIG. 25.

Dependence of the residence time ταout (α = Na or Cl) obtained from the simulations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The curves are drawn as guides for the eye. The insets depict the dependence of lnταout on (a) 1/T at 0 kbar or (b) p at 423 K. The curves are straight-line fits.

FIG. 25.

Dependence of the residence time ταout (α = Na or Cl) obtained from the simulations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The curves are drawn as guides for the eye. The insets depict the dependence of lnταout on (a) 1/T at 0 kbar or (b) p at 423 K. The curves are straight-line fits.

Close modal

Under isobaric conditions, the residence time shows an Arrhenius temperature dependence τout=τ0outexp(Ea/RT), where R is the molar gas constant [inset of Fig. 25(a)]. The activation energy Ea for the exchange of water molecules is larger for Na+ than for Cl (Table IV). Under isothermal conditions, the residence times do not show a strong dependence on pressure. From transition state theory, it is expected that

lnτoutpT=ΔVRT,
(12)

where the volume of activation ΔV is the difference in partial molar volumes of the reactants and transition state.74 An approximately linear dependence between ln τout and p is found for each ion [inset of Fig. 25(b)], leading to a small activation volume (Table IV).

TABLE IV.

The activation energy Ea at 0 kbar and the activation volume ΔV at 423 K for the exchange of water molecules at the Na+ and Cl ions.

IonEa (kJ mol−1)ΔV (cm3 mol−1)
Na+ 23.9(1.5) 0.15(3) 
Cl 17.1(0.7) −0.30(4) 
IonEa (kJ mol−1)ΔV (cm3 mol−1)
Na+ 23.9(1.5) 0.15(3) 
Cl 17.1(0.7) −0.30(4) 

The self-diffusion coefficient of a given species in the simulations was found from the long-time behavior of the mean-square displacement for that species using the relation

D=limt16t|ri(t)ri(0)|2,
(13)

where ri(t) is the position of particle i at time t. Figure 26 shows the temperature and pressure dependence of the self-diffusion coefficients for the Na+ and Cl ions and water molecules.

FIG. 26.

Dependence of the self-diffusion coefficients Dα (α = Na, Cl, or D2O) obtained from the simulations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The curves are drawn as guides for the eye. The insets show the dependence of ln Dα on (a) 1/T at 0 kbar or (b) p at 423 K. The curves are straight-line fits.

FIG. 26.

Dependence of the self-diffusion coefficients Dα (α = Na, Cl, or D2O) obtained from the simulations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The curves are drawn as guides for the eye. The insets show the dependence of ln Dα on (a) 1/T at 0 kbar or (b) p at 423 K. The curves are straight-line fits.

Close modal

Under isobaric conditions, the diffusion coefficients show an Arrhenius temperature dependence D = D0 exp(−Ea/RT) [inset of Fig. 26(a)]. The activation energies for the different species take comparable values (Table V) and are similar to the values found for the exchange of water molecules at the Na+ and Cl ions (Table IV). Under isothermal conditions, it was assumed that the ions or water molecules move by a jump diffusion mechanism75 involving a transition state, such that the pressure dependence of the mobility is associated with a volume of activation, as in Eq. (12). The data were analyzed accordingly, using  ln D/∂p|T = −ΔV/RT [inset of Fig. 26(b)], and the emergent activation volumes take similar values for all the species (Table V).

TABLE V.

The activation energy Ea at 0 kbar and the activation volume ΔV at 423 K for self-diffusion of the Na+ and Cl ions and water molecules.

SpeciesEa (kJ mol−1)ΔV (cm3 mol−1)
Na+ 22.9(1.0) 1.23(4) 
Cl 21.2(1.2) 1.41(7) 
D219.3(5) 1.51(7) 
SpeciesEa (kJ mol−1)ΔV (cm3 mol−1)
Na+ 22.9(1.0) 1.23(4) 
Cl 21.2(1.2) 1.41(7) 
D219.3(5) 1.51(7) 

Figure 27 maps the dependence of the Cl–O–D bond angle on the Cl–O distance. At ambient conditions, there is a peak near 3.1 Å and 8° that broadens as the state conditions are changed, extending to a maximum angle around 35°. The peak position shifts to ∼3.0 Å and 10° with increasing pressure as a second less-intense feature develops around 3.6 Å and 50°. From this breakdown, it follows that the peak at small angles in the Cl–O–D bond-angle distribution (Fig. 11) corresponds to short Cl–O distances associated with the first peak in gClO(r) and therefore originates from the formation of Cl⋯D–O hydrogen bonds [Fig. 12(b)]. With increasing pressure, the broad shoulder that develops around 50° corresponds to longer Cl–O distances and, therefore, to non-hydrogen-bonded Cl–water conformations.

FIG. 27.

Map showing the simulated dependence of the Cl–O–D bond angle [Fig. 12(b)] on the Cl–O distance. (a) T = 298 K, p = 0 kbar; (b) T = 323 K, p = 0 kbar; (c) T = 373 K, p = 0 kbar; (d) T = 423 K, p = 0 kbar; (e) T = 423 K, p = 10 kbar; and (f) T = 423 K, p = 33.8 kbar. The multiplicities are normalized to give a maximum value of unity.

FIG. 27.

Map showing the simulated dependence of the Cl–O–D bond angle [Fig. 12(b)] on the Cl–O distance. (a) T = 298 K, p = 0 kbar; (b) T = 323 K, p = 0 kbar; (c) T = 373 K, p = 0 kbar; (d) T = 423 K, p = 0 kbar; (e) T = 423 K, p = 10 kbar; and (f) T = 423 K, p = 33.8 kbar. The multiplicities are normalized to give a maximum value of unity.

Close modal

Figure 28 shows the contributions to the coordination number n̄ClO from hydrogen-bonded vs non-hydrogen-bonded configurations. From Fig. 27, oxygen and chlorine are deemed to be hydrogen bonded if rClO < 3.6 Å and the Cl–O–D bond angle <35. The n̄ClO values for hydrogen-bonded neighbors are in accord with the n̄ClD values obtained from the diffraction experiments using the Ti–Zr cell. The large increase in n̄ClO with pressure originates primarily from non-hydrogen-bonded Cl–water conformations at larger distances.

FIG. 28.

Dependence of the coordination numbers n̄ClO and n̄OO obtained from the simulations for hydrogen-bonded vs non-hydrogen-bonded configurations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The results for the aqueous solution of 5m NaCl in D2O are compared to those obtained for pure D2O. The coordination numbers for hydrogen-bonded configurations are given by the solid symbols with solid curves, and the coordination numbers for non-hydrogen-bonded configurations are given by the open symbols with broken curves. Also shown are the measured n̄ClD values from the Ti–Zr cell experiment (magenta solid downward triangles with vertical error bars) (Table III).

FIG. 28.

Dependence of the coordination numbers n̄ClO and n̄OO obtained from the simulations for hydrogen-bonded vs non-hydrogen-bonded configurations on the (a) temperature at 0 kbar vs (b) pressure at 423 K. The results for the aqueous solution of 5m NaCl in D2O are compared to those obtained for pure D2O. The coordination numbers for hydrogen-bonded configurations are given by the solid symbols with solid curves, and the coordination numbers for non-hydrogen-bonded configurations are given by the open symbols with broken curves. Also shown are the measured n̄ClD values from the Ti–Zr cell experiment (magenta solid downward triangles with vertical error bars) (Table III).

Close modal

Figure 29 maps the dependence of the O–O–D bond angle on the O–O distance. At ambient conditions, there is a peak near 2.7 Å and 8° and a second less-intense peak around 2.8 Å and 95°. Two additional features occur at longer O–O distances, one near 3.3 Å and 70° and the other at ≳3.6 Å and 52°. With increasing temperature, the peak at ∼2.7 Å broadens but remains distinct as the other features broaden and coalesce. With increasing pressure, the peak at ∼2.7 Å remains as a region of large intensity develops around 3.1 Å and 65°.

FIG. 29.

Map showing the simulated dependence of the O–O–D bond angle (Fig. 24) on the O–O distance. (a) T = 298 K, p = 0 kbar; (b) T = 323 K, p = 0 kbar; (c) T = 373 K, p = 0 kbar; (d) T = 423 K, p = 0 kbar; (e) T = 423 K, p = 10 kbar; and (f) T = 423 K, p = 33.8 kbar. The multiplicities are normalized to give a maximum value of unity.

FIG. 29.

Map showing the simulated dependence of the O–O–D bond angle (Fig. 24) on the O–O distance. (a) T = 298 K, p = 0 kbar; (b) T = 323 K, p = 0 kbar; (c) T = 373 K, p = 0 kbar; (d) T = 423 K, p = 0 kbar; (e) T = 423 K, p = 10 kbar; and (f) T = 423 K, p = 33.8 kbar. The multiplicities are normalized to give a maximum value of unity.

Close modal

The features at short O–O distances of 2.7–2.8 Å with small and large O–O–D bond angles originate from configurations in which a water molecule is either donating or accepting a hydrogen bond, respectively (Fig. 24). The feature at the O–O–D bond angle of 52° corresponds to longer O–O distances, at which the water molecules are not deemed to be hydrogen bonded (see below). From this breakdown, it follows that the peak at small angles in the O–O–D bond-angle distribution (Fig. 23) corresponds to short O–O distances and therefore originates from the formation of O⋯D–O hydrogen bonds in acceptor conformations (Fig. 24). The peak at 52° is at half the intra-molecular D–O–D bond angle and originates from non-hydrogen-bonded configurations. Its sharpness stems from the employment of rigid water molecules in the molecular dynamics simulations, as discussed in the supplementary material (Figs. S10 and S11).

Figure 28 shows the contributions to the coordination number n̄OO from hydrogen-bonded vs non-hydrogen-bonded configurations and compares the results to those obtained for pure D2O. Water molecules are deemed to be hydrogen bonded if the O–O–D bond angle <30 and the O–O distance <3.5 Å.68 At ambient conditions, the number of hydrogen-bonded water molecules in the solution is less than the value of ≃3.7 found for pure water. For both systems, the large increase in n̄OO with increasing pressure originates primarily from non-hydrogen-bonded water–water conformations, i.e., from water molecules at longer distances that are pushed into the first shell of hydrogen-bonded nearest-neighbors and thereby contribute toward the first peak in gOO(r). The sum of the hydrogen-bonded and non-hydrogen-bonded O–O coordination numbers gives an overall value n̄OO 12 for the solution vs n̄OO 14.5 for liquid D2O at the highest pressure [Fig. 28(b)]. The latter compares to the value n̄OO 12.6(4) measured for liquid H2O at 33 kbar and 500 K using x-ray diffraction.21 At room temperature, the measured O–O coordination number for liquid D2O is systematically larger than for H2O in the pressure range of up to ∼12 kbar.21 

As the packing of the anionic chlorine and oxygen species increases with pressure and the Cl–O and O–O coordination numbers grow (Fig. 28), the gClO(r) (Fig. 9) and gOO(r) (Fig. 22) functions show the development of pronounced ordering that extends to large r-values. In comparison, there are only modest changes to gNaO(r) (Fig. 13) and the Na–O coordination number (Fig. 10). The pressure-induced anionic ordering manifests itself in the dynamics by decreasing self-diffusion coefficients [Fig. 26(b)] but does not have a marked effect on the residence times [Fig. 25(b)].

The Cl–O ordering contributes toward the development with increasing pressure of a sharp first-peak at Q1 ∼ 2.45 Å−1 in the measured and simulated ΔFCl(Q) functions (Fig. 5): The Cl–O correlations receive a large weighting (Table II) and the wavelength of the oscillations in gClO(r) is given by 2π/Q1. There is a concomitant contribution to Q1 from the Cl–D correlations as indicated by the large-r oscillations in gClD(r) (Fig. 9), where this function describes the correlations between the chloride ions and the deuterium atoms of the water molecules and also receives a large weighting in the equation for ΔFCl(Q) (Table II). The sharpening of the first peaks in SClO(Q) and SClD(Q) with increasing pressure, which originates from the anion ordering, is shown in Fig. S12.

The in–out cage-correlation function70 

Cinout(t)=Θ1niout0,tΘ1niin0,t
(14)

depends on both the in and out motion of the water molecules. In the diffusive regime, Cin–out(t) ∝ exp(−t/τin–out), which defines the lifetime τin–out. For a concerted (interchange) exchange mechanism, one water molecule leaves a cage as another enters such that a distinct intermediate state is not formed. In this case, Cin–out(t) = Cout(t) and the ratio τout/τin–out = 1. For the opposite case of non-concerted exchange, the motion between incoming and outgoing water molecules is uncorrelated such that Cinout(t)=Cin(t)Cout(t)=Cout(t)2exp(2t/τout) and τout/τin–out = 2.

Figure 30 compares the cage-correlation functions for the exchange of water at the Na+ ion for several state points. The temperature and pressure dependence of the ratio τNaout/τNainout is shown in Fig. 31. At ambient conditions, a concerted water exchange mechanism is indicated by the similarity between the CNaout(t) and CNainout(t) functions, which give τNaout/τNainout 1.2. Under isobaric conditions, the in and out motion becomes non-concerted as the temperature is increased. Under isothermal conditions, the ratio τNaout/τNainout shows a broad minimum with increasing pressure. At the most extreme conditions, uncorrelated motion of the incoming and outgoing water molecules is indicated by the similarity between CNainout(t) and [CNaout(t)]2, which gives τNaout/τNainout 2.

FIG. 30.

Comparison of the simulated cage-correlation functions CNainout(t), CNaout(t), and CNaout(t)2 describing the exchange of water molecules at the Na+ ion for several state points.

FIG. 30.

Comparison of the simulated cage-correlation functions CNainout(t), CNaout(t), and CNaout(t)2 describing the exchange of water molecules at the Na+ ion for several state points.

Close modal
FIG. 31.

Dependence of the ratio τNaout/τNainout for water exchange at the Na+ ion obtained from the simulations on the temperature at 0 kbar vs pressure at 423 K.

FIG. 31.

Dependence of the ratio τNaout/τNainout for water exchange at the Na+ ion obtained from the simulations on the temperature at 0 kbar vs pressure at 423 K.

Close modal

Neutron diffraction is combined with molecular dynamics simulations to provide a comprehensive account of the structural and dynamical response of a concentrated solution of NaCl in water to extremes of temperature and pressure. The results show that the method of neutron diffraction with isotope substitution can now be used to give site-specific information on the ion solvation in aqueous solutions under conditions that extend into the gigapascal pressure regime. The simulations give a good account of the measured equation of state at low pressures and the empirical equation of state of Mantegazzi et al.5 at high pressures (Fig. 1). They also give a good account of the neutron diffraction results, including the pressure dependence of ΔFCl(Q) (Fig. 5). The simulated Cl–O coordination numbers for hydrogen-bonded configurations are in accord with the measured Cl–D coordination numbers (Fig. 28), although the simulated Cl–D coordination numbers for hydrogen-bonded plus non-hydrogen-bonded configurations are larger (Fig. 10).

It is found that the anionic species (chlorine and oxygen) do not pack uniformly with increasing density. Instead, pressure-induced anion ordering is observed, which manifests itself in the measured and simulated ΔFCl(Q) functions by the appearance of a sharp first-peak at ∼2.45 Å−1. This ordering originates from oscillations in gClO(r) and gOO(r) that extend to large distances in r-space and is accompanied by a large growth in the nearest-neighbor Cl–O and O–O coordination numbers, where the latter approaches n̄OO 12 at the highest pressure (Fig. 28). This enhancement does not emanate, however, from a large change to the number of Cl⋯D–O or O⋯D–O hydrogen bonds. Instead, there is a growth in the number of non-hydrogen-bonded configuration as the next nearest-neighbors push into the first coordination shell of the chloride ion or oxygen. The roles of hydrogen-bonded vs non-hydrogen-bonded configurations are elucidated with the aid of bond distance resolved bond-angle maps (Figs. 27 and 29).

Ion association manifests itself in terms of contact ion pairs, most of which are charge-neutral NaCl0 units and solvent-separated Na+–Na+ or Cl–Cl pairs. On the basis of dielectric continuum theory, the increase in ion pairing caused by the reduction in the dielectric constant of the pure solvent ɛ with increasing temperature is expected to be reversed by the enhancement in ɛ with increasing pressure: ɛ undergoes a cyclic variation under the conditions employed in the present work (Table I). Indeed, the fraction of solvent-separated Na+–Na+ pairs that are bridged by two water molecules decreases with increasing temperature (Fig. 16) before increasing by a similar amount with increasing pressure (Fig. 18). This observation appears, however, to be fortuitous, as befits a concentrated as opposed to dilute solution. For instance, the fraction of NaCl0 units increases with temperature but is relatively insensitive to increasing pressure (Fig. 15). The dependence of ion pairing on the force field chosen for the ions76 has not been investigated for the state conditions covered in the present work.

The ion–water binding times and self-diffusion coefficients were explored in the molecular dynamics simulations. The temperature dependence of this dynamics shows Arrhenius behavior. The pressure dependence was interpreted by assuming a transition state with an activation volume. Unlike the diffusion coefficients, the residence time of water molecules in the solvation cage of the Na+ and Cl ions shows only a weak pressure dependence. At ambient conditions, there is a concerted mechanism for the exchange of water at the Na+ ion that is disrupted as the state conditions are changed. It would be interesting to explore the impact on the dynamics of using a force field in which partial ion charges are combined with the TIP4P/2005 model for water.77,78

The observed evolution in the structural and dynamical properties of aqueous NaCl solutions under extreme conditions is important for understanding their role in phenomena that include the mechanisms of mass transport in geochemical processes (Sec. I) and the sequestration of CO2 in deep saline aquifers where the temperature, pressure, and fluid composition affect the dissolution and eventual fate of the gas.79 

See the supplementary material for the following: Fig. S1 shows the predictions of the molecular dynamics simulations for the partial pair-correlation functions that contribute toward ΔGCl(r). The n, the temperature and pressure dependence of the running coordination numbers (i) n̄ClD and n̄ClO (Fig. S2), (ii) n̄NaO and n̄NaD/2 (Fig. S3), and (iii) n̄NaCl (Fig. S4) are shown. Figure S5 gives the Na–Cl–Na and Cl–Na–Cl bond-angle distributions, and Fig. S6 gives the Na–O–Na bond-angle distributions. Figure S7 shows the temperature and pressure dependence of the running coordination number n̄OO. Figure S8 illustrates the temperature dependence of the cage-correlation functions CNaout(t) and CClout(t), and Fig. S9 illustrates the pressure dependence of these functions. Figure S10 gives a schematic showing the shorter and longer O–D inter-molecular distances for a configuration of two water molecules. Figure S11 shows a breakdown of the O–O–D bond-angle distribution into its contributions from hydrogen-bonded vs non-hydrogen-bonded configurations of the water molecules. Finally, Fig. S12 gives the contributions toward ΔFCl(Q) from the weighted Cl–O and Cl–D partial structure factors.

We thank Phil Mason for his advice on the sample preparation; Keiron Pizzey, Alain Bertoni, Claude Payre, and Jean-Luc Laborier for their help with the diffraction experiments; and Adrian Barnes and George Neilson for the provision of the Ti–Zr cell and chlorine isotopes. We also thank Yoshiki Ishii for his contributions toward molecular dynamics simulations using a different force field. The Bath group received support from the EPSRC via Grant No. EP/J009741/1. A.P. acknowledges funding and support from the Institut Laue Langevin (ILL) and the University of Bath (Collaboration Agreement No. ILL-1353.1). A.Z. was supported by a Royal Society-EPSRC Dorothy Hodgkin Research Fellowship. M.S. acknowledges HPC resources granted by the HPCaVe Centre at Sorbonne Université.

The authors have no conflicts to disclose.

R.F.R. and A.Z. prepared the samples; A.P., P.S.S., B.A., and H.E.F. performed the Ti–Zr cell experiment; and R.F.R., P.S.S., S.K., and H.E.F. performed the PE press experiment. A.P., R.F.R., and A.Z. analyzed the diffraction data. M.S. designed and performed the molecular dynamics simulations. M.S., A.Z., and P.S.S. analyzed the results. P.S.S. conceived the project and wrote the paper with input from all the co-authors.

The data that support the findings of this study are openly available in the University of Bath Research Data Archive at https://doi.org/10.15125/BATH-01054.80 The measured neutron diffraction datasets for the PE press and Ti–Zr cell experiments are available from Refs. 81 and 82, respectively.

2.
E. D.
Sloan
,
C. A.
Koh
, and
A. K.
Sum
,
Energies
3
,
1991
(
2010
).
3.
A.
Koschinsky
,
D.
Garbe-Schönberg
,
S.
Sander
,
K.
Schmidt
,
H.-H.
Gennerich
, and
H.
Strauss
,
Geology
36
,
615
(
2008
).
4.
C. E.
Manning
,
Earth Planet. Sci. Lett.
223
,
1
(
2004
).
5.
D.
Mantegazzi
,
C.
Sanchez-Valle
, and
T.
Driesner
,
Geochim. Cosmochim. Acta
121
,
263
(
2013
).
6.
E.
Roedder
, in , Physics and Chemistry of the Earth Vol. 13–14, edited by
D. T.
Rickard
and
F. E.
Wickman
(
Pergamon
,
Oxford
,
1981
), pp.
9
39
.
7.
R. J.
Bodnar
,
Geochim. Cosmochim. Acta
58
,
1053
(
1994
).
8.
R. J.
Bodnar
, in
Fluid Inclusions: Analysis and Interpretation
, Short Course Series Vol. 32, edited by
I.
Samson
,
A.
Anderson
, and
D.
Marshall
(
Mineralogical Association of Canada
,
Ottawa, Canada
,
2003
), Chap. 4, pp.
81
99
.
9.
S. K.
Trumbo
,
M. E.
Brown
, and
K. P.
Hand
,
Sci. Adv.
5
,
eaaw7123
(
2019
).
10.
R. C.
Newton
and
C. E.
Manning
,
Geofluids
10
,
58
(
2010
).
11.
E. H.
Oelkers
and
H. C.
Helgeson
,
Science
261
,
888
(
1993
).
12.
T. M.
Seward
, in , Physics and Chemistry of the Earth Vol. 13–14, edited by
D. T.
Rickard
and
F. E.
Wickman
(
Pergamon
,
Oxford
,
1981
), pp.
113
132
.
13.
D. A.
Sverjensky
,
B.
Harrison
, and
D.
Azzolini
,
Geochim. Cosmochim. Acta
129
,
125
(
2014
).
14.
H. C.
Helgeson
,
D. H.
Kirkham
, and
G. C.
Flowers
,
Am. J. Sci.
281
,
1249
(
1981
).
15.
G. D.
Miron
,
A. M. M.
Leal
, and
A.
Yapparova
,
Geofluids
2019
,
5750390
.
18.
A.
Zeidler
,
P. S.
Salmon
, and
L. B.
Skinner
,
Proc. Natl. Acad. Sci. U. S. A.
111
,
10045
(
2014
).
19.
A.
Zeidler
, “
Ordering in amorphous binary systems
,” Ph.D. thesis,
University of Bath
,
UK
,
2009
.
20.
T.
Strässle
,
A. M.
Saitta
,
Y.
Le Godec
,
G.
Hamel
,
S.
Klotz
,
J. S.
Loveday
, and
R. J.
Nelmes
,
Phys. Rev. Lett.
96
,
067801
(
2006
).
21.
G.
Weck
,
J.
Eggert
,
P.
Loubeyre
,
N.
Desbiens
,
E.
Bourasseau
,
J.-B.
Maillet
,
M.
Mezouar
, and
M.
Hanfland
,
Phys. Rev. B
80
,
180202(R)
(
2009
).
22.
Y.
Katayama
,
T.
Hattori
,
H.
Saitoh
,
T.
Ikeda
,
K.
Aoki
,
H.
Fukui
, and
K.
Funakoshi
,
Phys. Rev. B
81
,
014109
(
2010
).
23.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
24.
T.
Yagasaki
,
M.
Matsumoto
, and
H.
Tanaka
,
J. Chem. Theory Comput.
16
,
2460
(
2020
).
25.
J. E.
Enderby
,
S.
Cummings
,
G. J.
Herdman
,
G. W.
Neilson
,
P. S.
Salmon
, and
N.
Skipper
,
J. Phys. Chem.
91
,
5851
(
1987
).
26.
L. X.
Dang
,
G. K.
Schenter
,
V.-A.
Glezakou
, and
J. L.
Fulton
,
J. Phys. Chem. B
110
,
23644
(
2006
).
27.
M.
Galib
,
M. D.
Baer
,
L. B.
Skinner
,
C. J.
Mundy
,
T.
Huthwelker
,
G. K.
Schenter
,
C. J.
Benmore
,
N.
Govind
, and
J. L.
Fulton
,
J. Chem. Phys.
146
,
084504
(
2017
).
28.
T.
Yamaguchi
,
N.
Fukuyama
,
K.
Yoshida
, and
Y.
Katayama
,
J. Phys. Chem. Lett.
12
,
250
(
2021
).
29.
H. E.
Fischer
,
A. C.
Barnes
, and
P. S.
Salmon
,
Rep. Prog. Phys.
69
,
233
(
2006
).
30.
P. S.
Salmon
and
A.
Zeidler
,
J. Phys.: Condens. Matter
27
,
133201
(
2015
).
31.
D. M.
North
,
J. E.
Enderby
, and
P. A.
Egelstaff
,
J. Phys. C: Solid State Phys.
1
,
784
(
1968
).
32.
H. H.
Paalman
and
C. J.
Pings
,
J. Appl. Phys.
33
,
2635
(
1962
).
33.
A. K.
Soper
and
P. A.
Egelstaff
,
Nucl. Instrum. Methods
178
,
415
(
1980
).
34.
T. E.
Faber
and
J. M.
Ziman
,
Philos. Mag.
11
,
153
(
1965
).
35.
36.
37.
A. K.
Soper
,
G. W.
Neilson
,
J. E.
Enderby
, and
R. A.
Howe
,
J. Phys. C: Solid State Phys.
10
,
1793
(
1977
).
38.
E.
Lorch
,
J. Phys. C: Solid State Phys.
2
,
229
(
1969
).
39.
P. S.
Salmon
,
J. Phys.: Condens. Matter
18
,
11443
(
2006
).
40.
A.
Marin-Laflèche
,
M.
Haefele
,
L.
Scalfi
,
A.
Coretti
,
T.
Dufils
,
G.
Jeanmairet
,
S. K.
Reed
,
A.
Serva
,
R.
Berthin
,
C.
Bacon
,
S.
Bonella
,
B.
Rotenberg
,
P. A.
Madden
, and
M.
Salanne
,
J. Open Source Software
5
,
2373
(
2020
).
41.
J. L. F.
Abascal
,
E.
Sanz
, and
C.
Vega
,
Phys. Chem. Chem. Phys.
11
,
556
(
2009
).
42.
G. J.
Martyna
,
D. J.
Tobias
, and
M. L.
Klein
,
J. Chem. Phys.
101
,
4177
(
1994
).
43.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
,
J. Chem. Phys.
97
,
2635
(
1992
).
44.
H. C.
Andersen
,
J. Comput. Phys.
52
,
24
(
1983
).
45.
P. S.
Salmon
,
G. S.
Moody
,
Y.
Ishii
,
K. J.
Pizzey
,
A.
Polidori
,
M.
Salanne
,
A.
Zeidler
,
M.
Buscemi
,
H. E.
Fischer
,
C. L.
Bull
,
S.
Klotz
,
R.
Weber
,
C. J.
Benmore
, and
S. G.
MacLeod
,
J. Non-Cryst. Solids: X
3
,
100024
(
2019
).
46.
A.
Zeidler
,
P. S.
Salmon
,
R. A.
Martin
,
T.
Usuki
,
P. E.
Mason
,
G. J.
Cuello
,
S.
Kohara
, and
H. E.
Fischer
,
Phys. Rev. B
82
,
104208
(
2010
).
47.
H. E.
Fischer
,
G. J.
Cuello
,
P.
Palleau
,
D.
Feltin
,
A. C.
Barnes
,
Y. S.
Badyal
, and
J. M.
Simonson
,
Appl. Phys. A
74
,
S160
(
2002
).
48.
A.
Zeidler
,
P. S.
Salmon
,
H. E.
Fischer
,
J. C.
Neuefeind
,
J.
Mike Simonson
, and
T. E.
Markland
,
J. Phys.: Condens. Matter
24
,
284126
(
2012
).
49.
I.
Howell
and
G. W.
Neilson
,
J. Chem. Phys.
104
,
2036
(
1996
).
50.
B.
Annighöfer
,
A.
Polidori
,
A.
Zeidler
,
H. E.
Fischer
, and
P. S.
Salmon
,
High Pressure Res.
37
,
529
(
2017
).
51.
A.
Polidori
, “
Structure of disordered materials: From geological fluids to network glasses
,” Ph.D. thesis,
University of Bath
,
UK
,
2017
.
52.
A.
Zeidler
,
M.
Guthrie
, and
P. S.
Salmon
,
High Pressure Res.
35
,
239
(
2015
).
53.
P. S.
Salmon
,
A.
Zeidler
, and
H. E.
Fischer
,
J. Appl. Cryst.
49
,
2249
(
2016
).
54.
S.
Klotz
,
T.
Strässle
,
A. L.
Cornelius
,
J.
Philippe
, and
V.
Pomjakushin
,
J. Phys. D: Appl. Phys.
44
,
055406
(
2011
).
55.
W. G.
Marshall
and
D. J.
Francis
,
J. Appl. Cryst.
35
,
122
(
2002
).
56.
R. F.
Rowlands
, “
The role of structural disorder in crystalline, glassy and liquid materials
,” Ph.D. thesis,
University of Bath
,
UK
,
2015
.
57.
C. W. F. T.
Pistorius
,
E.
Rapoport
, and
J. B.
Clark
,
J. Chem. Phys.
48
,
5509
(
1968
).
58.
R. W.
Potter
 II
and
D. L.
Brown
, Geological Survey Bulletin 1421-C,
United States Department of the Interior
,
Washington
,
1977
; available at https://pubs.usgs.gov/bul/1421c/report.pdf.
59.
M.
Uematsu
and
E. U.
Frank
,
J. Phys. Chem. Ref. Data
9
,
1291
(
1980
).
61.
A.
Zeidler
,
P. S.
Salmon
,
H. E.
Fischer
,
J. C.
Neuefeind
,
J. M.
Simonson
,
H.
Lemmel
,
H.
Rauch
, and
T. E.
Markland
,
Phys. Rev. Lett.
107
,
145501
(
2011
).
62.
G.
Walford
,
J. H.
Clarke
, and
J. C.
Dore
,
Mol. Phys.
33
,
25
(
1977
).
63.
P. S.
Salmon
and
P. B.
Lond
,
J. Phys.: Condens. Matter
4
,
5249
(
1992
).
64.
S. E.
Okan
and
P. S.
Salmon
,
J. Phys.: Condens. Matter
6
,
3839
(
1994
).
65.
S.
Cummings
,
J. E.
Enderby
,
G. W.
Neilson
,
J. R.
Newsome
,
R. A.
Howe
,
W. S.
Howells
, and
A. K.
Soper
,
Nature
287
,
714
(
1980
).
66.
D. H.
Powell
,
A. C.
Barnes
,
J. E.
Enderby
,
G. W.
Neilson
, and
P. S.
Salmon
,
Faraday Discuss. Chem. Soc.
85
,
137
(
1988
).
67.
C.
Zhang
,
F.
Giberti
,
E.
Sevgen
,
J. J.
de Pablo
,
F.
Gygi
, and
G.
Galli
,
Nat. Commun.
11
,
3037
(
2020
).
68.
A.
Luzar
,
J. Chem. Phys.
113
,
10663
(
2000
).
69.
R. D.
Shannon
,
Acta Crystallogr., Sect. A: Found. Adv.
32
,
751
(
1976
).
70.
E.
Rabani
,
J. D.
Gezelter
, and
B. J.
Berne
,
J. Chem. Phys.
107
,
6867
(
1997
).
71.
O.
Pauvert
,
M.
Salanne
,
D.
Zanghi
,
C.
Simon
,
S.
Reguer
,
D.
Thiaudière
,
Y.
Okamoto
,
H.
Matsuura
, and
C.
Bessada
,
J. Phys. Chem. B
115
,
9160
(
2011
).
72.
P. S.
Salmon
,
J. Phys. C: Solid State Phys.
20
,
1573
(
1987
).
73.
P. S.
Salmon
,
W. S.
Howells
, and
R.
Mills
,
J. Phys. C: Solid State Phys.
20
,
5727
(
1987
).
74.
L.
Helm
and
A. E.
Merbach
,
Chem. Rev.
105
,
1923
(
2005
).
75.
C. T.
Chudley
and
R. J.
Elliott
,
Proc. Phys. Soc.
77
,
353
(
1961
).
76.
L. A.
Patel
,
T. J.
Yoon
,
R. P.
Currier
, and
K. A.
Maerzke
,
J. Chem. Phys.
154
,
064503
(
2021
).
77.
A. L.
Benavides
,
M. A.
Portillo
,
V. C.
Chamorro
,
J. R.
Espinosa
,
J. L. F.
Abascal
, and
C.
Vega
,
J. Chem. Phys.
147
,
104501
(
2017
).
78.
I. M.
Zeron
,
J. L. F.
Abascal
, and
C.
Vega
,
J. Chem. Phys.
151
,
134504
(
2019
).
79.
G. P. D.
De Silva
,
P. G.
Ranjith
, and
M. S. A.
Perera
,
Fuel
155
,
128
(
2015
).
80.
P. S.
Salmon
and
A.
Zeidler
(
2021
), “
Data sets for structure and dynamics of aqueous NaCl solutions at high temperatures and pressures
,” University of Bath Research Data Archive, .
81.
P. S.
Salmon
,
H. E.
Fischer
,
S.
Klotz
,
K. J.
Pizzey
,
R. F.
Rowlands
, and
A.
Zeidler
(
2013
), “
Structure of geological fluids
,” Institut Laue-Langevin, Grenoble, France, https://doi.ill.fr/10.5291/ILL-DATA.6-02-502.
82.
P. S.
Salmon
,
B.
Annighöfer
,
H. E.
Fischer
,
S.
Klotz
,
A.
Polidori
,
R. F.
Rowlands
, and
A.
Zeidler
(
2015
), “
Structural transformation in water at high pressure
,” Institut Laue-Langevin, Grenoble, France, https://doi.ill.fr/10.5291/ILL-DATA.6-02-547.

Supplementary Material