Toward In Silico Prediction of CO2 Diffusion in Champagne Wines.

Carbon dioxide diffusion is the main physical process behind the formation and growth of bubbles in sparkling wines, especially champagne wines. By approximating brut-labeled champagnes as carbonated hydroalcoholic solutions, molecular dynamics (MD) simulations are carried out with six rigid water models and three CO2 models to evaluate CO2 diffusion coefficients. MD simulations are little sensitive to the CO2 model but proper water modeling is essential to reproduce experimental measurements. A satisfactory agreement with nuclear magnetic resonance (NMR) data is only reached at all temperatures for simulations based on the OPC and TIP4P/2005 water models; the similar efficiency of these two models is attributed to their common properties such as low mixture enthalpy, same number of hydrogen bonds, alike water tetrahedrality, and multipole values. Correcting CO2 diffusion coefficients to take into account their system-size dependence does not significantly alter the quality of the results. Estimates of viscosities deduced from the Stokes-Einstein formula are found in excellent agreement with viscometry on brut-labeled champagnes, while theoretical densities tend to underestimate experimental values. OPC and TIP4P/2005 water models appear to be choice water models to investigate CO2 solvation and transport properties in carbonated hydroalcoholic mixtures and should be the best candidates for any MD simulations concerning wines, spirits, or multicomponent mixtures with alike chemical composition.

Introduction

Champagne wines are multicomponent aqueous solutions composed of ethanol (12.5% v/v), dissolved carbon dioxide (10–12 g L–1), sugars (≲50 g L–1), a broad variety of ions (e.g., K+, Ca2+, and Cl–), and a multitude of complex organic compounds.[1] Under standard tasting conditions, after uncorking a bottle of champagne, the supersaturation of the liquid phase with diffusing CO2 molecules results in the formation of bubbles by heterogeneous nucleation at the vicinity of cavities that may be salt crystals, glass scratches, or tiny vegetal pieces like cellulose fibers.[2] More precisely, the surface of a glass is always scattered with cellulose fibers coming from the environment. When a glass is poured with champagne, a gas pocket is trapped within the hollow (called lumen) of the hydrated cellulose fiber. This gas pocket grows due to the diffusion of CO2 molecules through the wall of the cellulose fiber until its size is large enough for enabling the release of a CO2 bubble at the fiber edge, as depicted in Figure a–d. The subsequent bubble dynamics is governed by the ability of bulk CO2 molecules to penetrate into the newly born bubble. This additional amount of CO2 makes the bubble grow, accelerate through buoyancy, and rise up to the free surface of the liquid, as illustrated in Figure e.[3,4]

Figure 1

Successive steps of the life cycle of a CO2 bubble in a glass of champagne, from its formation by heterogeneous nucleation in a cellulose fiber (a–d) to its release, growth, and rise in the champagne bulk (e).

CO2 diffusion is therefore the main physical process behind the formation and growth of bubbles in champagne wines, and more generally sparkling beverages. From the experimental point of view, the diffusion coefficients of CO2 in champagnes can be derived from nuclear magnetic resonance (NMR) spectroscopy measurements,[5−7] by means of the Stejskal–Tanner equation[8] that relates the intensity of NMR spectra to the diffusion coefficient of species in the sample or by applying the Stokes–Einstein formula provided that a value of the dynamic viscosity of the liquid and an estimate of the CO2 hydrodynamic radius are available.[9] In particular, 13C NMR measurements performed on a brut-labeled champagne (concentration of sugars below 12 g L–1) at 293 K yielded DCO(293 K) = 1.41 × 10–9 m2 s–1, a value similar to that obtained for a beer and another sparkling wine, but higher than for sodas and lower than for fizzy water.[5] The same order of magnitude was obtained by Autret et al. by nondestructive NMR measurements carried out on sealed bottles containing two different brut-labeled champagnes at 295 K, DCO(295 K) = (1.3 ± 0.1) × 10–9 m2 s–1, despite a higher noise on NMR spectra in these experiments.[6] More recently, accurate series of 13C NMR measurements on brut-labeled champagnes at temperatures ranging from 4 °C (fridge temperature) to 20 °C (ambient temperature) lead to values of CO2 diffusion coefficients extending from DCO(277 K) = 0.795 ± 0.014 × 10–9 m2 s–1 to DCO(293 K) = 1.287 ± 0.004 × 10–9 m2 s–1,[7] in agreement with experimental data by Autret et al.

From the theoretical point of view, the molecular modeling of multicomponent systems like champagnes may seem involved but these mixtures can be regarded as carbonated hydroalcoholic solutions in first approximation. In this context, Perret et al. carried out molecular dynamics simulations of a carbonated water/ethanol mixture that respected brut-labeled champagne proportions in the isothermal–isobaric ensemble.[10] Water molecules were modeled either by the 3-site SPC/E model[11] or the 5-site TIP5P model,[12] CO2 molecules by the popular EPM2 model,[13] and ethanol (EtOH) molecules by default parameters of the CHARMM27 force field.[14] By assuming champagnes as homogeneous and isotropic liquids on average, they were able to get CO2 diffusion coefficients over the whole experimental temperature range, namely, from 277 to 293 K. Values obtained at temperatures above 285 K with the SPC/E water model were in close agreement with NMR measurements performed on carbonated hydroalcoholic solutions and brut-labeled champagnes.[7] However, the agreement was much more questionable at low temperatures (T < 285 K) where theoretical CO2 diffusion coefficients strongly underestimated the experimental value at T = 277 K and overestimated it at T = 281 K. Moreover, CO2 diffusion coefficients obtained with the TIP5P water model overestimated experimental data over the whole temperature range, a result later confirmed by alike studies conducted by Lv et al.,[15] although the temperature dependence seemed qualitatively correct. Finally, convergence issues due to the relatively short duration of the production runs (i.e., 1 ns) motivated the authors to employ replica exchange dynamics, a parallel approach that might not be needed to get accurate values of CO2 diffusion coefficients in a hydroalcoholic solution.

In the present work, a comprehensive study of CO2 diffusion coefficients in carbonated hydroalcoholic solutions is undertaken as a function of temperature in an attempt to identify the most suitable molecular models to describe the diffusion of CO2 in brut-labeled champagnes. Six water models and three CO2 models are compared with each other before discussing deviations from experimental expectations in terms of enthalpy, number of hydrogen (H) bonds, tetrahedral arrangement of water molecules, and water multipole moments. Recommendations to build a reliable model for carbonated hydroalcoholic solutions representative of champagnes, and more generally, sparkling wines, are eventually supplied as concluding remarks.

Results and Discussion

Influence of Water Models

Molecular dynamics (MD) simulations based on classical force fields are practical tools to evaluate CO2 diffusion coefficients in water,[16,17] brines for applications in CO2 capture and sequestration,[18] and sparkling wines.[7,10,15] In the latter field of research, a particular emphasis is brought to the interactions between CO2 molecules and the other species of the mixture since CO2 is responsible for the production of bubbles in sparkling beverages. Such a carbonated mixture being mainly composed of water (∼95% of the quantity of matter), the water model in use should have a significant influence on the motion of molecules within the liquid. Figure depicts CO2 diffusion coefficients obtained for six rigid water models coupled with CO2 molecules described by the popular EPM2 model and EtOH molecules modeled by the OPLSaa force field.[19] The temperature dependence of all simulations is qualitatively similar, CO2 diffusion coefficients steadily increasing with temperature, and the results are not qualitatively influenced by the approach followed to evaluate diffusion coefficients. However, TIP4P/2005[20] and OPC[21] are the only water models in almost quantitative agreement with experimental data extracted from 13C NMR measurements on carbonated hydroalcoholic samples at all temperatures. Other water models yield diffusion coefficients that evolve in the vicinity of former TIP5P calculations that were known to overestimate the experimental values.[7] Although OPC provides the best agreement with experimental results to date, SPC/E replica exchange MD simulations carried out on shorter time scales (1 ns instead of 10 ns here) with more CO2 molecules (500 instead of 200), and different force field parameters for EtOH molecules (CHARMM27 force field instead of OPLSaa force field), were noticeably accurate at temperatures T ≥ 285 K, a behavior that is not reproduced in the present study with the SPC/E model. The convergence of our results being checked with respect to the number of CO2 molecules and simulation time, we attribute this deviation to variations in ethanol parameters.

Figure 2

CO2 diffusion coefficients in a carbonated hydroalcoholic solution as a function of temperature for six water models (circles) coupled with the EPM2 CO2 model. Diffusion coefficients are derived from the linearization of CO2 MSDs (red circles) and the integration of CO2 VACFs (green circles). These results are compared with experimental data (black upward triangles) and previous MD calculations (black squares) based on the SPC/E (black solid line) and TIP5P (black dashed line) water models.[7] The legend indicated in the top-left figure applies to the six figures. Reprinted (Adapted or Reprinted in part) with permission from ref (7). Copyright 2014 American Chemical Society.

The overestimated CO2 diffusion coefficients in carbonated hydroalcoholic mixtures containing TIP5P or SPC/E water molecules should be evidenced by energetic or structural properties of the liquid such as its enthalpy, its number of H bonds, or the average water tetrahedral order parameters (TOPs). According to Figure a, the lowest enthalpy is obtained for the mixture containing OPC water molecules, whereas the highest one is obtained for that containing TIP5P water molecules, which announces more cohesion of the former mixture and a subsequent slower molecular diffusion, in agreement with CO2 diffusion coefficients plotted in Figure . Other models roughly abide by the previous observations, although we would have intuitively expected TIP4P/2005 and TIP4P-Ew to yield lower enthalpies than TIP5P/2018 and SPC/E, respectively, if the previous rule were rigorously respected. The temperature dependence of enthalpy can be supplemented by the count of H bonds. Figure b shows that OPC and TIP4P/2005 have the maximum average number of H bonds at all temperatures, followed by TIP4P-Ew, SPC/E, TIP5P/2018, and TIP5P, which gets ∼3–6% fewer H bonds than other models. These observations are compatible with enthalpies, although we cannot exactly correlate the two quantities since energies depend not only on the number of H bonds but also on the well depths and equilibrium distances of all interactions involved in the mixture. As interactions between water molecules contribute to more than 90% of H bonds,[10] the microscopic arrangement of water molecules can provide additional information on the molecular networks in place in our simulations. The average radial part of water TOPs, ⟨Sk⟩, depicted in Figure c is a measure of the variance of O–O bonds between a central water molecule and its four nearest neighbors, a value of zero corresponding to the perfect tetrahedron.[22] Temperature little influences this parameter that remains roughly constant at values in between 1.73 × 10–3 for OPC and 2.31 × 10–3 for TIP5P. OPC and TIP4P/2005 are the models with little variance, making their molecular networks less flexible, and the ranking of water models is exactly that we would have deduced from CO2 diffusion coefficients if a correlation had been established between the flexibility of the water network and CO2 diffusion. The qualitative trends of the average angular part of water TOPs, ⟨Sg⟩, plotted in Figure d are similar for all of the models except TIP5P, which is subjected to a sharper increase. Tetrahedrality is slowly lost as temperature increases, and the water network built from OPC and TIP4P/2005 has the highest tetrahedral characters at all temperatures. The agreement with 13C NMR data of CO2 diffusion coefficients derived from MD simulations using the OPC and TIP4P/2005 water models demonstrate that a low enthalpy, high number of H bonds, and low value of TOPs are the required conditions to hope for any proper modeling of CO2 diffusion in carbonated hydroalcoholic solutions. Moreover, a number of past studies, including some of ours, considered the TIP5P water model as a candidate to investigate molecular diffusion in carbonated beverages,[7,10,15] a habit that should be discarded according to these results.

Figure 3

Energy and structural properties of carbonated hydroalcoholic mixtures for six water models coupled with the EPM2 CO2 model as a function of temperature. (a) Enthalpy, (b) number of H bonds, (c) average radial part ⟨Sk⟩ of water TOPs, and (d) average angular part ⟨Sg⟩ of water TOPs.

To pinpoint more finely possible reasons for the efficiency of OPC and TIP4P/2005, we compare the strategies followed to build the six water models. These strategies may differ from each other by the computational details (e.g., box size, cutoff distances of Coulomb and Lennard–Jones (LJ) interactions, and long-range corrections) of benchmark Monte Carlo or MD simulations, by the target properties the model should accurately reproduce, and by the fitting process used to optimize the parameters of the model. Most computational details cannot unambiguously explain differences observed in Figure since cutoff radii are similar for the six models (∼0.8–0.9 nm) and the model converged from simulations performed in the biggest simulation box (TIP5P/2018 with N = 2069) does not yield the best results. Moreover, the addition of long-range corrections to nonbonding interactions is known to produce discrepancies, such as lower densities, in TIP5P simulations[23] but this drawback should not apply to TIP4P-Ew, TIP4P/2005, TIP5P/2018, and OPC that include long-range corrections during the convergence process of their parameters. All of the models, except OPC, impose the geometry of the water molecule (OH bond length and H–O–H angle) and use target water properties like energies of vaporization, liquid densities, or the temperature of maximum density to optimize the water partial charges and the O–O LJ parameters. The originality of OPC lies in the special care born by Izadi et al.[21] to optimize the three lower-order water multipole moments from target bulk properties to get a water geometry able to improve predictions of solvation free energies of small molecules. The average root-mean-squared (rms) deviation between the multipole moments obtained for OPC and for the five other water models is illustrated in Figure . SPC/E, TIP5P, and TIP5P/2018 values deviate from OPC ones by at least 0.13 and they are significantly modified by the account of octopoles. On the contrary, TIP4P-Ew and TIP4P/2005 exhibit much smaller rms deviations to OPC (ΔξOPC ≲ 0.09) and are not sensitive to the addition of octopoles. A proper electrostatic description through accurate water multipole moments is probably required to improve the agreement with experiments but the fact that TIP4P-Ew has an rms deviation to OPC smaller than TIP4P/2005 while yielding less accurate CO2 diffusion coefficients emphasizes that a proper description of multipole moments is only part of the solution. In carbonated hydroalcoholic mixtures, the more satisfactory behavior of TIP4P/2005, based on the same water geometry as TIP4P-Ew, may come from parameter adjustments on both liquid water and ice properties, which makes this model applicable to a broader class of systems. In particular, TIP4P/2005 is known to reproduce phase diagrams much more accurately than TIP4P-Ew.[20]

Figure 4

Polar diagram for the rms error ΔξOPC when dipoles and quadrupoles (d+q) or dipoles, quadrupoles, and octopoles (d+q+o) are included in the definition of ΔξOPC. The value of ΔξOPC matches the distance to the diagram center.

Influence of Carbon Dioxide Models

Although OPC and TIP4P/2005 are the water models leading to the best agreement with NMR experiments when combined with EPM2 CO2 molecules, it might be relevant to evaluate the validity of this conclusion as CO2 is described by other molecular models. As an example, Figure gathers theoretical CO2 diffusion coefficients deduced from MD simulations carried out with OPC and three CO2 models, namely, EPM2, TraPPE,[24] and MSM-ZD,[25] by constraining bonds or leaving them free to vibrate. A close agreement with the experimental curve is reached for the EPM2 and TraPPE models whatever the constraints on bonds. Interestingly, the quality of the results is slightly lowered for MSM-ZD while this model was originally devised to better reproduce H2O–CO2 interactions on a wide range of temperatures and pressures. One possible explanation is that the authors were particularly concerned by geochemical applications commonly modeled by the SPC/E water model. Combining MSM-ZD CO2 and SPC/E water might improve the description of the mixture but we do not expect the occurrence of accurate CO2 diffusion coefficients because water self-diffusion is already overestimated by 8% in pure SPC/E water[11,26] and the addition of EtOH in the mixture might yield an additional overestimation of CO2 diffusion coefficients, as illustrated in Figure when the EPM2 model is used and supported by Figure where the properties of mixtures containing SPC/E water molecules differ from their OPC or TIP4P/2005 counterparts. It is worth noting that the small influence of the CO2 model on our results may come from the low concentration of carbon dioxide in the carbonated hydroalcoholic solutions considered here that respect the typical proportions of champagnes. The same conclusion should therefore hold true for all of the sparkling wines.

Figure 5

CO2 diffusion coefficients in a carbonated hydroalcoholic solution as a function of temperature for three CO2 models coupled with the OPC water model, by constraining bonds (red open circles) or leaving them free to vibrate (blue filled circles). Experimental data are represented by black upward triangles.[7] Reprinted (Adapted or Reprinted in part) with permission from ref (7). Copyright 2014 American Chemical Society.

CO2 Diffusion in Brut-Labeled Champagnes

Due to their low concentration in sugars (<12 g/L), brut-labeled champagne wines can be considered as carbonated hydroalcoholic solutions in first approximation as confirmed by recent 13C NMR measurements.[7] CO2 diffusion coefficients, viscosities, and densities extracted from TIP4P/2005 and OPC simulations are compared with experimental data obtained for carbonated hydroalcoholic and brut-labeled champagne samples in Figure . The dynamic viscosity of a mixture is experimentally determined by multiplying its kinematic viscosity by its liquid density. From the theoretical point of view, viscosity can be obtained by computing the transverse-current autocorrelation functions (TCAF) but a simpler approach consists in deducing its value from the Stokes–Einstein formula provided that the diffusion coefficient and hydrodynamic radius of one species of the mixture are known. The agreement with experiments is excellent at all temperatures when the Stokes–Einstein relation is used, TIP4P/2005 calculations only underestimating the experimental viscosity at 277 K by 0.36 × 10–3 Pa·s. This confirms the practical interest of the Stokes–Einstein relation to evaluate viscosities provided that accurate diffusion coefficients are available. On the contrary, viscosities evaluated from TCAF underestimate systematically experimental values obtained for carbonated hydroalcoholic solutions by 10–20% when water molecules are described by the TIP4P/2005 and OPC models. However, these viscosities have the proper order of magnitude and would yield very similar corrections to CO2 diffusion coefficients as those predicted by viscosities estimated from the Stokes–Einstein relation (see the Supporting Information). Incidentally, theoretical CO2 diffusion coefficients are in good agreement with experimental diffusion coefficients obtained in brut-labeled champagnes, as expected due to the closeness between the experimental curves specific to carbonated hydroalcoholic solutions and brut-labeled champagnes, and correcting these diffusion coefficients to take into account their dependence on the system size slightly degrades the agreement with experiments but the results remain within the uncertainties of the original MD calculations. Finally, we must recognize that OPC and TIP4P/2005 water models do not nicely describe all of the properties of carbonated hydroalcoholic mixtures representative of brut-labeled champagnes. For instance, the density is underestimated by the two models, although the experimental trend is maintained as the temperature increases; the best agreement with experimental densities is obtained for TIP4P/2005.

Figure 6

Experimental diffusion coefficients, viscosities, and densities of carbonated hydroalcoholic solutions (CHS) and brut-labeled champagnes (BC) compared with recommended theoretical data as a function of temperature. Theoretical diffusion coefficients corrected to take into account system-size dependence are depicted with dashed lines. Theoretical viscosities are derived from the calculation of transverse-current autocorrelation functions (TCAF) and from the Stokes–Einstein relation where the CO2 hydrodynamic radius is identified with its rms radius (SE-RMS). Reprinted (Adapted or Reprinted in part) with permission from ref (7). Copyright 2014 American Chemical Society.

Conclusions

The ability of a water model to properly reproduce water properties on a wide range of thermodynamic conditions is not necessarily representative of its ability to describe the properties of a multicomponent mixture, even if this mixture is mainly composed of water. As an example, the ancient TIP3P model[27] is known to predict more accurate hydration free energies of small neutral organic molecules than TIP5P or TIP4P-Ew.[28,29] In the present study, the efficiency of six water models and three carbon dioxide models liable to be good candidates for describing CO2 diffusion in carbonated hydroalcoholic solutions is investigated. Although CO2 models were found to have little effect on results, the choice of the water model is of paramount importance to reproduce experimental data extracted from 13C NMR measurements. The OPC and TIP4P/2005 models were found to yield the more accurate CO2 diffusion coefficients at all temperatures relevant for applications on champagne wines, while the TIP5P results were less reliable. This conclusion is not altered by correcting diffusion coefficients for system-size effects even if the agreement with the experiments is slightly degraded. Mixtures built from OPC and TIP4P/2005 water molecules were found to exhibit the same number of H bonds, low enthalpy, and high water tetrahedrality, which should account for the small diffusivity of CO2 compared with mixtures built from other water models. Moreover, OPC and TIP4P/2005 have very similar values for their lower-order multipole moments, although TIP4P/2005 parameters were not specifically optimized to improve the accuracy of these quantities. For a deeper understanding of the reasons why OPC or TIP4P/2005 reproduce so nicely the experimental results compared with other models, like TIP5P or SPC/E, might require the calculation of CO2 solvation energies or the critical analysis of other force field parameters for ethanol, a work that is beyond the scope of this study.

Liquid viscosities estimated from CO2 diffusion coefficients and hydrodynamic radii by applying the Stokes–Einstein relation are in excellent agreement with experimental data obtained for carbonated hydroalcoholic mixtures. Since brut-labeled champagnes can be approximated as carbonated hydroalcoholic mixtures due to their low concentration in sugar, it is also reasonable to get a good agreement between our theoretical data on carbonated hydroalcoholic mixtures and the corresponding experimental data on brut-labeled champagnes, provided that water molecules are properly modeled. These results open new avenues in the sparkling alcoholic beverage industry where OPC and TIP4P/2005 could be used as reference water models to tackle problems such as the evaporation of aerosols on top of glasses in analogy with marine aerosols,[30,31] the CO2 diffusion through the wall of cellulose fibers,[32] or the influence of ethanol and sugar concentrations on solvation and transport properties in the mixture. Besides these enological questions, hydroalcoholic mixtures are also common solvents in chemistry. Water/methanol mixtures are for instance typical solvents in electrospray ionization (ESI) experiments for applications in proteomics or genomics,[33,34] and MD simulations are sometimes employed to investigate the fragmentation process of multicharged hydroalcoholic droplets.[35] However, the water model most accurate to predict transport properties such as diffusion coefficients or viscosities is not necessarily the most suitable for modeling desolvation processes, although TIP4P/2005 was already proved useful to describe the ejection of proteins from charged droplets.[36] More work in this direction is needed, and wines or spirits are ideal prototype systems to perform both simulations and experiments on various thermodynamical conditions and ethanol concentrations with possible industrial collaborations.

Methods

Force Fields

Six non-polarizable rigid water models are considered in our calculations: the 3-site SPC/E model,[11] the 4-site TIP4P-Ew,[37] TIP4P/2005,[20] and OPC[21] models, and the 5-site TIP5P[12] and TIP5P/2018[38] models. The TIP4P-Ew, TIP4P/2005, OPC, and TIP5P/2018 models were primarily selected because of their satisfactory description of water self-diffusion at 298 K. These four models yield water self-diffusion coefficients that differ from the benchmark experimental values (DHexp1 = 2.23 × 10–9 m2 s–1 and DHexp2 = 2.299 × 10–9 m2 s–1)[39,40] by no more than 9%; the best agreements are obtained for OPC (DH = 2.30 ± 0.02 × 10–9 m2 s–1)[21] and TIP5P/2018 (DH = 2.34 ± 0.02 × 10–9 m2 s–1), two models whose parameters are converged using water diffusivity as a target property. Simulations with the SPC/E and TIP5P water models are only carried out for the sake of comparison with former theoretical studies.[7,10,15]

Three molecular models are considered to describe CO2: the EPM2 model adjusted to reproduce the liquid–vapor coexistence curve and critical properties of pure CO2,[13] the TraPPE force field[24] devised to describe binary and ternary mixtures involving CO2, N2, and alkanes, and the more recent model introduced by Zhang and Duan[25] devised for possible use in conjunction with SPC/E water molecules to tackle industrial and geochemical problems where a proper description of CO2–H2O interactions is required. Parameters of the latter model hardly depart from those of the MSM model proposed by Murthy et al.,[41,42] and we have therefore called this model MSM-ZD in this paper.

EtOH is parameterized on the basis of the OPLSaa force field,[19] unlike former theoretical studies that employed the CHARMM27 force field.[14] This choice was motivated by the native parameterization of EtOH in OPLSaa.

Molecular Dynamics Simulations

MD simulations are performed with GROMACS open-source software (2018 versions)[43−47] in the NVT and NPT ensembles at five temperatures representative of champagne storage and tasting conditions, namely, temperatures between 277 and 293 K by steps of 4 K. Although brut-labeled champagnes are multicomponent mixtures, they can be modeled in first approximation as carbonated hydroalcoholic solutions composed of 4 × 104 water molecules, 1760 EtOH molecules, and 200 CO2 molecules. Glycerol, lactic acid, tartaric acid, and sugars are the next most abundant molecules as evidenced by their typical concentrations reported in Table . Concentrations of small ions or macromolecules (e.g., proteins and amino acids) are too low to be included in a simulation box composed of 4 × 104 water molecules or to have any influence on transport properties like CO2 diffusion. Therefore, focusing on the three most abundant molecules in brut champagnes (i.e., water, EtOH, and CO2), a cubic box of a side length of ∼11.1 nm subject to periodic boundary conditions can be built. LJ and electrostatic interactions are truncated at 1.5 nm, and smooth particle-mesh Ewald (SPME) summation techniques are applied for long-range electrostatic interactions. LJ pair well depth and diameters between unlike atoms are inferred from geometric means, as advocated by the OPLSaa force field, and bonds are generally constrained.

Table 1

Typical Composition of Champagne Wines: Concentrations and Corresponding Numbers of Molecules in a Simulation Box Composed of 4 × 104 Water Molecules

species	concentration	number
water	main species	40 000
ethanol	∼12.5% vol/vol	1760
carbon dioxide	10–12 g L – 1	188–224
glycerol	∼5 g L–1	45
tartaric acid	2.5–4 g L – 1	14–22
lactic acid	∼4 g L–1	37
sugars	0–50 g L–1	0–229
VOCsa	∼0.7 g L–1	1–4
proteins	5–10 mg L – 1	0
polysaccharides	∼0.2 g L–1	0
polyphenols	∼0.1 g L–1	0
amino acids	0.8–2 mg L – 1	0
lipids	∼10 mg L–1	0
K+	0.2–0.45 g L – 1	4–9
Ca2+	60–120 mg L – 1	1–2
Mg2+	50–90 mg L – 1	1–3
SO42–	∼0.2 g L–1	2
Cl–	∼10 mg L–1	0

Volatile organic compounds. Reproduced from ref (1) with permission from the Royal Society of Chemistry.

A typical MD calculation is a three-step simulation composed of a 1 ns NVT equilibration run followed by a 19 ns NPT equilibration run at a pressure of 1 bar, and an additional 10 ns NPT production run at the same pressure. The temperature and pressure are maintained with a Nosé–Hoover thermostat (with a time constant of 0.5 ps) and a Parinello–Rahman barostat (with a time constant of 0.2 ps), respectively. Simulations in the isothermal–isobaric ensemble need the knowledge of the isothermal compressibility of the mixture as a function of temperature. Such data being unavailable in the literature, water isothermal compressibilities derived from density analysis achieved by Vedamuthu et al.[48,49] were used instead. Atomic positions and velocities are stored every picosecond during the production run, trajectories are visualized with VMD software (version 1.9.3),[50] and high-quality pictures of molecular systems are generated with Protein Imager.[51] Temperatures, pressures, densities, and diffusion coefficients are averaged over the whole production run. Restricting to the first 100 ps of the production run the averages on the number of hydrogen bonds (≳7 × 104 by snapshot) and water TOPs (4 × 104 by snapshot) was found sufficient to reach convergence.

Transport Properties in Carbonated Hydroalcoholic Solutions

Although the determination of accurate diffusion coefficients in multicomponent liquids like champagnes might require the implementation of sophisticated methods,[52] the properties of these beverages (i.e., isotropy and homogeneity, no chemical reaction on the simulation timescale, and high abundance of water in the mixture) make the use of an effective formula similar to Fick’s law for binary liquids possible.[10,53] In particular, the probability density of diffusing CO2 molecules is Gaussian and their diffusion coefficient can be derived from the linear fitting of their mean-squared displacement (MSD) in a three-dimensional space since MSD(t) = 6Dt at long times. This approach is mathematically identical to integrating the velocity autocorrelation function (VACF) of the diffusing molecules,[54] a method often referred to as the Green–Kubo formula for diffusion. However, the value of CO2 diffusion coefficients is known to depend on the system size when periodic boundary conditions (PBC) are applied. Yeh and Hummer[55] proposed the following correction to self-diffusion coefficients to compensate for this shortcoming of PBCwhere D0 is the corrected self-diffusion coefficient, DPBC is the original self-diffusion coefficient extracted from MD simulations, T is the temperature, η is the shear viscosity of the solvent, L is the length of the cubic cell, and ξ is a constant equal to 2.837297. The shear viscosity can be evaluated by calculating the transverse-current autocorrelation functions (TCAF) of the liquid and subsequently fitting the wavenumber-dependent viscosities η(k) to the function η(k) = η(0) + ak2 where η(0), TCAF viscosity plotted in Figure , and a are real fitting parameters.[56] An alternate way to estimate the viscosity consists in applying the Stokes–Einstein formula that relates the dynamic viscosity of a liquid to the diffusion coefficient and hydrodynamic radius of a solvated species. In this work, the CO2 hydrodynamic radius is identified with the rms distance of CO2 atoms to their molecular center of mass (another definition based on the CO2 radius of gyration is provided in the Supporting Information). For whatever method considered to compute viscosity, CO2 diffusion coefficients only increase from 0.02 × 10–9 to 0.04 × 10–9 m2 s–1 when the corrections for system-size dependence of eq are applied to calculations with the recommended water models, namely, OPC and TIP4P/2005. Although the inclusion of such small corrections would not alter any conclusion of the present paper, we included them in the final comparison with experimental data for the sake of completeness (see Figure ). We have also checked that dividing the number of molecules by 4 in the simulation box and averaging diffusion coefficients over four trajectories do not drastically alter our results.

Water Properties

The relevance of water models to describe carbonated hydroalcoholic solutions is discussed in terms of the number of H bonds, tetrahedral order parameters (TOPs), and multipole moments. A H bond is assumed to occur when the Od–Oa distance between the donor and acceptor oxygens remains below 0.35 nm and the H–Od–Oa angle does not exceed 35°, as advocated by Chandler.[57] TOPs are evaluated by dividing them into an angular component Sg and a radial component Sk,[22] both components reaching zero for the perfect tetrahedral geometry. Other definitions exist that impose tetrahedrality for a value of 1[58] but this does not alter the interpretation of the results. Lower-order multipole moments (dipoles, quadrupoles, and octopoles) are derived from their general Cartesian expressions[59] by borrowing the notations proposed by Izadi et al. to build the OPC water model.[21] An rms error ΔξOPCi, where i refers to a water model, is introduced to evaluate the deviation between multipole moments corresponding to model i and those obtained with the OPC modelwhere μz is the dipole, Q0 and Q are two quadrupole components, and Ω0 and ΩT are two octopole components. If octopoles are not included in eq , the factor 1/5 is to be replaced by 1/3.