Effect of Greenhouse Gases Dissolved in Seawater

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.


Introduction
The reason for the global warming of the earth's climate seems undoubtedly attributable to the increasing of greenhouse gases, especially carbon dioxide (CO 2 ) in the atmosphere. The consumption of fossil fuels since the industrial revolution has released a large amount of CO 2 into the air, the concentration of which is expected to exceed 400 ppm in 2015 [1]. Therefore, technologies for reducing CO 2 emissions are a matter of the utmost importance, and have widely been investigated, including CO 2 absorbing processes into liquids [2]. In addition to the soluble property of CO 2 in water, the increasing concentration of CO 2 in the atmosphere and the warmer climate may accelerate the dissolution of CO 2 into seawater, which can simply be modeled as a sodium chloride (NaCl) aqueous solution [3].
However, few molecular dynamics (MD) studies on CO 2 and NaCl aqueous solutions are available in the literature [4]. Although many experimental studies on CO 2 dissolved into a NaCl aqueous solution have been published, the solubility of CO 2 and/or the stability of the solution are their main purpose [4][5][6]. To the best of our knowledge, the thermal conductivity of the CO 2 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 CO 2 absorbed into a NaCl aqueous solution. The CO 2 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 CO 2 molecules, however, react with water molecules to create bicarbonate (HCO 3´) and carbonate (CO 3 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 CO 3 2´i ons increases in alkaline regions, whereas CO 2 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 HCO 3´, i.e., the following reaction is postulated [3]: 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 HCO 3´i ons as a more realistic model [8].
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].
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 of the sea [9]. The dissolution of the methane into seawater and into the atmosphere also seems to affect the thermal environment of the earth. Many MD studies have been reported on the methane and water system; however, the nucleation or the structural properties in the mixture with pure water are their main purpose [10,11]. In this study, we wish to show the results of a MD simulation of the methane and NaCl aqueous solution system, as a model of methane dissolved into seawater. The change of structure, the coordination number, and the thermal conductivity will be discussed under the various 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.

Results and Discussion
In this section, first the results of MD simulation are described for the CO2 and NaCl aqueous solution system; then for the system that includes water molecules, Na + , Cl − , H + , and HCO3 − ions; and finally for the methane and NaCl aqueous solution system.

The CO2 and NaCl Aqueous Solution
In the study for CO2 and NaCl aqueous solution, the equilibrium MD (EQMD) is used to obtain the structure, the coordination number, the mean square displacement, and the specific heat. The fundamental procedure of EQMD is same as that in our previous works dealing with molten salts [13,14]. As mentioned in the previous section, the thermal conductivity is evaluated using the nonequilibrium 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 calculation is executed in the number of constituent molecules, the temperature, and the pressure (NTP) constant condition [18][19][20]. The applied pressure varies from 0.5 to 100 MP, or 100 × 10 P, which is equivalent to an ocean depth of 40 to 10,000 m. MD is carried out for 50,000 steps with 0.1 fs time interval. The concentration of NaCl in the water is the same as that of seawater, or 1.1 mol % NaCl. The concentration of CO2 is adjusted to that of a saturated solution. The used molecular numbers for the structure calculation at 275 K are listed in Table 1. The parameters of the system are monitored for 50,000 steps under 30 MP at 275 K to guarantee convergence ( Figure 1). About 5000 molecules (15,000 atoms) are used for the thermal conductivity calculation so as to contribute to rapid convergence.   The pair distribution functions, g ij (r), and the integrated coordination number, n ij (r), have been obtained from the MD results, and are shown in Figure 2a- To examine the structural change by pressure more clearly, we have calculated the coordination number, or the value of n ij (r) at the first minimum of g ij (r). The obtained coordination numbers are listed in Table 2. As a matter of fact, the pressure dependence of the coordination numbers, nNaO w (r) and nClO w (r), are in the margin of error, although a slight decreasing tendency can be observed. The decreasing coordination number for nC CO2 O w (r) is obviously seen. These facts might be evidence of the cage formation around the CO 2 molecules, because the cage structure yields the space around CO 2 molecules, thereby decreasing the coordination numbers.  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 enhancement of the first peak heights of gij(r) of CO2-H2O, Na + -H2O, Cl − -H2O and H2O-H2O can be observed under higher pressures, which may correspond to the cage formation by H2O molecules. To examine the structural change by pressure more clearly, we have calculated the coordination number, or the value of nij(r) at the first minimum of gij(r). The obtained coordination numbers are listed in Table 2. As a matter of fact, the pressure dependence of the coordination numbers, nNaOw(r) and nClOw(r), are in the margin of error, although a slight decreasing tendency can be observed. The decreasing coordination number for nCCO2Ow(r) is obviously seen. These facts might be evidence of the cage formation around the CO2 molecules, because the cage structure yields the space around CO2 molecules, thereby decreasing the coordination numbers.   The interference functions, I ij (q), which is also obtainable from the neutron diffraction experiment, can be calculated from the Fourier transformation of the pair distribution function g ij (r) as expressed in Equation (10) [21]. The evaluated I ij (q)s from g ij (r)s obtained by MD are shown in Figure 3a,b. A sharp peak of I C-Ow (q) at 2 [Å´1] at 283 K, 0.5 MP in Figure 3a is extremely enhanced at 275 K, 100 MP in Figure 3b. Although the peak heights of I Na-Ow (q) and I Cl-Ow (q) are similar to that of I C-Ow (q) at 283 K, 0.5 MP in Figure 3a, their heights decrease and their widths are broadened at 275 K, 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 CO 2 increases.  The interference functions, Iij(q), which is also obtainable from the neutron diffraction experiment, can be calculated from the Fourier transformation of the pair distribution function gij(r) as expressed in Equation (10) [21]. The evaluated Iij(q)s from gij(r)s obtained by MD are shown in 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 275 K, 100 MP in Figure 3b. Although the peak heights of INa-Ow(q) and ICl-Ow(q) are similar to that of IC-Ow(q) at 283 K, 0.5 MP in Figure 3a, their heights decrease and their widths are broadened at 275 K, 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. To examine the pressure dependence of the transport properties, the diffusion coefficients of constituent molecules and ions are calculated, which are obtainable from the inclination of the mean square displacement (MSD) using Equation (11). The obtained MSD for 30 MP is shown in Figure 4. Although small oscillations still remain in the diffusion coefficient for Cl − , the inclinations of MSDs seem to converge until 5 × 10 s (5 ps) and the slopes are kept for longer periods, to a certain extent. The pressure dependence of the diffusion coefficients, Di, is listed in Table 3. The decreasing tendency of the diffusion coefficients also suggests the structure formation in the higher-pressure regions.  To examine the pressure dependence of the transport properties, the diffusion coefficients of constituent molecules and ions are calculated, which are obtainable from the inclination of the mean square displacement (MSD) using Equation (11). The obtained MSD for 30 MP is shown in Figure 4. Although small oscillations still remain in the diffusion coefficient for Cl´, the inclinations of MSDs seem to converge until 5ˆ10´1 2 s (5 ps) and the slopes are kept for longer periods, to a certain extent. The pressure dependence of the diffusion coefficients, D i , is listed in Table 3. The decreasing tendency of the diffusion coefficients also suggests the structure formation in the higher-pressure regions.  The interference functions, Iij(q), which is also obtainable from the neutron diffraction experiment, can be calculated from the Fourier transformation of the pair distribution function gij(r) as expressed in Equation (10) [21]. The evaluated Iij(q)s from gij(r)s obtained by MD are shown in 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 275 K, 100 MP in Figure 3b. Although the peak heights of INa-Ow(q) and ICl-Ow(q) are similar to that of IC-Ow(q) at 283 K, 0.5 MP in Figure 3a, their heights decrease and their widths are broadened at 275 K, 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. To examine the pressure dependence of the transport properties, the diffusion coefficients of constituent molecules and ions are calculated, which are obtainable from the inclination of the mean square displacement (MSD) using Equation (11). The obtained MSD for 30 MP is shown in Figure 4. Although small oscillations still remain in the diffusion coefficient for Cl − , the inclinations of MSDs seem to converge until 5 × 10 s (5 ps) and the slopes are kept for longer periods, to a certain extent. The pressure dependence of the diffusion coefficients, Di, is listed in Table 3. The decreasing tendency of the diffusion coefficients also suggests the structure formation in the higher-pressure regions.   Next, we calculate the thermal conductivity to investigate the thermal properties of 1.1 mol % NaCl aqueous solution with saturated CO 2 . To the best of the author's knowledge, this is the first report on the thermal conductivity of CO 2 and NaCl aqueous solution system obtained by MD. The perturbation is applied to the system in the thermal equilibrium at time t = 0. According to the Green-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 be derived by NEMD, postulating that the contribution from the atom in the same molecule is omitted from the perturbation current [22]. The thermal conductivity obtained by NEMD under various pressures is shown in Figure 5a, alongside the experimental data of pure water and 1 mol % NaCl aqueous solution [23,24]. The results of NEMD for the saturated concentration of CO 2 in 1.1 mol % NaCl aqueous solution deserves special mention: the thermal conductivity decreases above 80 MP, which forms a striking contrast with the positive pressure dependence of other thermal conductivity data on solutions, in which CO 2 is not included. This anomaly of thermal conductivity also signifies the structure change of the CO 2 absorbed NaCl aqueous solution under high pressure.  Next, we calculate the thermal conductivity to investigate the thermal properties of 1.1 mol % NaCl aqueous solution with saturated CO2. To the best of the author's knowledge, this is the first report on the thermal conductivity of CO2 and NaCl aqueous solution system obtained by MD. The perturbation is applied to the system in the thermal equilibrium at time t = 0. According to the Green-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 be derived by NEMD, postulating that the contribution from the atom in the same molecule is omitted from the perturbation current [22]. The thermal conductivity obtained by NEMD under various pressures is shown in Figure 5a, alongside the experimental data of pure water and 1 mol % NaCl aqueous solution [23,24]. The results of NEMD for the saturated concentration of CO2 in 1.1 mol % NaCl aqueous solution deserves special mention: the thermal conductivity decreases above 80 MP, which forms a striking contrast with the positive pressure dependence of other thermal conductivity data on solutions, in which CO2 is not included. This anomaly of thermal conductivity also signifies 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 the specific heat at the constant pressure. In (a,b), the horizontal axes show pressure. The abbreviates "aq.", "sol.", "exp." stand for "aqueous", "solution", and "experiment", respectively. The dotted lines 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 shown in Figure 5b. The significant increase of Cp under higher pressure can be seen, although the experimental and MD data for pure water and seawater show a decreasing tendency of Cp against pressure. These results suggest the possibility of heat storage in the depths of the sea. The pressure dependence of the specific heat at the constant pressure. In (a,b), the horizontal axes show pressure. The abbreviates "aq.", "sol.", "exp." stand for "aqueous", "solution", and "experiment", respectively. The dotted lines 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, C p , obtained by MD using Equation (13) is 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 C p against pressure. These results suggest the possibility of heat storage in the depths of the sea.

The System Including Water Molecule, Na + , Cl´, H + , and HCO 3´I ons
As stated in the previous subsection, the thermodynamic properties of CO 2 and NaCl aqueous solution show the anomalous features under high pressure. These facts prompt us to a further study. As mentioned in the previous section, the CO 2 molecule is ionized to form HCO 3´i n the neutral pH. In this study, we wish to show the results of MD on seawater saturated with HCO 3´i ons as a more realistic model. MD is performed in 1.1 mol % NaCl aqueous solution with saturated HCO 3´f rom 0.44 to 7.97 mol %, corresponding to 5-1200 atm [25]. Then, for the calculation, MD cell contains 2500 TIP4P, 28 Na + , 28 Cl´, and 11 to 219 H + and HCO 3´.
The pair distribution function, g ij (r), and the integrated coordination number, n ij (r), obtained from MD results are shown in Figures 6 and 7 under pressures from 5 to 1000 atm. As seen in Figure 6a, g COw (r) between C(HCO 3´) and O(water) has a pronounced peak at around 4 Å. The coordination number of water around a HCO 3´i s estimated to be 17, which agrees well with those in the literature [26]. The sharp first peaks of g NaOw (r) are found around 2.3 Å in Figure 6b, which shows the close distance between cations and water molecules. Figure 6a,b correspond to those of C(CO 2 )-O(water) and Na-O(water) in Figure 2a,b. To confirm the pressure dependence of the coordination number, which was seen in the Section 2.1, we calculate the coordination number to the first minimum of g ij (r). The obtained results are listed in Table 4. The negative pressure dependence of the coordination number is also observed, which is also the collaborating evidence of the structure formation around HCO 3´i ons. In Figure 7a,b, g CNa (r), C(HCO 3´) -Na + , g CH (r), and C(HCO 3´) -H + have pronounced two split peaks from 2.5 to 4.0 Å, which may be attributed to the asymmetric form of HCO 3´. As stated in the previous subsection, the thermodynamic properties of CO2 and NaCl aqueous solution show the anomalous features under high pressure. These facts prompt us to a further study. As mentioned in the previous section, the CO2 molecule is ionized to form HCO3 − in the neutral pH. In this study, we wish to show the results of MD on seawater saturated with HCO3 − ions as a more realistic model. MD is performed in 1.1 mol % NaCl aqueous solution with saturated HCO3 − from 0.44 to 7.97 mol %, corresponding to 5-1200 atm [25]. Then, for the calculation, MD cell contains 2500 TIP4P, 28 Na + , 28 Cl − , and 11 to 219 H + and HCO3 − .
The pair distribution function, gij(r), and the integrated coordination number, nij(r), obtained from MD results are shown in Figures 6 and 7 under pressures from 5 to 1000 atm. As seen in Figure  6a, gCOw(r) between C(HCO3 − ) and O(water) has a pronounced peak at around 4 Å. The coordination number of water around a HCO3 − is estimated to be 17, which agrees well with those in the literature [26]. The sharp first peaks of gNaOw(r) are found around 2.3 Å in Figure 6b, which shows the close distance between cations and water molecules. Figure 6a,b correspond to those of C(CO2)-O(water) and Na-O(water) in Figure 2a,b. To confirm the pressure dependence of the coordination number, which was seen in the Section 2.1, we calculate the coordination number to the first minimum of gij(r). The obtained results are listed in Table 4. The negative pressure dependence of the coordination number is also observed, which is also the collaborating evidence of the structure formation around HCO3 − ions. In Figure 7a     In Figure 8a, the diffusion coefficients of water molecule, HCO3 − and Na + , or DO(water), DC(HCO3-), and DNa, obtained from MSD defined by Equation (11), are plotted against pressure. The obtained values agree well those in the literature [27]. As seen in Figure 8a, all Di s decrease as the pressure increases. It is noteworthy that DC(HCO3-), and DNa show similar pressure dependence. These features of gij(r) s and Di s suggest that the complex [HCO3·(H2O)n] − is expected to be formed in the solution. Then, the clusters {Na + ·[HCO3·(H2O)n] − } and {H + ·[HCO3·(H2O)n] − } should be compounded to hold the local charge neutrality. Similar structures have also been found in the aqueous solutions. According to the ab initio MD study of Na + in aqueous solution, the n-coordinate hydration structures, such as Na(H2O)n + , have been found [28]. In an aqueous solution of CaCO3, Ca(HCO3)2(H2O)4 and Ca(HCO3)3(H2O)2 − are predicted to be stable [29].  In Figure 8a, the diffusion coefficients of water molecule, HCO 3´a nd Na + , or D O(water) , D C(HCO3-) , and D Na , obtained from MSD defined by Equation (11), are plotted against pressure. The obtained values agree well those in the literature [27]. As seen in Figure 8a, all D i s decrease as the pressure increases. It is noteworthy that D C(HCO3-) , and D Na show similar pressure dependence. These features of g ij (r) s and D i s suggest that the complex [   In Figure 8a, the diffusion coefficients of water molecule, HCO3 − and Na + , or DO(water), DC(HCO3-), and DNa, obtained from MSD defined by Equation (11), are plotted against pressure. The obtained values agree well those in the literature [27]. As seen in Figure 8a, all Di s decrease as the pressure increases. It is noteworthy that DC(HCO3-), and DNa show similar pressure dependence. These features of gij(r) s and Di s suggest that the complex [HCO3·(H2O)n] − is expected to be formed in the solution. Then, the clusters {Na + ·[HCO3·(H2O)n] − } and {H + ·[HCO3·(H2O)n] − } should be compounded to hold the local charge neutrality. Similar structures have also been found in the aqueous solutions. According to the ab initio MD study of Na + in aqueous solution, the n-coordinate hydration structures, such as Na(H2O)n + , have been found [28]. In an aqueous solution of CaCO3, Ca(HCO3)2(H2O)4 and Ca(HCO3)3(H2O)2 − are predicted to be stable [29].  The frequency dependent diffusion coefficient, D i pωq, can be derived from the velocity auto-correlation function (VAF), or xv i ptq¨v i p0qy using Equation (12). The obtained D CpHCO3 q pωq for HCO 3´u nder various pressures are observed in the THz or the infrared region as shown in Figure 8b. The peak of D CpHCO3 q pωq around 300 1/cm under low pressure is very close to the frequency of water caused by the translational cage effect. The peak position slightly sifts to the higher frequency around 500 1/cm, and the small hump around 1000 1/cm can be observed, which are comparable to the CO 2 -H 2 O intermolecular vibrational frequency [30].
In order to evaluate the lifetime of the complex [HCO 3¨( H 2 O) n ]´, we calculate the rotational correlation function, C 2 (t), of HCO 3´i on defined by Equation (19). The C 2 (t) is thought to be affected by the relaxation of the interaction between HCO 3´a nd the surrounding water molecules. The logarithm plot of the obtained C 2 (t)s in various pressure are shown in Figure 9. The lifetime can be evaluated from the inclination of the linear part of lnC 2 (t). The graph of lnC 2 (t) of HCO 3´a t 5 atm is extremely similar to that of water molecule. The oscillatory behavior at 3-4 ps may be interpreted as the "free-rotor" motion, which is observed in the dilute phase [15]. The estimated lifetime at 5 atm is 1.6 ps; on the other hand, those of at 100-1000 atm is 6.7˘1.3 ps, which is comparable to the relaxation time of H-bond in the aqueous carbonate solution [31]. The slight positive pressure dependence of the relaxation time also suggests the structure formation in the higher-pressure region. The frequency dependent diffusion coefficient, D (ω), can be derived from the velocity autocorrelation function (VAF), or 〈 ( ) • (0)〉 using Equation (12). The obtained D ( ) (ω) for HCO3 − under various pressures are observed in the THz or the infrared region as shown in Figure 8b. The peak of D ( ) (ω) around 300 1/cm under low pressure is very close to the frequency of water caused by the translational cage effect. The peak position slightly sifts to the higher frequency around 500 1/cm, and the small hump around 1000 1/cm can be observed, which are comparable to the CO2-H2O intermolecular vibrational frequency [30].
In order to evaluate the lifetime of the complex [HCO3·(H2O)n] − , we calculate the rotational correlation function, C2(t), of HCO3 − ion defined by Equation (19). The C2(t) is thought to be affected by the relaxation of the interaction between HCO3 − and the surrounding water molecules. The logarithm plot of the obtained C2(t)s in various pressure are shown in Figure 9. The lifetime can be evaluated from the inclination of the linear part of lnC2(t). The graph of lnC2(t) of HCO3 − at 5 atm is extremely similar to that of water molecule. The oscillatory behavior at 3-4 ps may be interpreted as the "free-rotor" motion, which is observed in the dilute phase [15]. The estimated lifetime at 5 atm is 1.6 ps; on the other hand, those of at 100-1000 atm is 6.7 ± 1.3 ps, which is comparable to the relaxation time of H-bond in the aqueous carbonate solution [31]. The slight positive pressure dependence of the relaxation time also suggests the structure formation in the higher-pressure region. For the next stage, according to the G-K formula, the shear viscosity is calculated using Equation (17). The shear viscosities obtained by MD are shown in Figure 10a with the experimental values for pure water, 0.6 M NaCl aqueous solution, and 0.6 M NaCl and 0.913 M CO2 aqueous solution in the literature [32,33]. The pressure dependence of the present result (0.6 M NaCl with saturated HCO3 − ) is positive, whereas those of experimental values are negative. This fact suggests that the interaction between HCO3 − and water, and/or other constituents increases as the pressure increases. A significant increasing tendency of viscosity with increasing mole fraction of dissolved CO2 has also 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 diffusion coefficient obtained by MD using the Stokes-Einstein (S-E) relation for a spherical particle, which is expressed as For the next stage, according to the G-K formula, the shear viscosity is calculated using Equation (17). The shear viscosities obtained by MD are shown in Figure 10a with the experimental values for pure water, 0.6 M NaCl aqueous solution, and 0.6 M NaCl and 0.913 M CO 2 aqueous solution in the literature [32,33]. The pressure dependence of the present result (0.6 M NaCl with saturated HCO 3´) is positive, whereas those of experimental values are negative. This fact suggests that the interaction between HCO 3´a nd water, and/or other constituents increases as the pressure increases. A significant increasing tendency of viscosity with increasing mole fraction of dissolved CO 2 has also been observed in the viscosity measurement of CO 2 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 diffusion coefficient obtained by MD using the Stokes-Einstein (S-E) relation for a spherical particle, which is expressed as where the parameters ξ and η stand for the friction constant and the shear viscosity, respectively. If the shear viscosity η is determined at a certain CO 2 concentration c 0 , then η(c) at any concentration c could be estimated using the following equation [35]: which is known as the Walden's rule. The calculated η(c), from Equation (3), is also plotted in Figure 10a, which agrees with the MD results to a certain extent. where the parameters ξ and η stand for the friction constant and the shear viscosity, respectively. If the shear viscosity η is determined at a certain CO2 concentration c0, then η(c) at any concentration c could be estimated using the following equation [35]: which is known as the Walden's rule. The calculated η(c), from Equation (3), is also plotted in Figure  10a, which agrees with the MD results to a certain extent. As mentioned in the previous subsection, we have calculated the thermal conductivity of 1.1 mol % NaCl aqueous solution with saturated CO2 by NEMD method [7,16]. In this study, we adopt the same method using the saturated HCO3 − ion in 1.1 mol % NaCl aqueous solution. As will be described in Section 3, the thermal conductivity λ is expressed as Equations (14)- (16). The obtained results of the thermal conductivity by NEMD are shown in Figure 10b. The experimental data of pure water and seawater are also shown in Figure 10b [23,36]. The negative pressure dependence of thermal conductivity is clearly seen in the MD result; on the other hand, those of the experimental data are positive.
As stated already, some anomalous results have been obtained by MD in the transport and thermal properties of 1.1 mol % NaCl aqueous solution saturated with HCO3 − . The experimental thermal conductivity data of electrolyte aqueous solutions show positive pressure dependence, and negative concentration dependence of electrolyte [37]. These phenomena have been explained to some extent by the extension of the additivity of the thermal conductivity by considering the interaction between components [37,38]. The results in this study might be influenced by the above-mentioned contradictory effects to the thermal conductivity, pressure and concentration. In addition, the results are also supposed to be attributed to the complex and/or the cluster formation in the solution.

The Methane and NaCl Aqueous Solution
Next, we will show the MD result of the methane and NaCl aqueous solution. The fundamental procedure of MD is the same as described in the previous subsections. The water (TIP4P) and the methane are treated as rigid body molecules. The concentration of NaCl is adjusted to be the same as that of seawater, 1.1 mol %. The number of CH4 in the MD cell is determined using the solubility data of CH4 in seawater [39]. The numbers of particles used in MD are listed in the Table 5. As mentioned in the previous subsection, we have calculated the thermal conductivity of 1.1 mol % NaCl aqueous solution with saturated CO 2 by NEMD method [7,16]. In this study, we adopt the same method using the saturated HCO 3´i on in 1.1 mol % NaCl aqueous solution. As will be described in Section 3, the thermal conductivity λ is expressed as Equations (14)- (16). The obtained results of the thermal conductivity by NEMD are shown in Figure 10b. The experimental data of pure water and seawater are also shown in Figure 10b [23,36]. The negative pressure dependence of thermal conductivity is clearly seen in the MD result; on the other hand, those of the experimental data are positive.
As stated already, some anomalous results have been obtained by MD in the transport and thermal properties of 1.1 mol % NaCl aqueous solution saturated with HCO 3´. The experimental thermal conductivity data of electrolyte aqueous solutions show positive pressure dependence, and negative concentration dependence of electrolyte [37]. These phenomena have been explained to some extent by the extension of the additivity of the thermal conductivity by considering the interaction between components [37,38]. The results in this study might be influenced by the above-mentioned contradictory effects to the thermal conductivity, pressure and concentration. In addition, the results are also supposed to be attributed to the complex and/or the cluster formation in the solution.

The Methane and NaCl Aqueous Solution
Next, we will show the MD result of the methane and NaCl aqueous solution. The fundamental procedure of MD is the same as described in the previous subsections. The water (TIP4P) and the methane are treated as rigid body molecules. The concentration of NaCl is adjusted to be the same as that of seawater, 1.1 mol %. The number of CH 4 in the MD cell is determined using the solubility data of CH 4 in seawater [39]. The numbers of particles used in MD are listed in the Table 5. From the MD results, the pair distribution functions, g ij (r)s, and the integrated coordination numbers, n ij (r)s, have been obtained for 10 to 100 MP, which are shown in Figure 11a-d. Although the pressure dependence of g ij (r) is not large, the slight change of the first peak height for CH 4 -CH 4 , and the depth for the first minimum for H 2 O-H 2 O can be observed, which may correspond to the cage formation in the solution. The water coordination number, n ij (r), of the first hydration shell around the solute is calculated to the first minimum of g ij (r) using Equation (9). The obtained water coordination numbers calculated under pressures of 10-100 MP are listed in Table 6. The slight decreasing tendency of water molecules around CH 4 has been detected, which might also be attributed to the cluster formation around CH 4 molecules.  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 the pressure dependence of gij(r) is not large, the slight change of the first peak height for CH4-CH4, and the depth for the first minimum for H2O-H2O can be observed, which may correspond to the cage formation in the solution. The water coordination number, nij(r), of the first hydration shell around the solute is calculated to the first minimum of gij(r) using Equation (9). The obtained water coordination numbers calculated under pressures of 10-100 MP are listed in Table 6. The slight decreasing tendency of water molecules around CH4 has been detected, which might also be attributed to the cluster formation around CH4 molecules. (c) (d) Figure 11. (a) gOwOw(r) and nOwOw(r); (b) gNaOw(r) and nNaOw(r); (c) gClOw(r) and nClOw(r) ; and (d) gCH4Ow(r) and nCH4Ow(r) under pressures of 10-100 MP. Figure 11. (a) g OwOw (r) and n OwOw (r); (b) g NaOw (r) and n NaOw (r); (c) g ClOw (r) and n ClOw (r) ; and (d) g CH4Ow (r) and n CH4Ow (r) under pressures of 10-100 MP. Finally, the pressure dependence of thermal conductivity of methane and NaCl aqueous solution is obtained by NEMD using Equations (14)- (16). The obtained results are shown in Figure 12 alongside the experimental data and MD data for NaCl aqueous solution and pure water. The negative pressure dependence of thermal conductivity in higher pressure is also observed. This result might be attributed to the structure change or the clathrate formation around the CH 4 molecule in the high-pressure region, which is consistent with the discussion regarding the decreasing of coordination number in the solution.  Finally, the pressure dependence of thermal conductivity of methane and NaCl aqueous solution is obtained by NEMD using Equations (14)- (16). The obtained results are shown in Figure 12 alongside the experimental data and MD data for NaCl aqueous solution and pure water. The negative pressure dependence of thermal conductivity in higher pressure is also observed. This result might be attributed to the structure change or the clathrate formation around the CH4 molecule in the high-pressure region, which is consistent with the discussion regarding the decreasing of coordination number in the solution.

Procedure
The essential numerical procedure of this MD simulation study is the same as our previous works on aqueous solutions [7,8,12,40]; however, the fundamental part of the procedure is described as follows for the reader's benefit.
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 In the equations hereafter, i and j stand for the constituent atoms; zi is the charge for the constituent species i; and e is the elementary charge. The used parameters are listed in Table 7.

Procedure
The essential numerical procedure of this MD simulation study is the same as our previous works on aqueous solutions [7,8,12,40]; however, the fundamental part of the procedure is described as follows for the reader's benefit.
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 In the equations hereafter, i and j stand for the constituent atoms; z i is the charge for the constituent species i; and e is the elementary charge. The used parameters are listed in Table 7. The interactions between Na + and Cl´, TIP4P-Na + , and TIP4P-Cl´are expressed as [41], The used parameters are listed in Table 8. The interactions between other solute ions and water molecules are also expressed in the charged L-J type potentials [42,43] as φ ij prq " z i z j e 2 r`A r 12´B r 6 (6) The parameters of the interactions between CO 2 molecule, HCO 3´i on, Na + , Cl´and water molecule, taken from the literature, are listed in Tables 9-12. Table 9. The potential parameters between CO 2 -CO 2 , CO 2 -Na + , CO 2 -Cl´, and CO 2 -TIP4P water.   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´1q 3 ă pthe total number of water molecules and solute ions and moleculesq ă n 3 Then, the relaxation of the first configuration using the Monte Carlo method is executed by the following procedure: (a) One ion or molecule is randomly selected. (b) 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. (c) 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 ∆r max . A random degree is applied for "∆Ω" and "∆θ". (d) The above procedures (a)-(c) are repeated for the defined number of steps until the ensemble average is calculated. (e) The potential energy difference ∆φ between the previous configuration {ξ} and the new configuration {ξ} 1 is calculated. (f) 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][19][20]. The very short time step is adopted so as to detect the fast movement of H + and water molecules.
The pair distribution function, g ij (r), can be firstly obtained from a time series data of coordinates of ions as [15] g ij prq "`V{N i N j˘N i ÿ k c kjˆr´∆ r 2 , r`∆ r 2˙{´4 πr 2 ∆r¯ (8) where N i and N j are the numbers of the ion species i and j, respectively. V is the volume of the cell. c ik 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, n ij (r), is defined as The interference functions for neutron I ij (q) are obtained from g ij (r), which is expressed as [21] I ij pqq " c i b i c j b j {˜ÿ k c k b k¸2 ż 8 0 4πr 2 ρ 0´gij prq´1¯s in prqq qr dr (10) where c i is the atomic fraction of the i-type atoms; ρ 0 is the average number density; and b i is the neutron scattering amplitude of the i-type atom. The diffusion coefficient for i-type atom, D i , can be obtained from MSD, which is defined as VAF is calculated to examine the dynamical and transport properties, which is expressed as xv i ptq¨v i p0qy. The frequency dependent diffusion coefficient, D i pωq, can be obtained from VAF as D i pωq " 1 3 ż 8 0 xv i ptq¨v i p0qy cosωt dt (12) where ω = 2π f ; f is the frequency, and D i (ω = 0) = D i . The thermodynamic properties are also important for the study of solutions. The specific heat at the constant pressure is expressed as 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], In Equation (14), J x (τ) stands for the x component of the perturbation current; V the volume of the system; k B the Boltzmann constant; and T the temperature of the system. The applied perturbation F ext is related to the average of the perturbation current <J x (τ)> t , xJ x y t " F ext k B T ż t 0 xJ x pτq J x p0qy dτ (15) Then, the relation between the thermal conductivity and the perturbation is expressed as [15] λ " 1 VF ext T lim tÑ8 xJ x y t (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], η " 1 kTV The reorientation of linear molecules is expressed as the time-correlation function defined as [15] C l ptq " xP l ru i ptq¨u i sy (18) where u i (t) is unit vector parallel to the principal axis of the molecule i. P l (x) is the Legendre polynomial. The rotational correlation function is the case of l = 2, which is expressed as