Diffusion in Various Carbonated Beverages - ACS Publications

Jan 9, 2018 - burst, which all depend on the diffusional behavior of gases.4,5 ...... Simulations through Multi-Level Parallelism from Laptops to Supe...
0 downloads 0 Views 817KB Size
Article pubs.acs.org/JPCB

Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

CO2 Diffusion in Various Carbonated Beverages: A Molecular Dynamics Study Ji Lv,†,‡ Kaixin Ren,† and Yakun Chen*,† †

Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun, Jilin Province 130023, People’s Republic of China ‡ Department of Chemistry, University of Minnesota, Minneapolis, Minnesota 55455, United States S Supporting Information *

ABSTRACT: Carbonated beverages are widely enjoyed in spare time, yet there remain many physical and chemical processes clouded at the molecular level. In this report, we employ molecular dynamics simulations to estimate the diffusion coefficients of CO2 and the molecular origin of its variations in three model systems with characteristic features of champagnes, sugar-based cola drinks, and club sodas. The computed diffusion coefficients of CO2 are in good agreement with experimental data. Analyses of hydrogen bonding and the solvent’s structural and dynamic properties reveal that the change of CO2 diffusion coefficient is closely associated with the diffusional behavior of the solvent water itself, as a result of changes in the number and strength of hydrogen bonding interactions among the species and solvent.



INTRODUCTION An essential component of carbonated beverages (e.g., champagnes, beers, sodas, and sparkling water) is the dissolved CO2 in the liquid mixtures. The size and stability of the bubbles play an important role in determining the degree of creaminess and smoothness, and therefore the level of enjoyment for consumption.1−3 Both properties are closely related to the life cycle of bubbles, including nucleation, growth, ripening, and burst, which all depend on the diffusional behavior of gases.4,5 Therefore, it is of interest to determine diffusion coefficients of various gases dissolved in aqueous solutions. This, however, is a challenging task experimentally for water−beverage mixtures because of the low solubility of the gaseous molecules.6 For example, in a study by Sharma and Adhikari, the uncertainty of experimental results on the diffusivity of N2 in aqueous solution was found to be as large as 40%.7 Recently, the diffusion coefficients of CO2 in mixed aqueous solutions have been computed by several researchers, although the origin of the effects of added ingredients on the diffusion coefficient is still not fully understood.8 Carbonated beverages are among the most popular drinks in the world. There are many interesting scientific and technological questions in their production, which have attracted extensive studies. Liger-Belair and his co-workers3,9−11 investigated the interplay between gaseous and dissolved CO2 in champagnes by theoretical and experimental methods. They paid attention to the nucleation, rise, and burst of CO2 bubbles and the diffusion of dissolved CO2 in champagnes. Nose et al.12−14 investigated the hydrogen bonding network in alcoholic beverages and the mechanism that it is affected by different solutes employing NMR experiments. It was found that different solutes (salts, acids, etc.) can either break or strengthen the network structure of water depending on the charge density of ions, acid strength (pKa), etc. Zhang and Xu5 © XXXX American Chemical Society

investigated the bubble growth in beer and champagne and noted that bubble growth is very slow in Guinness beer, resulting in smaller, more stable, and sinking bubbles. However, these studies focused on alcoholic drinks and paid less attention to soft drinks. Soft drinks include a variety of additive components, such as sugar and salt. These additive components greatly affect the properties of solution.15,16 In this work, we perform a molecular dynamics simulation study of carbon dioxide diffusion in mixed aqueous solutions as models for three types of carbonated beverages, including champagnes, cola drinks, and club sodas. The effects of ethanol, sucrose, or salt ions on the diffusion of CO2 in aqueous mixtures are examined, and analyses of the simulation results reveal factors that influence the gas diffusion in the mixed aqueous environments.



METHOD AND MODELS System Preparation. We constructed three systems in the present study to model carbonated beverages, with characteristic features resembling champagnes, cola drinks, and club sodas. As a control and validation of the computational procedures and force fields used, the diffusion coefficient of CO2 in water was also examined and compared with previous computational and experimental studies.17 In this control case,17 10 molecules of CO2 were introduced in a cubic box consisting of 104 water molecules. For champagne wines, we used a system identical to that in the work of Perret and his coworkers,9,17 which is a mixture of 50 CO2 molecules, 440 ethanol molecules, and 104 water molecules, corresponding to 12−13% vol/vol alcohol content and 12 g/L CO2. To model Received: October 23, 2017 Revised: January 6, 2018 Published: January 9, 2018 A

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the number of hydrogen bonds per molecule because different systems contain different organic and inorganic additives. The hydrogen bond lifetime τHB, the time interval a hydrogen bond breaks, provides a measure of the dynamic flexibility and strength of a hydrogen bond.21 The latter is a key factor that affects the diffusion rate of the solvent, as well as that of the solute CO2 in various systems. Simulation Details. All molecular dynamics simulations were performed using GROMACS v5.02.22 The CHARMM27 force field was used for carbon dioxide, sucrose, ethanol molecules, and bicarbonate ion.23 The five-point charge TIP5P model for water was used, which has been shown to yield good results for the diffusion of aqueous solution in previous studies, and the CO2 diffusion coefficient in dilute aqueous solution in comparison with experimental results.17,24 A shifting function was used to feather Lennard-Jones interactions to zero between 10 and 12 Å, whereas long-range electrostatic interactions were modeled by the particle mesh Ewald (PME) approach.25 A twostep equilibration procedure was carried out to equilibrate each system, first by 200 ps molecular dynamics simulations employing the NVT ensemble, followed by 1 ns equilibration using the NPT ensemble. The diffusion coefficients were determined through a total of 1 ns molecular dynamics simulations. Periodic boundary conditions were used in all calculations at 293 K and 1 bar. The Nosé−Hoover thermostat and Parinello−Rahman algorithm were used with a time constant of 0.3 ps when necessary.26,27 The leapfrog algorithm was used with an integration time step of 1 fs for water and 0.5 fs for other systems. The trajectories were recorded in every 1 ps for water and 0.5 ps for other systems.

cola drinks in which a typical 355 mL can contains about 39 g of sucrose and 2.2 g of CO2, a mixture of 28 CO2 molecules, 64 sucrose molecules, and 104 water molecules (0.36 M concentration) was constructed (for the ingredients of a typical cola drink, see Table S1 for details). Finally, we used a system containing roughly a 2:1 mixture of carbon dioxide and sodium bicarbonate (21 NaHCO3, or 0.12 M solution, and 10 CO2 molecules) in 104 water molecules to mimic club sodas. The salt NaHCO3 was used to neutralize the acidity resulting from flavor additives that are acidic in these drinks (for the ingredients of a typical club soda, see Table S2 for details). Of course, there are numerous other minor components, which are ignored in the present study but they could affect the diffusion property of CO2 in the real systems. Nevertheless, the present simplifications are expected to introduce relatively small errors in the calculated diffusion coefficients, capturing the dominant factors on CO2 diffusion as previously shown in similar studies.9 The initial coordinates for all systems were built with PACKMOL18 and were minimized by the steepest descent method to result in starting configurations for molecular dynamics simulations. Diffusion Coefficient. Two approaches are widely used to estimate the diffusion coefficient of a solute in a condensedphase medium. The first is based on the Einstein relationship D = lim

t →∞

⟨|r(t ) − r(0)2 |⟩ 6t

(1)

where ⟨|r(t) − r(0)2|⟩ is the mean-square displacement (MSD) of a particle s at time t with respect to time 0. Alternatively, the diffusion coefficient can be determined from the velocity autocorrelation function (VACF) by D=

1 3

∫0



RESULTS AND DISCUSSION Diffusion Coefficient. The MSDs of CO2 in different systems (Figure 1) show good linear dependence with time,



dt ⟨v(t )v(0)⟩

(2)

where ⟨v(t)v(0)⟩ is the VACF. Both methods are mathematically equivalent for simple fluids where the medium is isotropic and homogeneous.19 Bonhommeau and his co-workers9,17 have demonstrated that eq 1 is appropriate for modeling CO2 diffusion in champagnes. In this study, we also found that it is sufficient to use these approaches to investigate CO2 diffusion in dilute aqueous sugared or ionic solution. In addition, we have examined the system-size dependence of the calculated diffusion coefficients in water by doubling the number of water molecules to 2 × 104 water molecules, and reducing it to 5 × 103 (Table S3). It was found that the computed values are converged with little variations using the three box sizes. This finding is consistent with that found in a previous study of N2 in water using both eqs 1 and 2, showing that the diffusion coefficients of infinitely dilute solute do not vary significantly with the change of the system size.7 Statistical Uncertainty. We used the multiorigin algorithm to calculate the diffusion coefficient.17 Then, the error estimation was used to assess whether the simulation was long enough for the ergodic hypothesis to be valid. According to the central limit theorem, if D is normally distributed, then we think that the sample is sufficiently large. Hydrogen Bond Analysis. Hydrogen bond analysis was performed following a geometric definition,20 in which a hydrogen bond is defined for an interaction where the donor− acceptor distance (OD−OA) is less than 3.5 Å and the angle OD−H−OA is greater than 145°. We use the total number of hydrogen bonds in a given system for comparison instead of

Figure 1. Computed mean-square displacement (MSD) of the center of mass of CO2 in various simulation systems at 293 K, corresponding to aqueous solutions of CO2 (blue), 12% v/v ethanol−water mixture corresponding to contents of typical champagnes (red), 0.36 M sucrose solution as in a typical can of cola drink (cyan), and 0.12 M sodium bicarbonate solution as in club sodas (purple).

except toward the end where the averages have greater uncertainties, demonstrating a Fickian behavior. According to the central limit theorem, the dynamics of particles following classic random walks must be Fickian and Gaussian.28 However, being Fickian does not necessarily result from a Gaussian or non-Gaussian distribution of particle displacement, thus not guaranteeing the use of the MSD to estimate the diffusion coefficient.28−31 However, several molecules present in the current study, particularly the relatively large sucrose, might lead to a non-Gaussian diffusion of CO2 in the aqueous B

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Overall, the diffusion coefficient of CO2 decreases significantly upon addition of the other additives. For champagnes, our result (1.55 × 10−9 m2/s) is close to the recent theoretical value obtained by Bonhommeau et al.9 (1.64 × 10−9 m2/s) and their previous study17 (1.61 × 10−9 m2/s). However, it is slightly overestimated, in comparison with the recent experimental data obtained for the champagnes9 (1.29 × 10−9 m2/s). Because champagnes contain some other compounds, such as sugars, protein, etc., our data is close to the experimental data obtained for the hydroalcoholic samples9 (1.36 × 10−9 m2/s). For cola drinks and club sodas, to the best of our search, we did not find enough theoretical or experimental data other than those listed in Table 2. Although

solutions. To investigate the diffusion features of CO2, we calculated the van Hove correlation function Gs(r, t).32 In each case, Gs(r, t) can be satisfactorily fitted by a Gaussian function (see Figure S1 for details), strongly indicating that it is adequate to use eq 1 to estimate the diffusion coefficients in this study. The distributions of the diffusion coefficient of CO2 for all four systems on a series of time intervals are reported in Figure S2. It is noted that, despite that the simulation time is only 1 ns long for each system, it is already sufficiently long, as the distributions follow the normal distribution quite well. Therefore, the CO2 diffusion coefficients obtained in the current study should be considered converged. Listed in Table 1 are the computed diffusion coefficients of CO2 in water and in various carbonated beverage models along

Table 2. Computed Viscosity (10−3 Pa·s), Hydrogen Bond Lifetime (ps), and Average Number of Hydrogen Bonds in Different Carbonated Systems

Table 1. Computed and Experimental Diffusion Coefficients of CO2 (×10−9 m2/s) in H2O and in Different Carbonated Systems D

system water

12% v/v EtOH

cola model club sodas a c

temperature (K)

(×10−9 m2/s) of CO2

method

293 293 293 293 293 293 298 293 293 293 293 293 293 293 293 293 293 293 298 (1 MPa)

2.33 ± 0.16 2.36 ± 0.09 2.1 ± 0.3 2.1 ± 0.08 1.67 1.85a 1.91 1.55 ± 0.11 1.64 ± 0.11 1.32 ± 0.11 1.61 ± 0.1 1.02 ± 0.11 1.36 ± 0.01b 1.29 ± 0.00c 1.41c 1.56 ± 0.09 1.35 1.61 ± 0.12 1.34

TIP5P + MSD fit TIP5P + MSD fit17 SPC/E + MSD fit53 SPC/E + VACF53 experiment33 experiment8 experiment54 TIP5P + MSD fit TIP5P + MSD fit9 SPC/E + MSD fit9 TIP5P + MSD fit17 SPC/E + MSD fit17 experiment9 experiment9 experiment8 TIP5P + MSD fit experiment8 TIP5P + MSD fit experiment55

Sparkling water (Natural Springs). Champagne samples.

b

system

all hydrogen bonds

viscositya

hydrogen bond lifetimeb

water 12% v/v EtOH cola model club sodas

18112 18453 18600 18395

0.835 1.092 1.157 1.186

8.52 10.95 9.75 12.21

a

The theoretical viscosities calculated by evaluating the transverse current autocorrelation functions (TCAFs),56 for which the gmx tcaf GROMACS tool was used. bHere, we calculated the hydrogen bond lifetime (τHB) on the basis of the hydrogen bond time correlation function (HBTCF),21 for which the gmx hbond GROMACS tool was used.

the experimental data for club soda was obtained at 1 MPa, some experimental results35,36 and our extra MD simulation of a system containing CO2, H2O, and NaHCO3 at 1 MP (1.59 × 10−9 m2/s) suggested that the CO2 diffusion coefficient is almost independent of the pressure when the pressure is in the range between 1 and 40 bar. Interestingly, the diffusion coefficients of CO2 in the three carbonated beverages are very close to each other, either from experimental measurements or from MD simulations, whereas the molecular features of the additive solutes are quite different: both ethanol and sucrose are neutral molecules and possess hydroxyl groups, though with different volumes and different numbers of hydroxyl groups; on the other hand, NaHCO3 is hydrolyzed into Na+ and HCO3− ions. This drives us to carry out further analyses to account for such a phenomenon. Viscosity. The diffusion coefficient of CO2 is related to the solvent’s viscosity, which is in turn affected by the number and strength of hydrogen bonds in an aqueous solution. Consequently, the change in hydrogen bonding interactions by additives such as alcohol, sugar, and inorganic salt is often used to characterize the variations of CO2 diffusion in carbonated beverage models.17 Consistent with previous findings,9,37 we found that the presence of a greater number of hydrogen bonds leads to higher viscosity (Table 2) for alcohol and sugar solutions. However, this trend is not followed in club sodas with the NaHCO3 salt additive. In the latter case, the computed liquid viscosity depends not only on the number of hydrogen bonds but also on the strength of the hydrogen bond. We also calculated the hydrogen bond lifetime τHB; it is one of the key factors at the molecular level affecting the diffusion speed of the water, and also the diffusion speed of CO2 in various systems. On the basis of the data in Table 2, one notices that, in some

Hydroalcoholic samples.

with the experimental data. Employing the TIP5P model for water, we obtained a diffusion coefficient of 2.32 × 10−9 m2/s for CO2 in aqueous solution, which is in good agreement with the computed value of 2.36 × 10−9 m2/s from the work of Perret et al.17 For comparison, it was found experimentally that CO2 diffuses in the range (1.60−1.77) × 10−9 m2/s.33 Although the diffusion coefficient for CO2 is slightly overestimated, a number of factors may contribute to the computational results, possibly including the TIP5P water model itself and the simulation protocols. An illustration of the accuracy of the TIP5P model is that the self-diffusion coefficient of the TIP5P34 model for water (1.74 × 10−9 m2/s) is in good agreement with experimental data (1.77 × 10−9 m2/s) at 288 K, although the discrepancy between computational results and the experimental values increases at higher temperatures. Since the main emphasis of the present study is to understand the trends and factors that affect CO2 diffusion in different carbonated solutions, the force fields and computational procedures used here are reasonable in the following analysis. C

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

aqueous solution, but it remains unchanged when sodium bicarbonate is added. Remarkably, in all three beverage models, none of the added components are found in direct contact with CO2, suggesting that the change in CO2 diffusion is due to indirect effects of alteration in liquid structure. Water Tetrahedral Order q. Several geometric order parameters have been defined to identify the local structural properties of water around the solutes.40 For example,41 in a study of concentrated aqueous NaOH solution, analysis of local order parameters revealed that the increased fluid viscosity in that system is associated with a greater number of water molecules with a trigonal bipyramidal arrangement around Na+. Here, we employ a tetrahedral order parameter, q, one of the most widely used order parameters,15,42,43 to characterize the hydrogen bonding structure about the solvent water molecules.44 Here, the tetrahedral order parameter is defined on the basis of the four nearest oxygen atoms surrounding a central water molecule

cases, the importance of hydrogen bond strength may outweigh the number. For example, the largest viscosity is found in the club sodas in which the strength of hydrogen bonding interactions is exceptionally strong, as reflected by the relatively large hydrogen bond lifetime τHB. In addition, the great strength of hydrogen bond in champagnes is also in agreement with recent results on the energetics of hydrophobic hydration of polypeptides.38 Clearly, the computed viscosity is not determined by a single factor, but the balance between the number and the strength of hydrogen bonding interactions in the system is critical. Hydrogen Bond Analysis. Figure 2 shows the average number of hydrogen bonds in the entire simulation system

q=1−

3 8

3

4

2 ⎛ 1⎞ ⎜cos θ + ⎟ jk ⎝ 3⎠ k=j+1

∑ ∑ j=1

(3)

where θjk is the angle between the central oxygen atom and its neighboring oxygen atoms j and k. In the present analysis of mixed solutions, q is determined by using the four nearest oxygen atoms, irrespective of if they are from water or a solute molecule. For a perfect tetrahedral configuration, eq 3 yields a q value of 1, whereas its average is zero for a completely random distribution of four nearest neighbors. It is known that the distribution of the tetrahedral order is typically characterized by two peaks, a high-q peak (ice-like, i.e., a relative structured configuration) and a low-q peak (fluid-like configurations). Figure 3 shows a bimodal distribution for each of the four systems studied at 293 K. In the aqueous solution, the positions

Figure 2. Computed average number of hydrogen bonds in an aqueous solution containing dissolved CO2 (blue), 12% v/v ethanol− water mixture corresponding to typical champagnes (red), 0.36 M sucrose solution as in a typical can of cola drink (cyan), and 0.12 M sodium bicarbonate solution as in club sodas (purple).

under conditions corresponding to the major components of various carbonated beverages. In each solution, understandably, the total number of hydrogen bonds is dominated by intersolvent hydrogen bond interactions. As different types of species (ethanol, sucrose, and NaHCO3) are added into the solution, the hydrogen bonding network of the fluid is noticeably disrupted. In particular, the total number of hydrogen bonds between water molecules is reduced by 4.3% in the champagnes and by 3.4% in the cola drinks (0.35 M) in comparison with the aqueous solution of CO2 (label 4 in Figure 2). Interestingly, in the NaHCO3 solution, a model for club sodas, we found that the number of hydrogen bonds, in fact, is increased by 0.5%. Conversely, the decrease in the number of hydrogen bonds of the solvent is compensated by the formation of a greater number of hydrogen bonds between water and the added components in the first two systems (label 5 in Figure 2). Overall, both the alcohol and sucrose additives gain more hydrogen bonds than those lost between solvent molecules, resulting in a net overall increase in the total number over that without the additives. The increase in the number of hydrogen bonds in the mixtures is consistent with the finding of higher viscosities.9 CO2 itself is not considered to form specific hydrogen bonds with water in the liquid, although recent studies showed dissolved carbon dioxide can produce observable changes in hydrogen bonding interactions of water.39 Nevertheless, we have counted the nearest neighbors in the proximity of carbon dioxide by the same criterion as the hydrogen-bond definition. Figure 2 (label 2) reveals that a greater number of water molecules are pushed closer to carbon dioxide in alcohol and sucrose solutions relative to that in

Figure 3. Distribution of the tetrahedral parameter q characterizing the hydrogen bonding environment about a central water molecule for an aqueous solution of CO2 (blue), 12% v/v ethanol−water mixture corresponding to typical champagnes (red), 0.36 M sucrose solution as in a typical can of cola drink (green), and 0.12 M sodium bicarbonate solution as in club sodas (purple).

of the high and low peaks are located at 0.8 and 0.5, respectively. For comparison, both ethanol and sucrose solutions show an enhanced distribution in the low-q value region accompanied by a decrease in the high-q region. There are also small shifts in the maximum positions at about 0.85 and 0.45. Of the two solutions, the change in the alcohol−water mixture is more pronounced than that in the sucrose solution, probably because of the steric effect. A similar steric effect was demonstrated by Matysiak et al.;38 they found that nonpolar groups of the polypeptide prevent the bulk water molecules D

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B from approaching the hydration water molecules, leading to an overcrowded environment and thereby causing distortions of hydrogen bonds among proximal water molecules. For the salt solution of NaHCO3, the q distribution exhibits an appreciable increase in the population of the high-q peak with a minimal shift in the peak position. The low-q peak is, similar to two other solutions, also shifted toward a smaller value. These features suggest that, for the ethanol and sucrose solutions, the ice-like, tetrahedral structural populations of the solvent are somewhat reduced, consistent with the finding that the number of hydrogen bonds is less than that without the additives. In contrast, the added salt causes a slight enhancement in the structured configurations, which may be called a structuremaking additive. Sciortino et al.45 have shown that there is some association between the relatively flexible dynamics of neat water and the presence of defects in the tetrahedral hydrogen bond network. In the current study, the structuremaking additives apparently cause stronger tetrahedrality than the water system, therefore giving reduced flexibility in the water dynamics. Water Reorientation Dynamics. Since the diffusion of the solute CO2 is related to water dynamics, we therefore computed the reorientation time correlation function of water molecules in these model systems using the following formula C2(t ) = ⟨P2[u(0) ·u(t )]⟩

diffusivity with that of water, we carried out additional MD simulations, in which 10 CO2 molecules were dissolved in a 0.1 M concentration of sodium perchlorate aqueous solution (18 NaClO4 in 10 000 water molecules). The computed diffusion coefficient of CO2 in NaClO4 aqueous solution is 2.37 × 10−9 m2/s, even greater than that in the reference system of aqueous solution (2.33 × 10−9 m2/s). We note that the associated countercation also plays a role in water dynamics but not in a direct way. For champagnes, the water reorientation dynamics of water also indicates an appreciable retardation (Figure 4). Rezus and

(4)

where P2 is the second-order Legendre polynomial and u(t) is the molecular orientation. For water molecules, the molecular orientation can be defined either along the OH bond vector or along the dipole moment vector, and in this study, we used the OH bond vector of water for a given instantaneous configuration at time t. C2(t) is accessible from time-resolved IR spectroscopy and from molecular dynamics simulations.46 Figure 3 shows that all solute additives in the carbon dioxide aqueous mixture reduce the reorientation dynamics of water, with the greatest effect found in club sodas. We attribute this finding to the ability of the HCO3− anion to form much stronger hydrogen bonds with water than the other neutral species (column 4 in Table 2). This is also aided by the formation of a greater number of hydrogen bonds and increased in tetrahedral ordering in this solution (Figure 3). At a qualitative level, the presence of HCO3− plays a role of strengthening the hydrogen bond network of water, making it less mobile and reducing CO2 diffusion. Laage and his coworkers investigated47 the effects of inorganic salts such as NaClO4, and Na2SO4 on the hydrogen bond reorientation, and found that these ions may accelerate or slow down water dynamics in dilute solutions depending on the strength of their interactions with water. Surprisingly, in diluted NaClO4 aqueous solution, unlike most other ions, water diffusion is accelerated because of its relatively weak hydrogen bonding interactions and a larger excluded volume, which is reflected by its smaller viscosity compared to pure water.48 The anomalous water diffusion in salt solution was also observed by Kim et al.49 Interestingly, Galamba50 found that the I− with a smaller charge density (large, univalent ion) exhibits an increased tetrahedral order compared with other halide anions. These results indicated that there is an association between tetrahedral order and reorientation dynamics. In contrast to the diluted NaClO4 aqueous solution, water diffusion is significantly slowed down in Na2SO4 aqueous solutions; an anion has a similar excluded volume but significantly stronger hydrogen bond strength with water. To establish a correlation of CO2

Figure 4. Water orientation time correlation function for an aqueous solution of CO2 (blue), 12% v/v ethanol−water mixture corresponding to typical champagnes (red), 0.36 M sucrose solution as in a typical can of cola drink (cyan), and 0.12 M sodium bicarbonate solution as in club sodas (purple).

Bakker51 have demonstrated that water molecules around the hydrophobic group of the solutes that retards the water. In the current study, ethanol molecules in champagnes are a hydrophobic part of the ethanol molecules as well. Such a rationalization is also supported by the research carried out by Laage and his co-workers,52 in which the origin of the water retardation is more attributed to the excluded volume of the solute that slows H-bond exchange. In addition, one may presume that these immobilized water molecules correspond to the presence of higher-order structures. However, according to our results (Figure 3), it may seem counterintuitive in that, even though ethanol turns out to be structure-destructing, the champagnes still possess slow water reorientation dynamics. Nevertheless, it was already found that, for hydrophobic groups, there is a not rigorous correspondence between tetrahedral order and water dynamics at the room temperature (but more prominent in low temperature).52 The cola drink exhibits the smallest water retardation among the three carbonated beverages, yet still higher than that in the benchmark system according to the water reorientation dynamics results. The mechanism is most likely also caused by the formation of the intramolecular hydrogen bond in the sucrose molecules, leading to hindered exchange of hydrogen bonds (Figure S3).



CONCLUSION In this study, we have carried out molecular dynamic simulations to determine the CO2 diffusion coefficients and E

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

by Classical Molecular Dynamics and 13C NMR Spectroscopy. J. Phys. Chem. Lett. 2014, 5 (24), 4232−4237. (10) Liger-Belair, G.; Parmentier, M.; Jeandet, P. Modeling the Kinetics of Bubble Nucleation in Champagne and Carbonated Beverages. J. Phys. Chem. B 2006, 110 (42), 21145−21151. (11) Liger-Belair, G.; Marchal, R.; Robillard, B.; Vignes-Adler, M.; Maujean, A.; Jeandet, P. Study of Effervescence in a Glass of Champagne: Frequencies of Bubble Formation, Growth Rates, and Velocities of Rising Bubbles. Am. J. Enol. Vitic. 1999, 50 (3), 317−323. (12) Nose, A.; Hojo, M.; Ueda, T. Effects of Salts, Acids, and Phenols on the Hydrogen-Bonding Structure of Water− Ethanol Mixtures. J. Phys. Chem. B 2004, 108 (2), 798−804. (13) Nose, A.; Hamasaki, T.; Hojo, M.; Kato, R.; Uehara, K.; Ueda, T. Hydrogen Bonding in Alcoholic Beverages (Distilled Spirits) and Water− Ethanol Mixtures. J. Agric. Food Chem. 2005, 53 (18), 7074− 7081. (14) Nose, A.; Hojo, M. Hydrogen Bonding of Water−Ethanol in Alcoholic Beverages. J. Biosci. Bioeneg. 2006, 102 (4), 269−280. (15) Lee, S. L.; Debenedetti, P. G.; Errington, J. R. A Computational Study of Hydration, Solution Structure, and Dynamics in Dilute Carbohydrate Solutions. J. Chem. Phys. 2005, 122 (20), 204511. (16) Laage, D.; Stirnemann, G.; Sterpone, F.; Rey, R.; Hynes, J. T. Reorientation and Allied Dynamics in Water and Aqueous Solutions. Annu. Rev. Phys. Chem. 2011, 62, 395−416. (17) Perret, A.; Bonhommeau, D. A.; Liger-Belair, G.; Cours, T.; Alijah, A. CO2 Diffusion in Champagne Wines: A Molecular Dynamics Study. J. Phys. Chem. B 2014, 118 (7), 1839−1847. (18) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30 (13), 2157−2164. (19) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, CA, 2002. (20) Kumar, R.; Schmidt, J.; Skinner, J. Hydrogen Bonding Definitions and Dynamics in Liquid Water. J. Chem. Phys. 2007, 126 (20), 204107. (21) Luzar, A.; Chandler, D. Hydrogen-Bond Kinetics in Liquid Water. Nature 1996, 379 (6560), 55. (22) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25. (23) Zeebe, R. E. On the Molecular Diffusion Coefficients of Dissolved CO2, HCO3− and CO32‑ and Their Dependence on Isotopic Mass. Geochim. Cosmochim. Acta 2011, 75 (9), 2483−2498. (24) Mahoney, M. W.; Jorgensen, W. L. Diffusion Constant of the TIP5P Model of Liquid Water. J. Chem. Phys. 2001, 114 (1), 363−366. (25) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N× Log (N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089−10092. (26) Evans, D. J.; Holian, B. L. The Nose−Hoover Thermostat. J. Chem. Phys. 1985, 83 (8), 4069−4074. (27) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182−7190. (28) Wang, B.; Kuo, J.; Bae, S. C.; Granick, S. When Brownian Diffusion is not Gaussian. Nat. Mater. 2012, 11 (6), 481. (29) Kwon, G.; Sung, B. J.; Yethiraj, A. Dynamics in Crowded Environments: Is Non-Gaussian Brownian Diffusion Normal? J. Phys. Chem. B 2014, 118 (28), 8128−8134. (30) Wang, B.; Anthony, S. M.; Bae, S. C.; Granick, S. Anomalous yet Brownian. Proc. Natl. Acad. Sci. U. S. A. 2009, 106 (36), 15160−15164. (31) Guan, J.; Wang, B.; Granick, S. Even Hard-Sphere Colloidal Suspensions Display Fickian yet Non-Gaussian Diffusion. ACS Nano 2014, 8 (4), 3331−3336. (32) Van Hove, L. Correlations in Space and Time and Born Approximation Scattering in Systems of Interacting Particles. Phys. Rev. 1954, 95 (1), 249−262.

the origin of its variation in three model systems with characteristic features of champagne, sugar-based soft drinks, and club sodas. The calculated CO2 diffusion coefficients are in reasonable agreement with the experimental values. The diffusion coefficient of CO2 is reduced when any of ethanol, sucrose, or NaHCO3 additives is introduced in the solution, though to different extents. Analyses of hydrogen bonding and the solvent’s structural and dynamic properties reveal that the change of the CO2 diffusion coefficient is closely associated with the diffusional behavior of the solvent water itself. This is also supported by analysis of a different CO2 accelerating system that includes sodium perchlorate as the additive.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b10469. The ingredients of a typical cola and a typical club soda; diffusion coefficient of CO2 in different system sizes; Van Hove correlation functions for these carbonated beverages (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Ji Lv: 0000-0002-3693-8979 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This research has received funding support from the National Natural Science Fundation of China (21673092 and 21533003). The authors also particularly appreciate Jiali Gao (Jilin University & University of Minnesota) for his instructive discussion, technical support, and careful proofreading.



REFERENCES

(1) Barker, G. S.; Jefferson, B.; Judd, S. J. The Control of Bubble Size in Carbonated Beverages. Chem. Eng. Sci. 2002, 57 (4), 565−573. (2) Liger-Belair, G.; Jeandet, P. Effervescence in a Glass of Champagne: A Bubble Story. Europhys. News 2002, 33 (1), 10−14. (3) Liger-Belair, G.; Polidori, G.; Jeandet, P. Recent Advances in the Science of Champagne Bubbles. Chem. Soc. Rev. 2008, 37 (11), 2490− 2511. (4) Liger-Belair, G. Modeling the Losses of Dissolved CO2 from Laser-Etched Champagne Glasses. J. Phys. Chem. B 2016, 120 (15), 3724−3734. (5) Zhang, Y.; Xu, Z. Fizzics” of Bubble Growth in Beer and Champagne. Elements 2008, 4 (1), 47−49. (6) Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.; Michaelides, A. Crystal Nucleation in Liquids: Open Questions and Future Challenges in Molecular Dynamics Simulations. Chem. Rev. 2016, 116 (12), 7078−7116. (7) Sharma, K.; Adhikari, N. P. Temperature Dependence of Diffusion Coefficient of Nitrogen Gas in Water: A Molecular Dynamics Study. Int. J. Mod. Phys. B 2014, 28 (14), 1450084. (8) Liger-Belair, G.; Prost, E.; Parmentier, M.; Jeandet, P.; Nuzillard, J.-M. Diffusion Coefficient of CO2Molecules as Determined by 13C NMR in Various Carbonated Beverages. J. Agric. Food Chem. 2003, 51, 7560−7563. (9) Bonhommeau, D. A.; Perret, A.; Nuzillard, J.-M.; Cilindre, C.; Cours, T.; Alijah, A.; Liger-Belair, G. r. Unveiling the Interplay Between Diffusing CO2 and Ethanol Molecules in Champagne Wines F

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (33) Himmelblau, D. Diffusion of Dissolved Gases in Liquids. Chem. Rev. 1964, 64 (5), 527−550. (34) Vega, C.; Abascal, J. L. F.; Conde, M. M.; Aragones, J. L. What Ice Can Teach Us About Water Interactions- A Critical Comparison of the Performance of Different Water Models. Faraday Discuss. 2009, 141, 251−276. (35) Belgodere, C.; Dubessy, J.; Vautrin, D.; Caumon, M. C.; Sterpenich, J.; Pironon, J.; Robert, P.; Randi, A.; Birat, J. P. Experimental Determination of CO2 Diffusion Coefficient in Aqueous Solutions under Pressure at Room Temperature via Raman Spectroscopy: Impact of Salinity (NaCl). J. Raman Spectrosc. 2015, 46 (10), 1025−1032. (36) Cadogan, S. P.; Maitland, G. C.; Trusler, J. M. Diffusion Coefficients of CO2 and N2 in Water at Temperatures between 298.15 and 423.15 K at Pressures up to 45 MPa. J. Chem. Eng. Data 2014, 59 (2), 519−525. (37) Molinero, V.; Ç aǧın, T.; Goddard, W. A., III Sugar, Water and Free Volume Networks in Concentrated Sucrose Solutions. Chem. Phys. Lett. 2003, 377 (3−4), 469−474. (38) Matysiak, S.; Debenedetti, P. G.; Rossky, P. J. Dissecting the Energetics of Hydrophobic Hydration of Polypeptides. J. Phys. Chem. B 2011, 115 (49), 14859−14865. (39) Wolke, C. T.; Fournier, J. A.; Dzugan, L. C.; Fagiani, M. R.; Odbadrakh, T. T.; Knorke, H.; Jordan, K. D.; McCoy, A. B.; Asmis, K. R.; Johnson, M. A. Spectroscopic Snapshots of the Proton-Transfer Mechanism in Water. Science 2016, 354 (6316), 1131−1135. (40) Duboue-Dijon, E.; Laage, D. Characterization of the Local Structure in Liquid Water by Various Order Parameters. J. Phys. Chem. B 2015, 119 (26), 8406−18. (41) Hellströ m, M.; Behler, J. Structure of Aqueous NaOH Solutions: Insights from Neural-Network-Based Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2017, 19 (1), 82−96. (42) Galamba, N. Water’s Structure Around Hydrophobic Solutes and the Iceberg Model. J. Phys. Chem. B 2013, 117 (7), 2153−2159. (43) Kumar, P.; Buldyrev, S. V.; Stanley, H. E. A tetrahedral Entropy for Water. Proc. Natl. Acad. Sci. U. S. A. 2009, 106 (52), 22130−22134. (44) Errington, J. R.; Debenedetti, P. G. Relationship between Structural Order and the Anomalies of Liquid Water. Nature 2001, 409 (6818), 318. (45) Sciortino, F.; Geiger, A.; Stanley, H. E. Effects of Defects on Molecular Mobility in Liquid Water. Nature 1991, 354 (6350), 218− 221. (46) Zhang, L.; Greenfield, M. L. Relaxation Time, Diffusion, and Viscosity Analysis of Model Asphalt Systems Using Molecular Simulation. J. Chem. Phys. 2007, 127 (19), 194502. (47) Stirnemann, G.; Wernersson, E.; Jungwirth, P.; Laage, D. Mechanisms of Acceleration and Retardation of Water Dynamics by Ions. J. Am. Chem. Soc. 2013, 135 (32), 11824−11831. (48) Janssen, L. Physical Properties of Aqueous Solutions Containing NaClO3 and/or NaClO4. J. Appl. Electrochem. 1995, 25 (3), 291−291. (49) Kim, J. S.; Wu, Z.; Morrow, A. R.; Yethiraj, A.; Yethiraj, A. SelfDiffusion and Viscosity in Eelectrolyte Solutions. J. Phys. Chem. B 2012, 116 (39), 12007−12013. (50) Galamba, N. Mapping Structural Perturbations of Water in Ionic Solutions. J. Phys. Chem. B 2012, 116 (17), 5242−5250. (51) Rezus, Y. L. A.; Bakker, H. J. Observation of Immobilized Water Molecules around Hydrophobic Groups. Phys. Rev. Lett. 2007, 99 (14), 148301. (52) Duboué-Dijon, E.; Fogarty, A. C.; Laage, D. Temperature Dependence of Hydrophobic Hydration Dynamics: From Retardation to Acceleration. J. Phys. Chem. B 2014, 118 (6), 1574−1583. (53) In Het Panhuis, M.; Patterson, C. H.; Lynden-Bell, R. M. A Molecular Dynamics Study of Carbon Dioxide in Water: Diffusion, Structure and Thermodynamics. Mol. Phys. 1998, 94 (6), 963−972. (54) Jähne, B.; Heinz, G.; Dietrich, W. Measurement of the Diffusion Coefficients of Sparingly Soluble Gases in Water. J. Geophys. Res. 1987, 92 (C10), 10767−10776.

(55) Zhang, W.; Wu, S.; Ren, S.; Zhang, L.; Li, J. The Modeling and Experimental Studies on the Diffusion Coefficient of CO2 in Saline Water. J. CO2 Util. 2015, 11, 49−53. (56) Palmer, B. J. Transverse-Current Autocorrelation-Function Calculations of the Shear Viscosity for Molecular Liquids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1994, 49 (1), 359− 366.

G

DOI: 10.1021/acs.jpcb.7b10469 J. Phys. Chem. B XXXX, XXX, XXX−XXX