Effect of Greenhouse Gases Dissolved in Seawater

Download PDF
2 downloads 0 Views 7MB Size Report
Comment
Dec 30, 2015 - is executed in the number of constituent molecules, the temperature, and the pressure (NTP) constant condition ..... 1.5565 ˆ 102. 8.3931.
Article

Effect of Greenhouse Gases Dissolved in Seawater Shigeki Matsunaga Received: 30 September 2015; Accepted: 14 December 2015; Published: 30 December 2015 Academic Editor: Malcolm D´Souza National Institute of Technology, Nagaoka College, Nishikatakai 888, Nagaoka 940-8532, Japan; [email protected]; Tel./Fax: +81-258-34-9252

Abstract: A molecular dynamics simulation has been performed on the greenhouse gases carbon dioxide and methane dissolved in a sodium chloride aqueous solution, as a simple model of seawater. A carbon dioxide molecule is also treated as a hydrogen carbonate ion. The structure, coordination number, diffusion coefficient, shear viscosity, specific heat, and thermal conductivity of the solutions have been discussed. The anomalous behaviors of these properties, especially the negative pressure dependence of thermal conductivity, have been observed in the higher-pressure region. Keywords: molecular dynamics; carbon dioxide; methane; hydrogen carbonate; seawater; structure; transport properties; thermal properties

1. Introduction The reason for the global warming of the earth’s climate seems undoubtedly attributable to the increasing of greenhouse gases, especially carbon dioxide (CO2 ) in the atmosphere. The consumption of fossil fuels since the industrial revolution has released a large amount of CO2 into the air, the concentration of which is expected to exceed 400 ppm in 2015 [1]. Therefore, technologies for reducing CO2 emissions are a matter of the utmost importance, and have widely been investigated, including CO2 absorbing processes into liquids [2]. In addition to the soluble property of CO2 in water, the increasing concentration of CO2 in the atmosphere and the warmer climate may accelerate the dissolution of CO2 into seawater, which can simply be modeled as a sodium chloride (NaCl) aqueous solution [3]. However, few molecular dynamics (MD) studies on CO2 and NaCl aqueous solutions are available in the literature [4]. Although many experimental studies on CO2 dissolved into a NaCl aqueous solution have been published, the solubility of CO2 and/or the stability of the solution are their main purpose [4–6]. To the best of our knowledge, the thermal conductivity of the CO2 absorbed NaCl aqueous solution by MD simulation has not been reported up to the present. Considering the situation, we wish to show the results of a MD simulation of CO2 absorbed into a NaCl aqueous solution. The CO2 concentration is postulated to be that of saturation. The pressure of the solution in the MD calculation corresponds to the depth of the ocean [7]. The dissolved CO2 molecules, however, react with water molecules to create bicarbonate (HCO3 ´ ) and carbonate (CO3 2´ ) ions in the seawater. The concentrations of these ions, or which ions are more abundant, depend on the acidity of the seawater. The concentration of CO3 2´ ions increases in alkaline regions, whereas CO2 molecules are dominant in acidic regions. Because the average pH of seawater is ordinarily between 7.9 and 9.0, the more abundant ion is HCO3 ´ , i.e., the following reaction is postulated [3]: (1) CO2 ` H2 O Ñ HCO3 ´ ` H+ which has the function of keeping the pH of seawater constant to a certain extent, and guarantees the charge neutrality of the whole system. As the continuous investigation, we will show the results of MD on seawater saturated with HCO3 ´ ions as a more realistic model [8]. Int. J. Mol. Sci. 2016, 17, 45; doi:10.3390/ijms17010045

www.mdpi.com/journal/ijms

Int. J. Mol. Sci. 2016, 17, 45

1 of 17

which has the function of keeping the pH of seawater constant to a certain extent, and guarantees the charge neutrality of the whole system. As the continuous investigation, we will show the results of MD on seawater saturated with HCO3− ions as a more realistic model [8]. Int. J. Mol. Sci. 2016, 17, 45 2 of 18 Although methane is another the greenhouse gas, the methane and water system has attracted more attention as an energy resource, i.e., the methane hydrate, which is formed in the low temperature and high pressure region at the bottom themethane sea [9]. The of has the attracted methane Although methane is another the greenhouse gas,ofthe and dissolution water system into seawater the resource, atmosphere alsomethane seems to affect the thermal environment the earth. more attention and as aninto energy i.e., the hydrate, which is formed in the low of temperature Many MD studiesregion have been onthe thesea methane water system; the nucleation or and high pressure at thereported bottom of [9]. Theand dissolution of the however, methane into seawater and the structural properties in the mixture pure water are their of main [10,11]. this study, into the atmosphere also seems to affect with the thermal environment the purpose earth. Many MDIn studies have we wish to show themethane results of a MD simulation of the methane and NaCl aqueous solution system, been reported on the and water system; however, the nucleation or the structural properties in as amixture model with of methane dissolved into seawater. change of structure, the coordination number, the pure water are their main purposeThe [10,11]. In this study, we wish to show the results of thermalofconductivity be discussed various corresponding to aand MDthe simulation the methane will and NaCl aqueous under solutionthe system, as apressures, model of methane dissolved different sea depths [12]. We believe the results of these investigations will be a part of the into seawater. The change of structure, the coordination number, and the thermal conductivity will fundamental data ofthe thevarious ocean environment. be discussed under pressures, corresponding to different sea depths [12]. We believe the results of these investigations will be a part of the fundamental data of the ocean environment. 2. Results and Discussion 2. Results and Discussion In this section, first the results of MD simulation are described for the CO2 and NaCl aqueous +, Cl −, H +, and In this section, first MD simulation describedNa for the CO NaCl aqueous solution system; then forthe the results systemof that includes waterare molecules, HCO 3− ions; and 2 and + , Cl´ , H+ , and HCO ´ ions; solution system; then for the system that includes water molecules, Na finally for the methane and NaCl aqueous solution system. 3 and finally for the methane and NaCl aqueous solution system. 2.1. The CO2 and NaCl Aqueous Solution 2.1. The CO2 and NaCl Aqueous Solution In the study for CO2 and NaCl aqueous solution, the equilibrium MD (EQMD) is used to obtain In the study for CO2 and number, NaCl aqueous solution, equilibriumand MDthe (EQMD) usedThe to the structure, the coordination the mean squarethe displacement, specificisheat. obtain the structure, the coordination number, the in mean displacement, andwith the specific fundamental procedure of EQMD is same as that oursquare previous works dealing molten heat. salts The fundamental procedure of EQMD is same as that in our previous works dealing with molten [13,14]. As mentioned in the previous section, the thermal conductivity is evaluated using the nonsalts [13,14]. As the previous section, thermal is evaluated using the equilibrium MDmentioned (NEMD) in [15,16]. The rigid bodythe models areconductivity used for liquid water using the non-equilibrium MD (NEMD) [15,16]. The rigid body models are used for liquid water using the transferable potential functions of 4-site model (TIP4P) [17] and for CO2 molecules. The MD transferable potential functions of 4-site model (TIP4P) [17] and for CO molecules. The MD calculation calculation is executed in the number of constituent molecules, the2temperature, and the pressure is executed in the number of constituent molecules, the temperature, and the pressure (NTP) (NTP) constant condition [18–20]. The applied pressure varies from 0.5 to 100 MP, or 100constant × 10 P, P, which condition [18–20]. The applied from 0.5 100isMP, or 100 ˆfor 10650,000 is equivalent which is equivalent to an oceanpressure depth ofvaries 40 to 10,000 m.toMD carried out steps with 0.1 fs to an interval. ocean depth 40 to 10,000 of m.NaCl MD isincarried out is forthe 50,000 with fs time interval. The time The of concentration the water samesteps as that of0.1 seawater, or 1.1 mol % concentration of NaCl in the water is the same as that of seawater, or 1.1 mol % NaCl. The concentration NaCl. The concentration of CO2 is adjusted to that of a saturated solution. The used molecular of CO2 isfor adjusted to thatcalculation of a saturated solution. Thein used molecular numbers of forthe thesystem structure numbers the structure at 275 K are listed Table 1. The parameters are calculation at 275 K aresteps listed in Table 1. The parameters of the system are monitored monitored for 50,000 under 30 MP at 275 K to guarantee convergence (Figure for 1). 50,000 Aboutsteps 5000 under 30 MP at 275 K to guarantee convergence 1). About 5000 molecules (15,000 atoms) are molecules (15,000 atoms) are used for the thermal(Figure conductivity calculation so as to contribute to rapid used for the thermal conductivity calculation so as to contribute to rapid convergence. convergence.

Figure 1. 1. The 50,000 steps under 30 30 MPMP at 275 K. K. T: Figure Themonitor monitorparameters parametersofofthe thesystem systemduring during 50,000 steps under at 275 temperature, P: pressure, V: volume, U: internal energy. T: temperature, P: pressure, V: volume, U: internal energy.

Int. J. Mol. Sci. 2016, 17, 45

3 of 18

Table 1. The number CO2 and NaCl molecules used in the aqueous solution at 275 K. Int. J. Mol. Sci. 2016, 17, 45

2 of 17

Pressure

Na+

Water

Cl´

CO2

112 − Cl 112 112 112 112 112 112 112

48 CO2 575 48 750 575 812 750 812

Table 1. The number CO2 and NaCl molecules used in the aqueous solution at 275 K.

0.5 MP (283 K) 10,000 112 Pressure 10,000 Water 112 Na+ 30 MP 600.5 MP 10,000 112 MP (283 K) 10,000 112 100 MP 30 MP 10,000 10,000 112 112 60 MP 10,000 112 100 MP 10,000 112

The pair distribution functions, gij (r), and the integrated coordination number, nij (r), have been obtained from the MD results, and are shown in Figure 2a–c under 0.5, 30, 60 and 100 MP. The slight The pair distribution functions, gij(r), and the integrated coordination number, nij(r), have been ´ –H O and H O–H O can enhancement of the heights ofshown gij (r) of Na+ –H 2 –H 2 O, under 2 O, 2100 MP. The 2 slight 2 obtained from thefirst MDpeak results, and are in CO Figure 2a–c 0.5, 30,Cl 60 and + − be observed under higher pressures, which may correspond to the cage formation by H O molecules. enhancement of the first peak heights of gij(r) of CO2–H2O, Na –H2O, Cl –H2O and H2O–H22O can be To examine structural change by pressure more clearly, have calculated coordination observedthe under higher pressures, which may correspond to thewe cage formation by H2the O molecules. To examine the structural change by pressure more we obtained have calculated the coordination number, or the value of nij (r) at the first minimum of gclearly, coordination numbers are ij (r). The of nijof (r)fact, at the first minimum of gij(r). The coordination numbers are w (r) listednumber, in Tableor2.the Asvalue a matter the pressure dependence ofobtained the coordination numbers, nNaO listed in Table 2. As a matter of fact, the pressure dependence of the coordination numbers, nNaO w(r) and nClOw (r), are in the margin of error, although a slight decreasing tendency can be observed. and nClOw(r), are in the margin of error, although a slight decreasing tendency can be observed. The The decreasing coordination number for nCCO2 Ow (r) is obviously seen. These facts might be evidence decreasing coordination number for nCCO2Ow(r) is obviously seen. These facts might be evidence of of the cage formation around the CO2 molecules, because the cage structure yields the space around the cage formation around the CO2 molecules, because the cage structure yields the space around CO2 molecules, thereby decreasing the coordination numbers. CO2 molecules, thereby decreasing the coordination numbers.

(a)

(b)

(c)

Figure 2. (a) gCCO2Ow(r) and nCCO2Ow(r) under 0.5 (283 K), 30, 60, and 100 MP at 275 K; (b) gNaOw(r)

Figure 2. (a) gCCO2 Ow (r) and nCCO2 Ow (r) under 0.5 (283 K), 30, 60, and 100 MP at 275 K; (b) gNaOw (r) and nNaOw(r); and (c) gClOw(r) and nClOw(r) under 0.5 (283 K), 30, 60, and 100 MP at 275 K. and nNaOw (r); and (c) gClOw (r) and nClOw (r) under 0.5 (283 K), 30, 60, and 100 MP at 275 K.

Int. J. Mol. Sci. 2016, 17, 45

4 of 18

2.45The Int. J. Mol. Sci.Table 2016, 17,

coordination numbers for nCCO2 Ow (r), nNaOw (r), and nClOw (r).

Int. J. Mol. Sci. 2016,Table 17, 45 2.Pressure nCnumbers nClO (r) w(r). The coordination for nCnNaO CO2Ow(r), nNaOw(r), andwnClO w (r) CO2 Ow (r)

3 of 17 3 of 17

MP 20.1 nC CO2 Ow(r)for nC nNaO nClO w(r) Table 2. 0.5 ThePressure coordination numbers CO2Ow w(r) (r), nNaO w(r), and nClOw(r).

30 MP 19.3 0.5 MP 20.1 5.5 ˘ 0.1 6.4 ˘ 0.2 60 MP 19.1 Pressure nCCO2 Ow(r) nNaOw(r) nClOw(r) 30 MP 19.3 100 MP 18.7 5.5 ± 0.1 6.4 ± 0.2 0.5 MP 20.1 60 MP 19.1 30 MP 19.3 100 MP 18.7 5.5 ± 0.1 6.4 ± 0.2 60 MP The interference functions, Iij (q), which19.1 is also obtainable from the neutron diffraction experiment, 100 MPIij(q), which 18.7 is also obtainable from the neutron diffraction The interference functions,

can be calculated from the Fourier transformation of the pair distribution function gij (r) as expressed in experiment, canThe be calculated thefrom Fourier transformation the are pairshown distribution function (r)sharp Equation (10) [21]. evaluatedfrom Iij (q)s gij (r)s obtained byofMD in Figure 3a,b.gijA The interference functions, Iij(q), which is also obtainable from the neutron diffraction as expressed in Equation (10) [21]. The evaluated I ij(q)s from gij(r)s obtained by MD are shown in ´ 1 peak experiment, of IC-Ow (q) can at 2be [Åcalculated ] at 283from K, 0.5 in Figure 3a is extremely enhanced at 275 K, 100 MP in theMP Fourier transformation of the pair distribution function gij(r) Figure 3a,b. A sharp peak of IC-Ow(q) at 2 [Å−1] at 283 K, 0.5 MP in Figure 3a is extremely enhanced at Figure 3b. Although the peak heights of I (q) and I (q) are similar to that of I (q) at 283 K, C-Ow Cl-Ow as expressed in Equation (10) [21]. TheNa-Ow evaluated Iij(q)s from gij(r)s obtained by MD are shown in 275 K, 100 MP in Figure 3b. Although the peak heights of INa-Ow(q) and ICl-Ow(q) are similar to that of −1 0.5 MP in Figure theirpeak heights decrease widths K, 100 enhanced MP in Figure Figure 3a,b. 3a, A sharp of IC-Ow (q) at 2and [Å their ] at 283 K, 0.5 are MPbroadened in Figure 3aatis275 extremely at 3b. IC-Ow(q) at 283 K, 0.5 MP in Figure 3a, their heights decrease and their widths are broadened at 275 K, 275 K, 100 MP in Figure 3b.also Although the peaktoheights of INa-Ow(q) and ICl-Owof (q)water are similar to that as of they These100 structure changes may be attributed the structure formation molecules, MP in Figure 3b. These structure changes may also be attributed to the structure formation of IC-Ow(q) at as 283the K, concentration 0.5 MP in Figureof3a, their heights decrease and their widths are broadened at 275 K, are enhanced CO 2 increases. water molecules, as they are enhanced as the concentration of CO2 increases. 100 MP in Figure 3b. These structure changes may also be attributed to the structure formation of water molecules, as they are enhanced as the concentration of CO2 increases.

Figure 3. (a) Iij(r) at 283 K, 0.5 MP and (b) Iij(r) at 275 K, 100 MP.

Figure 3. (a) Iij (r) at 283 K, 0.5 MP and (b) Iij (r) at 275 K, 100 MP. Figure 3. (a) Iij(r) at 283 K, and (b) Iijproperties, (r) at 275 K, 100 To examine the pressure dependence of0.5 theMP transport theMP. diffusion coefficients of

To examine the pressure dependence of the transport properties, the diffusionofcoefficients of constituent molecules and ions are calculated, which are obtainable from the inclination the mean To examine the pressure dependence of the transport properties, the diffusion coefficients of square displacement (MSD) using Equation (11). The obtained MSD for 30 MP is shown in Figure 4. constituent molecules and ions are calculated, which are obtainable from the inclination of the mean constituent molecules and ions are calculated, which are obtainable from −the inclination of the mean Although small oscillations still remain in the diffusion coefficient for Cl theMP inclinations of in MSDs square displacement (MSD) using Equation (11). The obtained MSD for, 30 is shown Figure 4. square displacement (MSD) using Equation (11). The obtained MSD for 30 ´ MP is shown in Figure 4. seem to converge until 5 × 10 s (5 ps) and the slopes are kept for longer periods, to a certain extent. Although small oscillations still remain in the diffusion coefficient for Cl , the inclinations of MSDs Although small oscillations still remain in the diffusion coefficient for Cl−, the inclinations of MSDs The pressure dependence of´the diffusion coefficients, Di, is listed in Table 3. The decreasing tendency 12 seem seem to converge until 5 ˆ510× 10 s (5s ps) and the forlonger longer periods, a certain extent. to converge until (5 ps) and theslopes slopesare are kept kept for periods, to atocertain extent. of the diffusion coefficients also suggests the structure formation in the higher-pressure regions. The pressure dependence of the diffusion listedinin Table 3. The decreasing tendency The pressure dependence of the diffusioncoefficients, coefficients, Dii,, isislisted Table 3. The decreasing tendency of theofdiffusion coefficients also suggests thehigher-pressure higher-pressure regions. the diffusion coefficients also suggeststhe thestructure structure formation formation ininthe regions.

Figure 4. The MSDs of constituent molecules and ions under pressure 30 MP and temperature 275 K. Figure 4. The MSDs of constituent moleculesand and ions ions under 3030 MP and temperature 275 K. Figure 4. The MSDs of constituent molecules underpressure pressure MP and temperature 275 K.

Int. J. Mol. Sci. 2016, 17, 45

5 of 18

Int. J. Mol. Sci. 2016, 17, 45

4 of 17

Table 3. The pressure dependence of Di (ˆ10´5 cm2 /s) at 275 K. Table 3. The pressure dependence of Di (×10−5 cm2/s) at 275 K. Pressure 2 2 Pressure CO CO 0.5 0.5 MP MP 30 MP 30 MP 60 MP 60 MP 100 MP

100 MP

1.42 1.42 1.24 1.24 1.17 1.17 0.91

0.91

+

Na + Na 1.28 1.28 0.70 0.70 0.55 0.55 0.52 0.52

´

− Water ClCl Water 1.58 2.652.65 1.58 1.13 1.651.65 1.13 0.94 1.47 0.94 1.47 0.70 1.30 0.70 1.30

Next, to investigate investigate the the thermal thermal properties properties of of 1.1 1.1 mol mol % % Next, we we calculate calculate the the thermal thermal conductivity conductivity to NaCl aqueous solution with saturated CO . To the best of the author’s knowledge, this is the first NaCl aqueous solution with saturated CO22. To the best of the author’s knowledge, this is the first report on the thethermal thermalconductivity conductivityofofCO CO NaCl aqueous solution system obtained by MD. 2 and report on 2 and NaCl aqueous solution system obtained by MD. The The perturbation is applied to the system in the thermal equilibrium at time t = 0. According to the perturbation is applied to the system in the thermal equilibrium at time t = 0. According to the GreenGreen-Kubo (G-K) formula, the thermal conductivity λ is obtained using Equations (14)–(16). Kubo (G-K) formula, the thermal conductivity λ is obtained using Equations (14)–(16). The thermal conductivity of the aqueous solution of the molecule containing a few atoms can also The thermal conductivity of the aqueous solution of the molecule containing a few atoms can be derived by NEMD, postulating that the contribution from the atom in the same molecule is omitted also be derived by NEMD, postulating that the contribution from the atom in the same molecule is from thefrom perturbation current [22]. [22]. The The thermal conductivity obtained omitted the perturbation current thermal conductivity obtainedbybyNEMD NEMDunder under various various pressures data of of pure pressures is is shown shown in in Figure Figure 5a, 5a, alongside alongside the the experimental experimental data pure water water and and 11 mol mol % % NaCl NaCl aqueous solution [23,24]. The results of NEMD for the saturated concentration of CO in 1.1 mol % % aqueous solution [23,24]. The results of NEMD for the saturated concentration of CO22 in 1.1 mol NaCl solution deserves deserves special special mention: mention: the the thermal thermal conductivity conductivity decreases decreases above above 80 80 MP, MP, NaCl aqueous aqueous solution which forms a striking contrast with the positive pressure dependence of other thermal conductivity which forms a striking contrast with the positive pressure dependence of other thermal conductivity data is not not included. This anomaly anomaly of of thermal thermal conductivity data on on solutions, solutions, in in which which CO CO22 is included. This conductivity also also signifies signifies the structure change of the CO absorbed NaCl aqueous solution under high pressure. 2 the structure change of the CO2 absorbed NaCl aqueous solution under high pressure.

(a)

(b)

Figure 5. (a) The pressure dependence of the thermal conductivity; (b) The pressure dependence of Figure 5. (a) The pressure dependence of the thermal conductivity; (b) The pressure dependence of the specific specific heat heat at at the the constant constant pressure. pressure. In the horizontal horizontal axes axes show show pressure. pressure. The the In (a,b), (a,b), the The abbreviates abbreviates “aq.”, “sol.”, “exp.” stand for “aqueous”, “solution”, and “experiment”, respectively. The dotted lines lines “aq.”, “sol.”, “exp.” stand for “aqueous”, “solution”, and “experiment”, respectively. The dotted in (a) and the orange and blue lines in (b) are drawn to guide the reader’s eyes. in (a) and the orange and blue lines in (b) are drawn to guide the reader’s eyes.

Furthermore, the specific heat at constant pressure, Cp, obtained by MD using Equation (13) is Furthermore, the specific heat at constant pressure, C p , obtained by MD using Equation (13) is shown in Figure 5b. The significant increase of Cp under higher pressure can be seen, although the shown in Figure 5b. The significant increase of C p under higher pressure can be seen, although the experimental and MD data for pure water and seawater show a decreasing tendency of Cp against experimental and MD data for pure water and seawater show a decreasing tendency of C p against pressure. These results suggest the possibility of heat storage in the depths of the sea. pressure. These results suggest the possibility of heat storage in the depths of the sea.

Int. J. Mol. Sci. 2016, 17, 45

6 of 18

2.2. The System Including Water Molecule, Na+ , Cl´ , H+ , and HCO3 ´ Ions As stated in the previous subsection, the thermodynamic properties of CO2 and NaCl aqueous Int. J. Mol. Sci. 2016, 17, 45 5 of 17 solution show the anomalous features under high pressure. These facts prompt us to a further study. ´ in the neutral pH. −, H+, and HCO As mentioned in theIncluding previous section, the Na CO+,2Clmolecule is ionized 2.2. The System Water Molecule, 3− Ions to form HCO3 In this study, we wish to show the results of MD on seawater saturated with HCO3 ´ ions as a more As stated in the previous subsection, the thermodynamic properties of CO2 and NaCl aqueous ´ from 0.44 realistic model. MD is performed in 1.1 mol % NaCl aqueous solution with saturated HCO solution show the anomalous features under high pressure. These facts prompt us to a further3study. − in contains to 7.97As mol %, corresponding to 5–1200 atm Then, foristhe calculation, MD 3cell mentioned in the previous section, the[25]. CO2 molecule ionized to form HCO the neutral2500 pH. TIP4P, + ´ + ´ thisClstudy, we11wish to show the results 28 Na In, 28 , and to 219 H and HCO3 of. MD on seawater saturated with HCO3− ions as a more realistic MD is performed NaCl aqueous solution with saturated 3− from The pair model. distribution function, in gij1.1 (r),mol and%the integrated coordination number,HCO nij (r), obtained 0.44 to 7.97 mol %, corresponding to 5–1200 atm [25]. Then, for the calculation, MD cell contains 2500 from MD results are shown in Figures 6 and 7 under pressures from 5 to 1000 atm. As seen TIP4P, 28 Na+, 28 Cl−, and 11 to 219 H+ and HCO3−. in Figure 6a, gCOw (r) between C(HCO3 ´ ) and O(water) has a pronounced peak at around 4 Å. The pair distribution function, gij(r), and the integrated coordination number, nij(r), obtained The coordination number of water around a HCO3 ´pressures is estimated to be 17, which agrees well with from MD results are shown in Figures 6 and 7 under from 5 to 1000 atm. As seen in Figure those 6a, in the literature [26]. The sharp first peaks of g (r) are found around Å in Figure 6b, − NaOw gCOw(r) between C(HCO3 ) and O(water) has a pronounced peak at around 4 Å. The2.3 coordination − whichnumber shows the close distance between cations and water molecules. Figure 6a,b correspond to those of water around a HCO3 is estimated to be 17, which agrees well with those in the literature [26]. The sharp first of gNaOw(r) are found around 2.3 confirm Å in Figure which shows the closeof the of C(CO andpeaks Na–O(water) in Figure 2a,b. To the6b, pressure dependence 2 )–O(water) distancenumber, between cations Figure 6a,b of C(CO2)–O(water) coordination which and waswater seen molecules. in the Section 2.1, wecorrespond calculate to thethose coordination number to the and Na–O(water) Figure 2a,b. To confirm pressure dependence of the coordination number, first minimum of gij (r).inThe obtained results arethe listed in Table 4. The negative pressure dependence of which was seen in the Section 2.1, we calculate the coordination number to the first minimum of gij(r). the coordination number is also observed, which is also the collaborating evidence of the structure The obtained results are listed in Table 4. The negative pressure dependence of the coordination formation around HCO3 ´ ions. In Figure 7a,b, gCNa (r), C(HCO3 ´ )–Na+ , gCH (r), and C(HCO3 ´ )–H+ number is also observed, which is also the collaborating evidence of the structure formation around have pronounced two split7a,b, peaks 2.5 to3−4.0 Å,+, which mayC(HCO be attributed to the asymmetric form HCO3− ions. In Figure gCNafrom (r), C(HCO )–Na gCH(r), and 3−)–H+ have pronounced two ´. of HCO 3 split peaks from 2.5 to 4.0 Å, which may be attributed to the asymmetric form of HCO3−.

(a)

(b) Figure 6. (a) gCOw(r) and nCOw(r) under pressures of 5–1000 atm; and (b) gNaOw(r) and nNaOw(r) under

Figure 6. (a) gCOw (r) and nCOw (r) under pressures of 5–1000 atm; and (b) gNaOw (r) and nNaOw (r) under pressures of 5–1000 atm. pressures of 5–1000 atm.

Int. J. Mol. Sci. 2016, 17, 45

7 of 18

Int. J. Mol. Sci. 2016, 17, 45

6 of 17

Table 4. The coordination numbers for nCOw (r), nNaOw (r), and nClOw (r). Table 4. The coordination numbers for nCOw(r), nNaOw(r), and nClOw(r). Int. J. Mol. Sci. 2016, 17, 45

Pressure n (r) n (r) n (r) Pressure COw nCOw(r) nNaOwNaOw (r) nClOw(r) ClOw 19.1 5.31 Table 54.atm The coordination numbers for5.31 nCOw (r), nNaOw nClOw(r). 5 atm 19.1 6.55(r), and6.55 100 atm 100 atm 18.3 4.65 6.49 18.3 4.65 6.49 COw(r) nNaOw(r) 500 atmPressure n 17.4 4.57 nClOw(r) 6.00 500 atm 17.4 4.57 6.00 19.1 5.314.43 6.55 800 atm 5 atm 17.3 5.51 800 atm 17.3 4.43 5.51 1000 atm100 atm 16.5 5.21 18.3 4.654.18 6.49 1000 atm 16.5 4.18 5.21 500 atm 17.4 4.57 6.00 800 atm 17.3 4.43 5.51 1000 atm 16.5 4.18 5.21

6 of 17

Figure 7. (a) gCNa(r) under pressures of 5–1000 atm; and (b) gCH(r) under pressures of 5–1000 atm.

Figure 7. (a) gCNa (r) under pressures of 5–1000 atm; and (b) gCH (r) under pressures of 5–1000 atm.

In Figure 8a, the diffusion coefficients of water molecule, HCO3− and Na+, or DO(water), DC(HCO3-), Figure 7. (a) gCNa(r) under pressures of 5–1000 atm; and (b) gCH(r) under pressures 5–1000 atm. + , of and DNa , obtained from MSD defined byofEquation (11), are plotted pressure. The obtained In Figure 8a, the diffusion coefficients water molecule, HCO3 ´against and Na or DO(water) , DC(HCO3-) , values agree well those in the literature [27]. As seen in Figure 8a, all D i s decrease as the pressure and D , obtained from MSD defined by Equation (11), are plotted against pressure. The obtained Na In Figure 8a, the diffusion coefficients of water molecule, HCO3− and Na+, or DO(water), DC(HCO3-) , increases. It is noteworthy that D C(HCO3-), and DNa show similar pressure dependence. These features values agree those in the [27]. As seen inare Figure 8a,against all Di spressure. decrease as obtained the pressure and DNa, well obtained from MSDliterature defined by Equation (11), plotted The of gij(r) s and Di s suggest that the complex [HCO3·(H2O)n]− is expected to be formed in the solution. increases. is noteworthy DC(HCO3-) show similar8a, pressure features valuesItagree well those that in the literature, and [27]. D As in Figure all Di s dependence. decrease as theThese pressure Naseen −} should be compounded to hold the Then, the clusters {Na+·[HCO3·(H2O)n]−} and {H+·[HCO3·(H2O) ´n]is of gij (r) s and D thatthat theDcomplex [HCO ¨ (H2 O) ] expected to be formed in the solution. increases. Iti iss suggest noteworthy C(HCO3-), and DNa 3show similar pressure dependence. These features n local charge neutrality. Similar structures also +been found in the aqueous solutions. According + ¨ [HCO ´ }have ´ }beshould gij(r)clusters s and D{Na i s suggest that3 ¨the complex [HCO{H 3·(H¨2[HCO O)n]− is3expected formedbe in compounded the solution. to Then,of the (H O) ] and ¨ (H2 O)nhydration ]to n 2 + to the ab initio MD study of Na in aqueous solution, the n-coordinate structures, such as +·[HCO3·(H2O)n]−} and {H+·[HCO3·(H2O)n]−} should be compounded to hold the Then, the clusters {Na hold Na(H the local charge neutrality. Similar structures have also been found in the aqueous solutions. 2O)n+, have been found [28]. In an aqueous solution of CaCO3, Ca(HCO3)2(H2O)4 and local charge neutrality. Similar structures have also been found in then-coordinate aqueous solutions. According + in According to the study aqueous solution, the hydration structures, Ca(HCO3)3(Hab 2O)initio 2− are MD predicted toofbeNa stable [29]. to the ab initio +MD study of Na+ in aqueous solution, the n-coordinate hydration structures, such as such as Na(H2+ O)n , have been found [28]. In an aqueous solution of CaCO3 , Ca(HCO3 )2 (H2 O)4 and Na(H2O)n , have been found [28]. In an aqueous solution of CaCO3, Ca(HCO3)2(H2O)4 and Ca(HCO3 )3 (H2 O)2 ´ −are predicted to be stable [29]. Ca(HCO3)3(H2O)2 are predicted to be stable [29].

Figure 8. (a) Di under pressures of 5–1000 atm; and (b) D

Figure 8. (a) Di under pressures of 5–1000 atm; and (b) D

(

) (ω) under

pressures of 5–1000 atm.

(ω) under pressures of 5–1000 atm.

( ) Figure 8. (a) Di under pressures of 5–1000 atm; and (b) DCpHCO ´ pωq under pressures of 5–1000 atm. q 3

Int. J. Mol. Sci. 2016, 17, 45

8 of 18

Int. J. Mol. Sci. 2016, 17, 45

7 of 17

The frequency dependent diffusion coefficient, (ω), D pωq, can be derived from the velocity The frequency dependent diffusion coefficient, D i can be derived from the velocity autoxv ptq ¨ v p0qy auto-correlation function (VAF), or using Equation (12). The obtained D ´ pωq for i correlation function (VAF), or 〈 ( i) ∙ (0)〉 using Equation (12). The obtained DCpHCO 3 )q(ω) for ( ´ HCO under various pressures are observed in the THz or the infrared region as shown in Figure 8b. 8b. HCO33− under various pressures are observed in the THz or the infrared region as shown in Figure pωqaround The peak of of D DCpHCO´(ω) around 300 1/cm under pressure is very the frequency of The peak 300 1/cm under lowlow pressure is very closeclose to theto frequency of water ( )3 q water caused by the translational cage effect. The peak position slightly sifts to the higher frequency caused by the translational cage effect. The peak position slightly sifts to the higher frequency around around 500 1/cm, and the small hump1000 around 1000 can be observed, are comparable 500 1/cm, and the small hump around 1/cm can1/cm be observed, which arewhich comparable to the COto 2– the CO –H O intermolecular vibrational frequency [30]. 2 2 H2O intermolecular vibrational frequency [30]. ´ we calculate the rotational In 3 ¨ 3(H In order order to to evaluate evaluate the thelifetime lifetimeofofthe thecomplex complex[HCO [HCO ·(H2 2O) O)nn]]−, , we calculate the rotational ´ correlation function, C (t), of HCO ion defined by Equation (19). The C (t) is thought to to be be affected by − 2 3 2 correlation function, C2(t), of HCO3 ion defined by Equation (19). The C2(t) is thought affected ´ and the surrounding water molecules. The logarithm the relaxation of the interaction between HCO − 3 HCO3 and the surrounding water molecules. The by the relaxation of the interaction between plot of the obtained C2 obtained (t)s in various are shown in Figure 9. The can be evaluated from logarithm plot of the C2(t)spressure in various pressure are shown in lifetime Figure 9. The lifetime can be ´ at 5 atm is extremely the inclination of the linear part of lnC (t). The graph of lnC (t) of HCO similar − evaluated from the inclination of the 2linear part of lnC2(t).2The graph 3of lnC2(t) of HCO3 at 5 atm is to that of water molecule. oscillatory at 3–4 ps may beatinterpreted “free-rotor” extremely similar to that ofThe water molecule.behavior The oscillatory behavior 3–4 ps mayasbethe interpreted as motion, which is observed in the dilute phase [15]. The estimated lifetime at 5 atm is 1.6 ps; on the other the “free-rotor” motion, which is observed in the dilute phase [15]. The estimated lifetime at 5 atm is hand, at 100–1000 atm is ˘ 1.3 ps, which is comparable to the relaxation time of H-bond in 1.6 ps;those on theofother hand, those of6.7 at 100–1000 atm is 6.7 ± 1.3 ps, which is comparable to the relaxation the carbonate solution carbonate [31]. The slight positive of the relaxation time timeaqueous of H-bond in the aqueous solution [31]. pressure The slightdependence positive pressure dependence of also suggests the structure formation in the higher-pressure region. the relaxation time also suggests the structure formation in the higher-pressure region.

− under pressures of 5–1000 atm with logC2(t) of water under 5 atm. Figure9.9.lnC lnC2(t) (t) of of HCO HCO3´ Figure 2 3 under pressures of 5–1000 atm with logC2 (t) of water under 5 atm.

For the next stage, according to the G-K formula, the shear viscosity is calculated using Equation For the next stage, according to the G-K formula, the shear viscosity is calculated using (17). The shear viscosities obtained by MD are shown in Figure 10a with the experimental values for Equation (17). The shear viscosities obtained by MD are shown in Figure 10a with the experimental pure water, 0.6 M NaCl aqueous solution, and 0.6 M NaCl and 0.913 M CO2 aqueous solution in the values for pure water, 0.6 M NaCl aqueous solution, and 0.6 M NaCl and 0.913 M CO2 aqueous solution literature [32,33]. The pressure dependence of the present result (0.6 M NaCl with saturated HCO3−) in the literature [32,33]. The pressure dependence of the present result (0.6 M NaCl with saturated is positive, whereas those of experimental values are negative. This fact suggests that the interaction HCO3 ´ ) is positive, whereas those of experimental values are negative. This fact suggests that the between HCO3− and water,´and/or other constituents increases as the pressure increases. A significant interaction between HCO3 and water, and/or other constituents increases as the pressure increases. increasing tendency of viscosity with increasing mole fraction of dissolved CO2 has also been A significant increasing tendency of viscosity with increasing mole fraction of dissolved CO2 has also observed in the viscosity measurement of CO2 saturated seawater at 303 to 333 K under constant been observed in the viscosity measurement of CO2 saturated seawater at 303 to 333 K under constant pressure 10 to 20 MPa [34]. To ensure the MD result, the shear viscosity is also estimated from the pressure to 20 MPa [34]. Tobyensure the MD result, the shear(S-E) viscosity is also the diffusion 10 coefficient obtained MD using the Stokes–Einstein relation for aestimated sphericalfrom particle, diffusion coefficient obtained by MD using the Stokes–Einstein (S-E) relation for a spherical particle, which is expressed as which is expressed as k T k T = B (2) D “= B “ (2)  6πη ξ 6πηr

Int. J. Mol. Sci. 2016, 17, 45

9 of 18

Int. J. Mol. Sci. 2016, 17, 45

8 of 17

where the parameters ξ and η stand for the friction constant and the shear viscosity, respectively. If the where the parameters ξ and ηat stand for the constant cand theη(c) shear respectively. If shear viscosity η is determined a certain COfriction at viscosity, any concentration c could 2 concentration 0 , then the shear viscosity η isfollowing determined at a certain be estimated using the equation [35]: CO2 concentration c0, then η(c) at any concentration c could be estimated using the following equation [35]: D pcq η pc0 q (3) ( ) “ η( ) D pc0 q = η pcq (3) ( ) η( ) which is known as the Walden’s rule. The calculated η(c), from Equation (3),(3), is also plotted in Figure 10a, which is known as the Walden’s rule. The calculated η(c), from Equation is also plotted in Figure which agreesagrees with the MD to a certain extent. 10a, which with theresults MD results to a certain extent.

Figure 10. (a) Viscosity under pressures of 5–1000 atm; and (b) Thermal conductivity under pressures Figure 10. (a) Viscosity under pressures of 5–1000 atm; and (b) Thermal conductivity under pressures of 5–1200 atm. of 5–1200 atm.

As mentioned in the previous subsection, we have calculated the thermal conductivity of As mentioned in the solution previouswith subsection, calculated the thermal conductivity of 1.1 mol % NaCl aqueous saturatedwe COhave 2 by NEMD method [7,16]. In this study, we 1.1 mol % NaCl aqueous solution with saturated CO by NEMD method [7,16]. In thissolution. study, we adopt the same method using the saturated HCO3−2ion in 1.1 mol % NaCl aqueous Asadopt will ´ the same method using3,the HCO3 ionλin mol % as NaCl aqueous solution. will be be described in Section thesaturated thermal conductivity is 1.1 expressed Equations (14)–(16). TheAs obtained described the thermal by conductivity is expressed as10b. Equations (14)–(16). The results of in theSection thermal3,conductivity NEMD are λshown in Figure The experimental dataobtained of pure results theseawater thermalare conductivity byinNEMD in The Figure 10b. The experimental dataof of water of and also shown Figureare 10bshown [23,36]. negative pressure dependence pure water and seawater are alsoseen shown in MD Figure 10b [23,36]. The negative pressure of thermal conductivity is clearly in the result; on the other hand, those of thedependence experimental thermal conductivity is clearly seen in the MD result; on the other hand, those of the experimental data are positive. data are Aspositive. stated already, some anomalous results have been obtained by MD in the transport and −. The As stated already, some results solution have been obtained by HCO MD 3in the transport and thermal properties of 1.1 molanomalous % NaCl aqueous saturated with experimental ´ . The experimental thermal properties of 1.1 mol % NaCl aqueous solution saturated with HCO thermal conductivity data of electrolyte aqueous solutions show positive pressure dependence, and 3 thermal of electrolyte aqueous[37]. solutions positive pressure dependence, and negativeconductivity concentrationdata dependence of electrolyte These show phenomena have been explained to some negative concentration dependence of electrolyte These phenomena been explained to some extent by the extension of the additivity of the [37]. thermal conductivity byhave considering the interaction between components results of in the thisthermal study might be influenced by the above-mentioned extent by the extension[37,38]. of the The additivity conductivity by considering the interaction contradictory effects to the thermal conductivity, pressure andbe concentration. Inthe addition, the results between components [37,38]. The results in this study might influenced by above-mentioned are also supposed toto bethe attributed the complexpressure and/or the cluster formationIn in addition, the solution. contradictory effects thermalto conductivity, and concentration. the results are also supposed to be attributed to the complex and/or the cluster formation in the solution. 2.3. The Methane and NaCl Aqueous Solution 2.3. The Methane and NaCl Aqueous Solution Next, we will show the MD result of the methane and NaCl aqueous solution. The fundamental Next, we MDasresult of theinmethane and NaCl aqueousThe solution. fundamental procedure of will MD show is thethe same described the previous subsections. water The (TIP4P) and the procedure of MD is the samebody as described inThe the concentration previous subsections. waterto(TIP4P) and the methane are treated as rigid molecules. of NaCl isThe adjusted be the same as that of seawater, 1.1as mol %. The of CH 4 inconcentration the MD cell is of determined using the data methane are treated rigid bodynumber molecules. The NaCl is adjusted tosolubility be the same as of CH 4 in seawater [39]. %. TheThe numbers ofof particles in MD listed in theusing Tablethe 5. solubility data that of seawater, 1.1 mol number CH4 in used the MD cell are is determined of CH4 in seawater [39]. The numbers of particles used in MD are listed in the Table 5.

Int. J. Mol. Sci. 2016, 17, 45 Int. J. Mol. Sci. 2016, 17, 45

10 of 18 9 of 17

Table5.5.The Thenumber numberofofCH CH4and andNaCl NaClmolecules moleculesused usedin inthe theaqueous aqueoussolution solutionatat275 275K. K. Table 4

Pressure Water +Na+ Cl´− Water Na Cl 10 MP 10,000 112 112 10 MP 112 30 MP 10,000 10,000112112 112 30 MP 10,000 112 112 60 MP 10,000 10,000112112 112 60 MP 112 100 MP10,000 10,000112112 112 100 MP 112

Pressure

CH4 CH 26 4 46 26 46 63 63 82 82

From the MD results, the pair distribution functions, gij(r)s, and the integrated coordination From the MD results, the pair distribution functions, gij (r)s, and the integrated coordination numbers, nij(r)s, have been obtained for 10 to 100 MP, which are shown in Figure 11a–d. Although numbers, nij (r)s, have been obtained for 10 to 100 MP, which are shown in Figure 11a–d. Although the the pressure dependence of gij(r) is not large, the slight change of the first peak height for CH4–CH4, pressure dependence of gij (r) is not large, the slight change of the first peak height for CH4 –CH4 , and and the depth for the first minimum for H2O–H2O can be observed, which may correspond to the the depth for the first minimum for H2 O–H2 O can be observed, which may correspond to the cage cage formation in the solution. The water coordination number, nij(r), of the first hydration shell formation in the solution. The water coordination number, nij (r), of the first hydration shell around the around the solute is calculated to the first minimum of gij(r) using Equation (9). The obtained water solute is calculated to the first minimum of gij (r) using Equation (9). The obtained water coordination coordination numbers calculated under pressures of 10–100 MP are listed in Table 6. The slight numbers calculated under pressures of 10–100 MP are listed in Table 6. The slight decreasing tendency decreasing tendency of water molecules around CH4 has been detected, which might also be of water molecules around CH4 has been detected, which might also be attributed to the cluster attributed to the cluster formation around CH4 molecules. formation around CH4 molecules.

(a)

(b)

(c)

(d)

Figure 11. 11. (a) (r) and nOwOw (b) gNaOw(r) and nNaOw(r); (c) gClOw(r) and nClOw(r) ; and (d) gCH4Ow(r) Figure (a)gOwOw gOwOw (r) and n(r); OwOw (r); (b) gNaOw (r) and nNaOw (r); (c) gClOw (r) and nClOw (r) ; and n CH4Ow(r) under pressures of 10–100 MP. and (d) gCH4Ow (r) and nCH4Ow (r) under pressures of 10–100 MP.

Int. J. Mol. Sci. 2016, 17, 45

11 of 18

Int. J. Mol. Sci. 2016, 17, 45

10 of 17

Table 6. Water coordination Table 6. Water coordination number number under under pressures pressures of of 10–100 10–100 MP. MP. + + Pressure Water Water Na Na Pressure

10 MP 10 MP 30 MP 30 MP 60 MP 60 MP 100 MP 100 MP

4.56 4.56 4.56 4.56 4.56 4.56 4.48 4.48

5.50 5.50 5.53 5.53 5.69 5.69 5.77 5.77

Cl Cl−´ 6.57 6.57 6.58 6.58 6.58 6.58 6.58 6.58

CHCH 4 4 20.58 20.58 20.25 20.25 19.95 19.95 19.65 19.65

Finally,the the pressure pressuredependence dependenceof of thermal thermal conductivity conductivity of of methane methane and and NaCl NaCl aqueous solution Finally, is obtained obtainedbybyNEMD NEMD using Equations (14)–(16). The obtained are inshown 12 is using Equations (14)–(16). The obtained results results are shown Figurein 12 Figure alongside alongside the experimental MD data for solution NaCl aqueous solution and negative pure water. The the experimental data and MDdata dataand for NaCl aqueous and pure water. The pressure negative pressure dependence of thermal conductivity higher pressure is also result dependence of thermal conductivity in higher pressure isinalso observed. This resultobserved. might beThis attributed might be attributed to the change or the clathrate formation around in thethe CHhigh-pressure 4 molecule in to the structure change or structure the clathrate formation around the CH4 molecule the high-pressure region, with which consistentregarding with thethe discussion regarding the decreasing of region, which is consistent theisdiscussion decreasing of coordination number in coordination the solution. number in the solution.

Figure 12. The pressure dependence of the thermal conductivity of CH4 and NaCl aqueous solution, Figure 12. The pressure dependence of the thermal conductivity of CH4 and NaCl aqueous solution, and NaCl aqueous obtained by molecular dynamics (MD), alongside the experimental data. and NaCl aqueous obtained by molecular dynamics (MD), alongside the experimental data.

3. Procedure 3. Procedure The essential numerical procedure of this MD simulation study is the same as our previous The essential numerical procedure of this MD simulation study is the same as our previous works works on aqueous solutions [7,8,12,40]; however, the fundamental part of the procedure is described on aqueous solutions [7,8,12,40]; however, the fundamental part of the procedure is described as as follows for the reader’s benefit. follows for the reader’s benefit. The water is treated as the rigid body model, TIP4P [17]. The potential function for TIP4P is The water is treated as the rigid body model, TIP4P [17]. The potential function for TIP4P is expressed in the charged Lennard-Jones (L-J) type potentials as expressed in the charged Lennard-Jones (L-J) type potentials as σ σ (4) "´ ¯ + ´ ¯ *  ( ) = 2 + 4ε 12 zi z j e σ σ 6 φij prq “ ` 4ε ` (4) r r In the equations hereafter, i and j stand for ther constituent atoms; zi is the charge for the constituent species i;hereafter, and e is the elementary The usedatoms; parameters listedfor inthe Table 7. In the equations i and j stand forcharge. the constituent zi is theare charge constituent species i; and e is the elementary charge. The used parameters are listed in Table 7.

Int. J. Mol. Sci. 2016, 17, 45

12 of 18

Table 7. The potential parameters for TIP4P.

O–O O–XX O–H XX–XX XX–H H–H

ε (gA2 /fs2 )

zj

zi

The Atom Pair

0 0 0 ´1.04 ´1.04 0.52

0 ´1.04 0.52 ´1.04 0.52 0.52

1.0769 ˆ 0 0 0 0 0

10´28

σ (A) 3.153650 0 0 0 0 0

The interactions between Na+ and Cl´ , TIP4P-Na+ , and TIP4P-Cl´ are expressed as [41], φij prq “

z i z j e2 D C ` 9´ 6 r r r

(5)

The used parameters are listed in Table 8. Table 8. The potential parameters between Na+ , Cl´ and TIP4P water. The Atom Pair

zi

zj

C (10´19 JA9 )

D (10´19 JA6 )

Na+ –Na+ Cl´ –Cl´ Na+ –Cl´ Na+ –O Na+ –H Na+ –XX Cl´ –O Cl´ –H Cl´ –XX

1.0 ´1.0 1.0 1.0 1.0 1.0 ´1.0 ´1.0 ´1.0

1.0 ´1.0 ´1.0 0 0.52 ´1.04 0 0.52 ´1.04

1.5565 ˆ 102 2.4808 ˆ 104 3.2015 ˆ 103 5.8553 ˆ 102 0 0 9.2132 ˆ 103 0 0

8.3931 2.6162 ˆ 102 4.6860 ˆ 10 2.3463 ˆ 10 0 0 1.3099 ˆ 102 0 0

The interactions between other solute ions and water molecules are also expressed in the charged L-J type potentials [42,43] as z i z j e2 A B ` 12 ´ 6 (6) φij prq “ r r r The parameters of the interactions between CO2 molecule, HCO3 ´ ion, Na+ , Cl´ and water molecule, taken from the literature, are listed in Tables 9–12. Table 9. The potential parameters between CO2 –CO2 , CO2 –Na+ , CO2 –Cl´ , and CO2 –TIP4P water. The Molecule Pair

The Atom Pair

zi

zj

A (kcal A12 /mol)

B (kcal A6 /mol)

CO2 –CO2

C 2–C 2 C 2–O 2 O 2–O 2

0.4578 0.4578 ´0.2289

0.4578 ´0.2289 ´0.2289

3.2481 ˆ 106 1.1110 ˆ 106 3.8000 ˆ 105

1.1680 ˆ 103 8.1233 ˆ 102 5.6498 ˆ 102

CO2 –Na+

C 2–Na+ O 2–Na+

0.4578 ´0.2289

1.0000 1.0000

7.9212 ˆ 105 3.3217 ˆ 105

8.3121 ˆ 102 5.3912 ˆ 102

CO2 –Cl´

C 2–Cl´ O 2–Cl´

0.4578 ´0.2289

´1.0000 ´1.0000

3.4465 ˆ 106 1.1789 ˆ 106

1.5642 ˆ 103 1.0879 ˆ 103

CO2–TIP4P

C 2–O 1 C 2–XX C 2–H O 2–O 1 O 2–XX O 2–H

0.4578 0.4578 0.4578 ´0.2289 ´0.2289 ´0.2289

0.0000 ´1.04 0.5200 0.0000 ´1.04 0.5200

1.0834 ˆ 106 0 1.9677 ˆ 105 3.7057 ˆ 105 0 6.7305 ˆ 104

7.6092 ˆ 102 0 2.3881 ˆ 102 5.2922 ˆ 102 0 1.6609 ˆ 102

Int. J. Mol. Sci. 2016, 17, 45

13 of 18

Table 10. The potential parameters between H+ –H+ , HCO3 ´ –HCO3 ´ , TIP4P–H+ , TIP4P–HCO3 ´ , Na+ –H+ , and Na+ –HCO3 ´ . The Atom Pair

zi

zj

A (kcal A12 /mol)

B (kcal A6 /mol)

H+ –H+

H+ –H+

1.0000

1.0000

1.7199 ˆ 104

3.2337 ˆ 101

HCO3 ´ –HCO3 ´

C 2–C 2 C 2–O 2 C 2–O 1 C 2–H O 2–O 2 O 2–O 1 O 2–H O 1–O 1 O 1–H H–H

1.2149 1.2149 1.2149 1.2149 ´0.8727 ´0.8727 ´0.8727 ´0.9424 ´0.9424 0.4194

1.2149 ´0.8727 ´0.9424 0.4194 ´0.8727 ´0.9424 0.4194 ´0.9424 0.4194 0.4194

1.1713 ˆ 106 5.3596 ˆ 105 5.3596 ˆ 105 1.5060 ˆ 105 2.3212 ˆ 105 2.3212 ˆ 105 6.3567 ˆ 104 2.3212 ˆ 105 6.3567 ˆ 104 1.7199 ˆ 104

6.6752 ˆ 102 4.5224 ˆ 102 4.5224 ˆ 102 1.5134 ˆ 102 2.9808 ˆ 102 2.9808 ˆ 102 9.8477 ˆ 101 2.9808 ˆ 102 9.8477 ˆ 101 3.2337 ˆ 101

TIP4P–H+

O 1–H+ XX–H+ H–H+

0 ´1.0400 0.5200

1.0000 1.0000 1.0000

6.3567 ˆ 104 0 1.7199 ˆ 104

9.8477 ˆ 101 0 3.2337 ˆ 101

TIP4P–HCO3 ´

O 1–C 2 O 1–O 2 O 1–O 1 O 1–H XX–C 2 XX–O 2 XX–O 1 XX–H H–C 2 H–O 2 H–O 1 H–H

0 0 0 0 ´1.0400 ´1.0400 ´1.0400 ´1.0400 0.5200 0.5200 0.5200 0.5200

1.2149 ´0.8727 ´0.9424 0.4194 1.2149 ´0.8727 ´0.9424 0.4194 1.2149 ´0.8727 ´0.9424 0.4194

5.3596 ˆ 105 2.3212 ˆ 105 2.3212 ˆ 105 6.3567 ˆ 104 0 0 0 0 1.5060 ˆ 105 6.3567 ˆ 104 6.3567 ˆ 104 1.7199 ˆ 104

4.5224 ˆ 102 2.9808 ˆ 102 2.9808 ˆ 102 9.8477 ˆ 101 0 0 0 0 1.5134 ˆ 102 9.8477 ˆ 101 9.8477 ˆ 101 3.2337 ˆ 101

Na+ –H+

Na+ –H+

1.0000

1.0000

8.9597 ˆ 104

1.7676 ˆ 102

Na+ –HCO3 ´

Na+ –C 2 Na+ –O 2 Na+ –O 1 Na+ –H

1.0000 1.0000 1.0000 1.0000

1.2149 ´0.8727 ´0.9424 0.4194

7.9212 ˆ 105 3.3217 ˆ 105 3.3217 ˆ 105 8.9597 ˆ 104

8.3121 ˆ 102 5.3912 ˆ 102 5.3912 ˆ 102 1.7676 ˆ 102

The Molecule Pair

Table 11. The potential parameters between Cl´ –H+ , Cl´ –HCO3 ´ , and H+ –HCO3 ´ . The Molecule Pair

The Atom Pair

zi

zj

A (kcal A12 /mol)

B (kcal A6 /mol)

Cl´ –H+

Cl´ –H+

´1.0000

1.0000

2.0880ˆ105

3.1983 ˆ 102

Cl´ –HCO3 ´

Cl´ –C 2 Cl´ –O 2 Cl´ –O 1 Cl´ –H

´1.0000 ´1.0000 ´1.0000 ´1.0000

1.2149 ´0.8727 ´0.9424 0.4194

3.4465 ˆ 106 1.1789 ˆ 106 1.1497 ˆ 106 2.0880 ˆ 105

1.5642 ˆ 103 1.0879 ˆ 103 1.0191 ˆ 103 3.1983 ˆ 102

H+ –HCO3 ´

H+ –C 2 H+ –O 2 H+ –O 1 H+ –H

1.0000 1.0000 1.0000 1.0000

1.2149 ´0.8727 ´0.9424 0.4194

1.9677 ˆ 105 6.7305 ˆ 104 6.5635 ˆ 104 1.1921 ˆ 104

2.3881 ˆ 102 1.6609 ˆ 102 1.5558 ˆ 102 4.8828 ˆ 10

Int. J. Mol. Sci. 2016, 17, 45

14 of 18

Table 12. The potential parameters between CH4 –CH4 , TIP4P–CH4 , Na+ –CH4 , and Cl´ –CH4 . The Molecule Pair

The Atom Pair

zi

zj

A (kcal A12 /mol)

B (kcal A6 /mol)

CH4 –CH4

C 1–C 1 C 1–H H–H

´0.3744 ´0.3744 0.0936

´0.3744 0.0936 0.0936

3.2481 ˆ 106 1.9677 ˆ 105 1.1921 ˆ 104

1.1680 ˆ 103 2.3881 ˆ 102 4.8828 ˆ 101

CH4 –TIP4P

C 1–O 1 H–O 1 C 1–XX H–XX C 1–H H–H

´0.3744 0.0936 ´0.3744 0.0936 ´0.3744 0.0936

0 0 ´1.0400 ´1.0400 0.5200 0.5200

1.0834 ˆ 106 6.5635 ˆ 104 0 0 1.9677 ˆ 105 1.1921 ˆ 104

7.6092 ˆ 102 1.5558 ˆ 102 0 0 2.3881 ˆ 102 4.8828 ˆ 101

CH4 –Na+

C 1–Na+ H–Na+

´0.3744 0.0936

1.0000 1.0000

7.9212 ˆ 105 8.9598 ˆ 104

8.3121 ˆ 102 1.7676 ˆ 102

CH4 –Cl´

C 1–Cl´ H–Cl´

´0.3744 0.0936

´1.0000 ´1.0000

3.4465 ˆ 106 2.0880 ˆ 105

1.5642 ˆ 103 3.1983 ˆ 102

The solute ions and molecules and water molecules are at first randomly placed at the lattice points, which are obtained by dividing each side by the number “n” that satisfies the following relation, pn ´ 1q3 ă pthe total number of water molecules and solute ions and moleculesq ă n3

(7)

Then, the relaxation of the first configuration using the Monte Carlo method is executed by the following procedure: (a) (b)

(c)

(d) (e) (f)

One ion or molecule is randomly selected. One degree of freedom is selected randomly for the set of {ξ} = {r, Ω, θ}, where “r” is the translation of the center of mass, “Ω” is the rotation around the center of mass, and “θ” is the rotation around the bond axis of the molecule. The selected degree of freedom ξ of the selected ion or molecule is changed by ∆ξ, then the now set {ξ}1 is created. The increment of freedom “∆r” is randomly determined from 0 to ∆rmax . A random degree is applied for “∆Ω” and “∆θ”. The above procedures (a)–(c) are repeated for the defined number of steps until the ensemble average is calculated. The potential energy difference ∆φ between the previous configuration {ξ} and the new configuration {ξ}1 is calculated. The decision whether the new configuration {ξ}1 is adopted or not is made according to the following condition: If ∆φ < 0, then the new configuration {ξ}1 is adopted, If ∆φ > 0, then a uniform random number η is compared to the Boltzmann factor exp(-∆φ /kT), If η ď exp(-∆φ/kT), then {ξ}1 is adopted as the new configuration, If η ą exp(-∆φ/kT), then {ξ}1 is not adopted.

The (a)–(f) procedures above are repeated to obtain the lower energy configuration, which we then adopt as the first configuration of the MD calculation. The cut of distance for the van der Waals interaction is 15 Å. The Ewald method is used for the calculation of the Coulomb interaction, in which the square of the cut off distance in the reciprocal lattice space is 27. The static properties of the solutions are calculated in the NTP constant condition for 100,000–500,000 steps with 0.1 fs being one time step [18–20]. The very short time step is adopted so as to detect the fast movement of H+ and water molecules.

Int. J. Mol. Sci. 2016, 17, 45

15 of 18

The pair distribution function, gij (r), can be firstly obtained from a time series data of coordinates of ions as [15] ˆ ˙ ´ Ni ¯ ` ˘ÿ ∆r ∆r gij prq “ V{Ni Nj ckj r ´ , r` { 4πr2 ∆r (8) 2 2 k

where Ni and N j are the numbers of the ion species i and j, respectively. V is the volume of the cell. cik is the number of ion k in the spherical shell of the thickness of ∆r at the distance r from the ion i. The distance dependent coordination number, or the integrated coordination number, nij (r), is defined as ÿ

nij prq “

pNj {Vqgij pn¨ ∆rq

(9)

năr{∆r

The interference functions for neutron Iij (q) are obtained from gij (r), which is expressed as [21] Iij pqq “ ci bi c j b j {

˜ ÿ

¸2 ż c k bk

8

0

k

´ ¯ sin prqq 4πr2 ρ0 gij prq ´ 1 dr qr

(10)

where ci is the atomic fraction of the i-type atoms; ρ0 is the average number density; and bi is the neutron scattering amplitude of the i-type atom. The diffusion coefficient for i-type atom, Di , can be obtained from MSD, which is defined as E 1 A |ri ptq ´ ri p0q|2 tÑ8 6t

Di “ lim

(11)

VAF is calculated to examine the dynamical and transport properties, which is expressed as xvi ptq ¨ vi p0qy. The frequency dependent diffusion coefficient, Di pωq, can be obtained from VAF as 1 Di pωq “ 3

ż8 xvi ptq ¨ vi p0qy cosωt dt

0

(12)

where ω = 2π f ; f is the frequency, and Di (ω = 0) = Di . The thermodynamic properties are also important for the study of solutions. The specific heat at the constant pressure is expressed as 1 dH 1 “ Cp “ M dT M

ˆ

˙ˇ dU dV ˇˇ `P dT dT ˇp

(13)

In order to calculate the thermal conductivity, NEMD method is used. In the NEMD method, the average of energy flux overtime is performed to avoid the margin of error caused by the integration of the energy flux autocorrelation function, which is used in the direct method or the EQ method [16]. The system reaches thermal equilibrium by EQMD, then the perturbation is applied at time t = 0. The thermal conductivity is expressed as the G-K formula [15], 1 λ“ Vk B T 2

ż8 xJx pτq Jx p0qy dτ

(14)

0

In Equation (14), Jx (τ) stands for the x component of the perturbation current; V the volume of the system; kB the Boltzmann constant; and T the temperature of the system. The applied perturbation Fext is related to the average of the perturbation current t , xJx yt “

Fext kB T

żt xJx pτq Jx p0qy dτ 0

(15)

Int. J. Mol. Sci. 2016, 17, 45

16 of 18

Then, the relation between the thermal conductivity and the perturbation is expressed as [15] λ“

1 lim xJx yt VFext T tÑ8

(16)

According to the G-K formula, the shear viscosity is expressed by the integration of the autocorrelation function of an off-diagonal element of the stress tensor in the long wave length limit k Ñ 0 [15], ż8 @ D 1 xz η“ lim σkxz ptq σ´ p0q (17) k kTV 0 kÑ0 The reorientation of linear molecules is expressed as the time-correlation function defined as [15] Cl ptq “ xPl rui ptq ¨ ui sy

(18)

where ui (t) is unit vector parallel to the principal axis of the molecule i. Pl (x) is the Legendre polynomial. The rotational correlation function is the case of l = 2, which is expressed as C2 ptq “ P2 pcosωtq “

E 1A E 1A 3 pcosωtq2 ´ 1 “ 3 tui ptq ¨ ui u2 ´ 1 2 2

(19)

In the case that the decay in polarization anisotropy C2 (t) is modeled in the standard exponential form C2 (t) = A exp(´t/τ), the relaxation time is obtained by the logarithm plot of C2 (t). The main part of MD is performed using SIGRESS ME package (Fujitsu) [44] at the supercomputing facilities in Kyushu University. The obtained results from MD simulation are summarized in the previous section. 4. Conclusions MD simulations have been performed for the system CO2 , HCO3 ´ , and CH4 dissolved in 1.1 mol % NaCl aqueous solution as a seawater model. The pressure dependence of structure, coordination number, diffusion coefficient, frequency distribution of ions, shear viscosity, heat capacity, and thermal conductivity have been examined. It is worth noting that the negative pressure dependence of thermal conductivity has been detected in these solutions, which is considered to be attributed to the structure change of the solutions under high pressure. Acknowledgments: The author thanks S. Tamaki for his helpful comments and encouragement on this study. This study was supported by JSPS KAKENHI Grant Number 15K05136, and the Mayekawa Houonkai Foundation. Part of the results in this study was obtained using the supercomputing facilities at Research Institute for Information Technology, Kyushu University. Conflicts of Interest: The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

References 1. 2. 3. 4. 5. 6.

Tans, P.; Keeling, R. Trends in Atmospheric Carbon Dioxide; NOAA Earth System Research Laboratory: Boulder, CO, USA. Available online: http://www.esrl.noaa.gov/gmd/ccgg/trends/ (accessed on 25 September 2015). Kohl, A.; Nielsen, R. Gas Purification, 5th ed.; Gulf Publishing Company: Houston, TX, USA, 1997. Caldeira, K.; Elderfield, H.; Hoegh-Guldberg, O.; Liss, P.; Riebesell, U.; Shepherd, J.; Turley, C.; Watson, A. Ocean Acidification Due to Increasing Atmospheric Carbon Dioxide; The Royal Society: London, UK, 2005. Huang, T. Molecular Dynamics Simulation of Carbon Dioxide in Aqueous Electrolyte Solution. Ph.D. Thesis, Swinburne University of Technology, Melbourne, Australia, 2012. Sabil, K.; Witkamp, G.J.; Peters, C.J. Phase equilibria of mixed carbon dioxide and tetrahydrofuran hydrates in sodium chloride aqueous solutions. Fluid Phase Equilibria 2009, 284, 38–43. [CrossRef] Yan, Y.; Chen, C.C. Thermodynamic modeling of CO2 solubility in aqueous solutions of NaCl and Na2 SO4 . J. Supercrit. Fluids 2010, 55, 623–634. [CrossRef]

Int. J. Mol. Sci. 2016, 17, 45

7. 8. 9. 10.

11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26.

27. 28. 29.

30.

31.

32.

17 of 18

Matsunaga, S. A molecular dynamics study of structure and thermal properties of carbon dioxide in sodium chloride aqueous solution. J. Phys. Conf. Ser. 2014, 490. [CrossRef] Matsunaga, S. Influence of HCO3 ´ ion on structure and transport properties of seawater. Trans. Mater. Res. Soc. Jpn. 2015, 40, 373–377. [CrossRef] Sloan, E.D.; Koh, C.A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press: Boca Raton, FL, USA, 2007. Tung, Y.T.; Chen, L.J.; Chen, Y.P.; Lin, S.T. Molecular dynamics study on the growth of structure I methane hydrate in aqueous solution of sodium chloride. J. Phys. Chem. B 2012, 116, 14115–14125. [CrossRef] [PubMed] Shvab, I.; Sadus1, R.J. Thermodynamic properties and diffusion of water + methane binary mixtures. J. Chem. Phys. 2014, 140. [CrossRef] [PubMed] Matsunaga, S. Effect of dissolution of methane in aqueous NaCl solution: A molecular dynamics study. JPS Conf. Proc. 2014, 1. [CrossRef] Matsunaga, S. Structural features in molten RbAg4 I5 by molecular dynamics simulation. Mol. Simul. 2013, 39, 119–122. [CrossRef] Matsunaga, S. Anomalous electrical properties in superionic (Agx Cu1´x )Br (x = 0.5): Ab initio study. Ionics 2015, 21, 161–166. [CrossRef] Hansen, J.P.; McDonald, I.R. Theory of Simple Liquids; Academic Press: London, UK, 1986. Yoshida, M.; Harada, A.; Takeuchi, A.; Matsunaga, K.; Matsubara, H. Perturbed molecular dynamics for calculating thermal conductivity of zirconia. Mol. Simul. 2004, 30, 953–961. [CrossRef] Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. [CrossRef] Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255–268. [CrossRef] Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511–519. [CrossRef] Andersen, H.C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72, 2384–2393. [CrossRef] Waseda, Y. The Structure of Non-Crystaline Materials; McGraw-Hill: New York, NY, USA, 1980. Petravic, J. Thermal conductivity of ethanol. J. Chem. Phys. 2005, 123, 1–7. [CrossRef] [PubMed] Abdulagatov, I.M.; Magomedov, U.B. Thermal conductivity of pure water and aqueous SrBr2 solutions at high temperatures and high pressures. High Temp. High Press 2003, 35–36, 149–168. [CrossRef] Nagasaka, Y.; Okaga, H.; Suzuki, J.; Nagashima, A. Absolute measurements of thermal conductivity of aqueous NaCl solutions at pressure up to 40 MPa. Ber. Bunsenges. Phys. Chem. 1983, 87, 859–866. [CrossRef] Duan, Z.; Sun, R. An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chem. Geol. 2003, 193, 257–271. [CrossRef] Dopieralski, P.D.; Burakowski, A.; Latajka, Z.; Olovsson, I. Hydration of NaHCO3 , KHCO3 , (HCO3 ´ )2 , HCO3 ´ and CO3 2´ from molecular dynamics simulation and speed of sound measurements. Chem. Phys. Lett. 2011, 507, 89–95. [CrossRef] Zeebe, R.E. On the molecular diffusion coefficients of dissolved CO2 , HCO3 ´ , and CO3 2´ and their dependence on isotopic mass. Geochim. Cosmochim. Acta 2011, 75, 2483–2498. [CrossRef] Rempe, S.B.; Pratt, L.R. The hydration number of Na+ in liquid water. Fluid Phase Equilibria 2001, 183, 121–132. [CrossRef] Tommaso, D.D.; de Leeuw, N.H. The onset of calcium carbonate nucleation: A density functional theory molecular dynamics and hybrid microsolvation/continuum study. J. Phys. Chem. B 2008, 112, 6965–6975. [CrossRef] [PubMed] Andersen, J.; Heimdal, J.; Mahler, D.W.; Nelander, B.; Wugt Larsen, R. THz absorption spectrum of the CO2 –H2 O complex: Observation and assignment of intermolecular van der Waals vibrations. J. Chem. Phys. 2014, 140, 1–5. [CrossRef] [PubMed] Kumar, P.P.; Kalinichev, A.G.; Kirkpatrick, R.J. Hydrogen-bonding structure and dynamics of aqueous carbonate species from car-parrinello molecular dynamics simulations. J. Phys. Chem. B 2009, 113, 794–802. [CrossRef] [PubMed] Harris, K.R.; Woolf, L.A. Temperature and volume dependence of the viscosity of water and heavy water at low temperatures. J. Chem. Eng. Data 2004, 49, 1064–1069. [CrossRef]

Int. J. Mol. Sci. 2016, 17, 45

33. 34. 35. 36. 37. 38. 39.

40. 41. 42. 43. 44.

18 of 18

Kumagami, A.; Yokoyama, C. Viscosities of aqueous NaCl solutions containing CO2 at high pressures. J. Chem. Eng. Data 1999, 44, 227–229. [CrossRef] Bando, S.; Takemura, F.; Nishio, M.; Hihara, E.; Akai, M. Viscosity of aqueous NaCl solutions with dissolved CO2 at (30 to 60) ˝ C and (10 to 20) MPa. J. Chem. Eng. Data 2004, 49, 1328–1332. [CrossRef] Atkins, P.W. Physical Chemistry, 4th ed.; Oxford University Press: Oxford, UK, 1990. Castelli, V.; Stanly, E.M.; Fischer, E.C. The thermal conductivity of seawater as a fuction of pressure and temperature. Deep Sea Res. 1974, 21, 311–319. Wang, P.; Anderko, A. Modeling thermal conductivity of electrolyte mixtures in wide temperature and pressure ranges: Seawater and its main components. Int. J. Thermophys. 2012, 33, 235–258. [CrossRef] Riedel, L. The heat conductivity of aqueous solutions of strong electrolytes. Chem. Ing. Tech. 1951, 23, 59–64. [CrossRef] Duan, Z.; Møller, N.; Greenberg, J.; Weare, J.H. The prediction of methane solubility in natural waters to high ionic strength from 0 to 250 ˝ C and from 0 to 1600 bar. Geochim. Cosmochim. Acta 1992, 56, 1451–1460. [CrossRef] Mataunaga, S.; Tamaki, S. Ionic conduction in electrolyte solution. J. Solut. Chem. 2014, 43, 1771–1790. [CrossRef] Zhengwei, P.; Ewig, C.S.; Hwang, M.J.; Waldman, M.; Hagler, A.T. Comparison of simple potential functions for simulating liquid water. J. Phys. Chem. A 1997, 101, 7243–7252. Jorgensen, W.L.; Laird, E.R.; Nguyen, T.B.; Tirado-Rives, J. Monte Carlo simulations of pure liquid substituted benzenes with OPLS potential functions. J. Comput. Chem. 1993, 14, 206–215. [CrossRef] Mayo, S.L.; Olafson, B.D.; Goddard, W.A., III. A Generic force field for molecular simulations. J. Phys. Chem. 1990, 94, 8897–8909. [CrossRef] FUJITSU Technical Computing Solution SCIGRESS. Available online: http://www.fujitsu.com/global/ solutions/business-technology/tc/sol/scigress/index.html (accessed on 25 September 2015). © 2015 by the author; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons by Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).