Mutual Diffusivities of Mixtures of Carbon Dioxide and Hydrogen and Their Solubilities in Brine: Insight from Molecular Simulations

H2-CO2 mixtures find wide-ranging applications, including their growing significance as synthetic fuels in the transportation industry, relevance in capture technologies for carbon capture and storage, occurrence in subsurface storage of hydrogen, and hydrogenation of carbon dioxide to form hydrocarbons and alcohols. Here, we focus on the thermodynamic properties of H2-CO2 mixtures pertinent to underground hydrogen storage in depleted gas reservoirs. Molecular dynamics simulations are used to compute mutual (Fick) diffusivities for a wide range of pressures (5 to 50 MPa), temperatures (323.15 to 423.15 K), and mixture compositions (hydrogen mole fraction from 0 to 1). At 5 MPa, the computed mutual diffusivities agree within 5% with the kinetic theory of Chapman and Enskog at 423.15 K, albeit exhibiting deviations of up to 25% between 323.15 and 373.15 K. Even at 50 MPa, kinetic theory predictions match computed diffusivities within 15% for mixtures comprising over 80% H2 due to the ideal-gas-like behavior. In mixtures with higher concentrations of CO2, the Moggridge correlation emerges as a dependable substitute for the kinetic theory. Specifically, when the CO2 content reaches 50%, the Moggridge correlation achieves predictions within 10% of the computed Fick diffusivities. Phase equilibria of ternary mixtures involving CO2-H2-NaCl were explored using Gibbs Ensemble (GE) simulations with the Continuous Fractional Component Monte Carlo (CFCMC) technique. The computed solubilities of CO2 and H2 in NaCl brine increased with the fugacity of the respective component but decreased with NaCl concentration (salting out effect). While the solubility of CO2 in NaCl brine decreased in the ternary system compared to the binary CO2-NaCl brine system, the solubility of H2 in NaCl brine increased less in the ternary system compared to the binary H2-NaCl brine system. The cooperative effect of H2-CO2 enhances the H2 solubility while suppressing the CO2 solubility. The water content in the gas phase was found to be intermediate between H2-NaCl brine and CO2-NaCl brine systems. Our findings have implications for hydrogen storage and chemical technologies dealing with CO2-H2 mixtures, particularly where experimental data are lacking, emphasizing the need for reliable thermodynamic data on H2-CO2 mixtures.


INTRODUCTION
−4 Given the intermittent production of renewable energy sources such as wind and solar, storing excess energy becomes essential to ensure consistent supply. 5−4 At present, worldwide hydrogen consumption stands at approximately 115 million metric tons (115 × 10 9 kg) per year 3 and is expected to increase 10fold (to the gigatonnes per year scale) to accommodate the expanding global energy needs. 3Storing a gigatonne of hydrogen would necessitate immense storage volumes of the order of 10 11 m 3 , even at significant pressures of 10 MPa and temperatures of 323.15 K. 3 Only geological formations such as depleted gas reservoirs and salt caverns can offer such large volumes for hydrogen storage. 33 Processes such as the mixing of hydrogen with the cushion gas via molecular diffusion, and its dissolution into brine can decrease the purity of the stored hydrogen.1,3,6 Thus, estimating the mutual diffusion coefficients characterizing gas mixing and the solubility of the resulting gas mixtures in brine is essential for enabling underground storage of hydrogen.2,3 Oldenburg and colleagues 7,8 showed that using carbon dioxide (CO 2 ) as a cushion gas facilitates efficient injection and production of natural gas due to its high compressibility at reservoir-relevant thermodynamic conditions.CO 2 is opted as the cushion gas in our study due to the similarities between hydrogen and natural gas storage. 4 .1. Fic Diffusivities of CO 2 -H 2 Mixtures.Mass transport in binary gas mixtures, in the absence of temperature and pressure gradients, is described by Fick's law, which linearly relates the diffusive flux to the concentration gradient of a species.9−11 An alternative description of the same phenomenon can be formulated by the Maxwell−Stefan (MS) framework, 12 which describes diffusive mass transport based on chemical-potential gradients.The proportionality constant in Fick's law, denoted by the mutual or Fick diffusion coefficient, 10,11 can exhibit a strong dependence on the temperature, pressure, and mixture composition.10,11 Analytic expressions, experimental measurements, semiempirical methods guided by theory, and numerical simulations have been used to calculate mutual diffusivities in binary gas mixtures.11,13−15 Chapman and Enskog 16−18 independently derived an analytic expression to predict the mutual diffusion coefficients of binary gas mixtures by solving the Boltzmann transport equation in the dilute gas limit in a framework commonly referred to as the kinetic theory of gases.14,16−18 Kinetic theory assumes gases are made of monatomic, spherical, nonpolar particles that undergo elastic collisions, i.e., there are no strong molecular interactions between them. The pediction of diffusivities by Chapman and Enskog has been shown to be highly accurate in the dilute gas regime.19−21 At moderate to high pressures (less than tens of MPa 11 ), the applicability of kinetic theory for dilute gases is in doubt due to its simplistic assumptions of intermolecular interactions.11,17,18 Many researchers have expanded kinetic theory to address scenarios involving dense gases, 14,17,22 notably the work of van Beijeren and Ernst, who derived diffusion coefficients and other properties for moderately dense, multicomponent mixtures of hard spheres, known as the revised Enskog theory (RET).22 Jervell and Wilhelmsen 23 recently presented an RET for multicomponent mixtures of Mie fluids (RET-Mie), which is claimed to accurately predict transport properties such as diffusion coefficients, thermal diffusivities, viscosities, and thermal conductivities.With parameters of the Mie-potentials fitted to equilibrium properties of the respective components, predictions of mutual diffusivities of noble gases for pressures ca.0.1 MPa for a wide temperature range (200−1400 K) are within 4% of the experimental data.24,25 RET-Mie predictions of self-diffusivities of H 2 are within 10% of the available experimental data up to 200 MPa and at temperatures between 171 and 372 K, and mutual diffusivities of noble gases are predicted within 4% of the experimentally determined values.Fick diffusivities of carbon dioxide-methane (CH 4 ) mixtures are within 20% of experimental measurements and molecular simulations for pressures up to 14.7 MPa. 26 Except for xenon in the vicinity of the critical point, the thermal conductivity of noble gases and their mixtures is reproduced within 10% of the experimental data at 300 K and pressures between 0.01 and 10 MPa.RET-Mie predictions of the viscosity are within ca. 10%of the experimental data for methane, nitrogen, and argon up to 300 bar for temperatures ranging between 233 and 523 K.At pressures up to 500 bar and temperatures from 200 to 800 K, the predictions are within ca.15% of the most accurate correlation for the viscosity of air.Larger inaccuracies in the predictions of mutual diffusivities of CO 2 -CH 4 mixtures reflect the inability of the theory to account for nonspherical shapes of molecules.In conclusion, kinetic theory is incapable of predicting the mutual diffusivities of CO 2 -H 2 mixtures at reservoir-relevant thermodynamic conditions due to its simplistic assumptions.Advanced theories such as the RET-Mie (or RET) are cumbersome to use and still require significant efforts to account for nonspherical molecules.14,23 Marrero and Mason 19 and Li et al. 27 provide an extensive review of experimental data on the mutual diffusivities in binary mixtures of CO 2 and H 2 , a summary of which can be found in Table 1.This table highlights the extensive literature containing data for a wide range of temperatures, albeit confined to pressures near the atmospheric pressure (0.1 Loschmidt 141 273 and 286 closed tube Loschmidt 142 273−289 closed tube Schmidt 143 288 closed tube Lonius 144 286−294 closed tube Boardman and Wild 145 288 closed tube Suetin 146 273 closed tube Wrestschko 297 closed tube Suetin and Ivakin 147 292 closed tube Ivakin and Suetin 148 292−473 closed tube von Obermayer 149 285−335 open tube von Obermayer 150 284−293 open tube von Obermayer 151 280−294 open tube Schafer et al. 152 252−308 two-bulb Lonsdale and Mason 153 259−358 two bulb MPa).Chou and Martin 28 measured the self-diffusion coefficients of infinitely diluted Carbon-14 labeled C 14 13 To the best of our knowledge, experimental data for the Fick diffusivities of binary mixtures of CO 2 and H 2 at pressures exceeding the atmospheric pressure have not been reported.
To facilitate rapid and precise predictions of Fick diffusivities in gas mixtures, numerous semiempirical expressions guided by the kinetic theory have been developed. 11,13These expressions aim to fill the gap resulting from the scarcity of experimental measurements, specifically for predicting binary diffusion coefficients in both dilute and dense gases. 13Correlations by Wilke-Lee 29 and Fuller et al. 30−32 are often used to predict diffusion coefficients in dilute gas mixtures.These correlations fail at high pressures, 11,13 just as Chapman and Enskog's analytic expression for mutual diffusivities, and cannot take into account the compositional dependence of the mixture.Several methods have been proposed to estimate diffusion coefficients in dense gas mixtures. 13Takahashi and Hongo 33 and Riazi and Whitson 34 suggested empirical equations based on the method of corresponding states to correlate the product of pressure and diffusion coefficient at low and high pressures. 13Takahashi and Hongo's correlation for the effect of pressure and temperature on the binary diffusion coefficient requires the evaluation of an empirical equation using look-up tables that are tabulated as functions of the critical temperature and pressure of the components in the mixture. 13Riazi and Whitson's correlation requires the viscosity of the mixture as an input and does not yield the correct value of mutual diffusivities in the limit of low pressure. 13No correlation is entirely satisfactory for predicting mutual diffusivities at both low and high pressures, primarily due to the simplistic assumptions of the theory of corresponding states. 13Although semiempirical methods are relatively easy to use and computationally fast, their accuracy depends on the extent and quality of the experimental data that was used for their development.Generally, one needs to exercise caution when using semiempirical methods outside the range used for its calibration.
Predictions of Fick diffusivities for binary mixtures of CO 2 -H 2 using semiempirical methods can be accurate between low to intermediate pressures but can be inaccurate due to neglecting strong molecular interactions emanating at large pressures (≳10 MPa).Molecular dynamics (MD) simulations, a potent tool in materials research, 35,36 operate at atomistic length and time scales, offering insights at the molecularscale. 35,36The formulation of accurate force fields for multiple components, the design of efficient algorithms, coupled with advancements in computational power, has made MD a go-to method for calculating the diffusivities of pure components, mixtures, and other transport properties of fluids. 35,36−42 In sharp contrast, the Fick diffusivities of binary gas mixtures computed from MD simulations have not received much attention.Vella 43 computed Fick diffusivities for CO 2 -CH 4 mixtures using MD simulations at temperatures and pressures corresponding to the dilute and dense regimes.The computed diffusivities (within statistical uncertainties), described faithfully by Fuller's correlation, 30−32 agreed with the experimental data between 293 and 500 K and 1.01325 MPa (dilute gas regime).Deviations exceeding 40% were observed between Fuller's correlation and the computed diffusivities at ca. 310 K and pressures above ca. 2 MPa (the dense gas regime).Due to the lack of experimental data in the dense gas regime, Vella 43 was unable to confirm the accuracy of the Fick diffusivities computed by MD simulations.Guevara-Carrion 26 compared the measured Fick diffusivities of very dilute CH 4 in supercritical CO 2 at temperatures of 292.55 and 332.85K and pressures of 9, 12.5, and 14.7 MPa.The mutual diffusivities changed from liquid-like to gas-like values within a small temperature range near the supercritical region across the so-called Widom line. 26Saric et al. 44 measured the mutual diffusivities of binary mixtures of H 2 diluted in CO 2 , along with other binary gas mixtures.Three mole fractions of H 2 , representing the infinite-dilution limit (mostly 0.5, 1.0, and 1.5 mol %), were examined between 290 and 350 K at 9 MPa.Saric et al. 44 proposed a simple equation to predict the Widom line of dilute supercritical CO 2 mixtures as a function of pressure, solute mole fraction, and critical properties of the mixture components.In a separate study, Saric et al. 45 examined the mutual diffusion of various hydrocarbons dispersed in supercritical carbon dioxide via MD simulations near the Widom line, between 290 and 345 K at 9 MPa.Fick diffusivities were strongly affected close to the critical region of the mixture due to the pronounced clustering of solute molecules.Saric et al. 45 pointed out that the validity of the Stokes−Einstein relation near the Widom line was questionable.A correlation function was proposed to facilitate quick calculations of Fick diffusivities for various binary mixtures.
Rosenfeld 46 proposed an empirical correlation connecting transport properties to the excess entropy of a fluid at a given temperature and pressure.−51 Despite its nonrigorous nature, the concept of excess entropy could serve as a valuable tool in connecting thermodynamics and dynamics for predicting transport properties in mixtures. 47,48o summarize, neither experimental data nor MD simulations have reported the Fick diffusivities of CO 2 -H 2 mixtures under reservoir conditions.Analytic expressions can be cumbersome and require effort to include nonspherical molecules. 14,17,23Semiempirical methods, though easy to use, are not accurate for the entire range of thermodynamic conditions of interest.Understanding the dissolution of H 2 in saline water, referred to as brine, at various thermodynamic conditions is necessary for estimating hydrogen purity during subsurface storage.This involves predicting solubilities by assessing equilibrium between vapor and liquid phases in a multicomponent system at different temperatures and pressures, known as vapor−liquid equilibria (VLE). 13,52,53

Industrial & Engineering Chemistry Research
Studies by Spycher et al. 54,55 and Mutailipu provide an overview of the VLE of various CO 2 -brine systems compiled from experimental studies.Spycher et al. 54 compared data on mutual solubilities in CO 2 -H 2 O mixtures from various experiments for pressures up to 50 MPa and temperatures below 373.15 K. Solubilities measured in the two-phase region, where a CO 2 -rich phase (usually gas) coexists with a H 2 O-rich liquid phase, were consistent across various studies. 54A thermodynamic model was developed by Spycher et al. 54 based on the Redlich−Kwong equation-of-state to enable quick predictions of mutual solubilities in CO 2 -H 2 O systems.This model predicts mutual solubilities with an accuracy of a few percent of the experimental data for pressures up to 60 MPa and a temperature ca.383.15 K.In a subsequent study, Spycher et al. 55 extended the thermodynamic model developed initially for CO 2 -H 2 O mixtures 54 to include the influence of chloride salts dissolved in the aqueous phase.The extended thermodynamic model to predict mutual solubilities included activity coefficients for aqueous CO 2 derived from studies developed primarily for sodium chloride (NaCl).The activity coefficients of NaCl provided by Duan and Sun 56 and Rumpf et al. 68 were shown to yield CO 2 solubilities in NaCl brine that were closest to the experiments by Duan and Sun, 56 Drummond, 69 Malinin and Saveleva, 70 Malinin and Kurovskaya, 71 Rumpf et al., 68 Nighswander et al., 72 Bando et al., 73 and Kiepe et al., 74 between pressures 3.5 and 6.8 MPa and between 293.15 K and 373.15 K.The extended solubility model provided by Spycher and Pruess, 55 including the effect of salts, was shown to be accurate to within the experimental uncertainty for solutions up to 6 molal (= moles of salt per kg of water) for NaCl and 4 molal for CaCl 2 .Spycher and Pruess 75 extended their solubility model 55 to cover temperatures up to 573 K (from 285 K) and pressures below 60 MPa.
Molecular simulations have been extensively applied to investigate VLE in CO 2 -H 2 O and CO 2 -NaCl brine systems.Vorholz et al. 63 studied the VLE of H 2 O between 323 and 573 K, CO 2 between 230 and 290 K, and CO 2 -H 2 O mixtures between 348 and 393 K using NVT and NPT versions of Gibbs Ensemble Monte Carlo 35,36 (GEMC) simulations.Three different H 2 O force fields, namely SPC, SPC-E, and TIP4P, the EPM2 force field for CO 2 , and combinations thereof were examined to find the optimal VLE description for binary mixtures of CO 2 -H 2 O.The best prediction for the VLE of CO 2 -H 2 O binary mixtures, up to pressures of 20 MPa, was achieved by combining the SPC and TIP4P force fields for water with the EPM2 force field for carbon dioxide.Liu et al. 60 used histogram-reweighting Grand-Canonical Monte Carlo (GCMC) 35,36 simulations to obtain the phase behavior of CO 2 -H 2 O mixtures for 323.15 and 723.15 K, and 0 and 100 MPa.Several different force fields for H 2 O and CO 2 were tested using the conventional Lorentz−Berthelot mixing rules 35,36 for the interaction between unlike atoms.Experimentally measured solubilities deviated from solubilities computed using force field based simulations for the entire range of temperatures and pressures spanning 323.15 to 723.15 K and 0 to 100 MPa, respectively.The work of Liu et al. 61 suggests the need for improved atomistic force fields to accurately predict the VLE of CO 2 -H 2 O systems.Liu et al. 61 investigated the phase behavior and interfacial tension of CO 2 -H 2 O-NaCl mixtures using MD simulations between 323.15 and 523.15 K, 0−60 MPa, and NaCl concentrations 1−4 molal.TraPPE 76,77 /TIP4P2005 78 performed the most effectively among various combinations of force fields evaluated for predicting the VLE of the H 2 O-CO 2 -NaCl system.The dependence of the interfacial tension on temperature, pressure, and salt concentration was also examined, revealing insights into the predictive capabilities of different force field combinations.Orozco et al. 62 conducted NPT GEMC simulations to obtain optimized intermolecular potential parameters to describe the phase behavior of CO 2 -H 2 O mixtures.The simulations covered a temperature range of 423−523 K and pressures from 20 to 80 MPa, addressing conditions pertinent to carbon capture and sequestration processes.Lobanova et al. 64 57,59 Ansari, 58 and Tori ́n-Ollarves and Trusler. 79Recognizing the lack of data for pressures exceeding 0.1 MPa, Chabab et al. 57 conducted experiments to measure the solubility of H 2 in brine, specifically NaCl brine, with molality ranging from 0 to 6 molal.Experiments were conducted between 323 and 373 K and 0.1−23 MPa.To enable quick and accurate predictions of H 2 solubilities in brine as functions of temperature, pressure, and molality, the e-PR-CPA  57 and Tori ́n-Ollarves and Trusler. 79The H 2 solubility data were used to parametrize and assess three different thermodynamic models: (1) the e-NRTL activity coefficient model, 86,87 (2) the Duan-type (noniterative) model, 88,89 and (3) the Søreide and Whitson's equation of state. 90These models were shown to be highly accurate in predicting the phase equilibria of H 2 -NaCl Industrial & Engineering Chemistry Research brine systems, with an average absolute deviation of less than 3%. 59 2 solubilities in H 2 O and NaCl brine have also been studied using molecular simulations.66,91 Rahbari et al. 66 reported H 2 solubilities in H 2 O calculated from Monte Carlo simulations for pressures between 1 and 100 MPa and temperatures between 283 and 423 K.The H 2 solubilities calculated using molecular simulations were in excellent agreement with the experimental data at high pressures by Wiebe and Gaddy. 92 Inanother article, Rahbari et al. 93 studied the VLE of hydrogen-water mixtures along the melting line of ice Ih, corresponding to temperatures between 264.21 and 272.4 K. Based on the computed solubility data of H 2 in H 2 O, the freezing point depression of water caused by NaCl was calculated for temperatures between 264.21 and 272.4 K. Lopez-Lazaro et al. 94 computed the solubilities of H 2 in aqueous NaCl solutions using molecular simulations for molalities up to a maximum of 2 molal between temperatures of 280 and 350 K. Recently, van Rooijen et al. 95 reported H 2 solubilities in aqueous NaCl solutions for temperatures between 298 and 363 K, pressures ranging between 0.1 and 100 MPa, and molalities between 0 and 6 molal, together with other thermodynamic properties such as interfacial tensions and transport properties of H 2 -H 2 O-NaCl systems. ThH 2 concentrations in NaCl brine, calculated using molecular simulations, 95 are larger compared to the predictions by Chabab et al. 57 and lower than the prediction by Tori ́n-Ollarves and Trusler 79 for NaCl molalities larger than 0.5 molal.Despite numerous studies exploring solubilities in CO 2 -NaCl brine and H 2 -NaCl brine systems for diverse temperatures and pressures, there is a lack of data regarding solubilities in CO 2 -H 2 -NaCl brine systems.
Recognizing the lack of information on binary Fick diffusivities in CO 2 -H 2 mixtures at high pressures, temperatures, and various compositions and the limited understanding of the VLE in ternary systems of CO 2 -H 2 -NaCl brine, molecular simulations are performed to address this gap.This study investigates the Fick diffusion coefficients computed from equilibrium MD simulations of H 2 -CO 2 mixtures.−100 Temperatures between 323.15 and 423.15K and pressures between 5 and 50 MPa are explored.We only examine ternary mixtures where the overall mole fractions of CO 2 , H 2 , and Brine, calculated using the total number of molecules in both the gas and liquid phases, are 1/6, 1/6, and 2/3, respectively.This paper is organized as follows: Section 2 provides details on the force fields and molecular simulation techniques.In Section 3, the self-diffusivities of H 2 and CO 2 are compared, along with the mutual diffusivities and solubilities obtained via molecular simulations.Comparisons to available experimental data and analytic expressions are included.Concluding remarks are presented in Section 4. All computed raw data from this study are tabulated in the Supporting Information.

Force Fields.
The TraPPE force field is used to simulate CO 2 as a rigid linear molecule. 76,77TraPPE compares favorably with experimental data for the VLE of pure CO 2 and its multicomponent mixtures for a wide range of temperatures, pressures, and compositions. 76,77The three-site Marx force field 101,102 is used to simulate H 2 as a rigid body for its accuracy in reproducing the bulk densities and fugacities of H 2 at pressures up to 100 MPa. 66,103Quantum effects emanating at low temperatures are insignificant for the temperatures considered in this work (323.15−423.15K). 104,105 Water is simulated as a rigid body using the TIP4P/μ force field, a modified version of the TIP4P/2005 model 78 developed by Rahbari and colleagues. 66The TIP4P/μ force field accurately predicts the VLE of H 2 O-H 2 mixtures between 323 and 423 K and 10−100 MPa 66 as it was partially trained to reproduce the excess chemical potential of liquid water.The Madrid-2019 force field 106 is used to model the Na + and Cl − ions, which are parametrized for the TIP4P/2005 H 2 O force field. 78Lennard-Jones (LJ) and electrostatic interactions are calculated for all species.All interaction parameters used in this study are listed in Table S1 of the Supporting Information.The Lorentz− Berthelot (LB) mixing rules 35,36 are applied in most cases barring a few exceptions, see Table S1 of the Supporting Information.
2.2.Molecular Dynamics.MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 107 (version 23 June 2022) to calculate the self-and Fick diffusivities in CO 2 -H 2 mixtures.Pressures between 5 and 50 MPa in 5 MPa increments, temperatures between 323.15 and 423.15K in 25 K increments, and hydrogen mole fractions between 0 and 1 in 0.1 increments were explored.At these thermodynamic conditions, H 2 -CO 2 mixtures exhibit either gas-or dense-gaslike behavior. 108Simulation boxes are periodic in all three directions.Equations of motion are integrated using the velocity-Verlet algorithm using a time step of 0.5 fs.Lennard-Jones and long-range electrostatic molecular interactions are calculated with a cutoff radius of 10 Å, supplemented by tail corrections for energy and pressure.All long-range electrostatic interactions are computed using the Particle−Particle Particle-Mesh method with a relative precision of 10 −6 .Neighbor lists for every molecule are generated each time step using a binning procedure by considering a radius equal to the force cutoff of 10 Å plus an additional "skin" of thickness 3 Å.
Transport properties such as viscosities, self-, MS, and Fick diffusivities are computed in the NVE ensemble using systems of 120 molecules.Small systems of 120 molecules are chosen to manage the high computational demands associated with the gas phase simulations for the calculation of (self-and MS) diffusivities for trajectories typically lasting 100 ns.Accurate computation of transport properties requires correction for finite-size effects, a topic addressed in Section 2.2.1.To ensure that simulations in the NVE ensemble yield the correct average pressure and temperature, the box sizes and total system energy (in kJ•mol −1 ) are determined from independent simulations of 2000 molecules.Small systems of 120 molecules are unsuitable for computations of densities due to the possibility of the interaction cutoff radii exceeding half of the box size, violating the nearest image convention.Long simulations in the NVT phase are necessary to obtain an average temperature due to the substantial fluctuations

Industrial & Engineering Chemistry Research
resulting from the small system size.The procedure and the raw data of the simulations for the larger systems are described and tabulated in Section S1.1 and Table S2 of the Supporting Information.

Calculation of Transport Properties.
Viscosities describe momentum transport, 109 while self-, MS, and Fick diffusivities elucidate mass transport. 10,109,110Self-diffusion of a molecule in a fluid is caused by Brownian motion and occurs irrespective of gradients in chemical potential, temperature, or pressure. 110In contrast, MS and Fick diffusion 10,109 describe the collective motion of molecules in response to gradients in chemical potential and concentration, respectively.Self-, MS, and Fick diffusion all contribute to mass transport via random or directed molecular motion, as explained in various texts. 10,109,110Viscosities, self-, MS, and Fick diffusivities are calculated in the NVE ensemble.The simulation box size corresponds to the density determined from systems of 2000 molecules, and the (scaled) total energy is imposed to obtain the desired temperature.Simulations consist of an initialization phase lasting 1 ns to achieve equilibration, followed by a production phase lasting between ca.30−150 ns.The duration of the production phase varies widely due to density variations, affecting computational demands.Simulations at higher densities (e.g., systems at 50 MPa and 323.15 K) are computationally cheaper than those at lower densities (e.g., systems at 5 MPa and 423.15 K).During the production phase, the OCTP plugin 111 in LAMMPS 107 is used to compute the shear viscosities and self-and MS diffusivities based on meansquared displacements.0][41][42]112,113 Recently, Vlugt and co-workers 114 showed that the thermodynamic factors (Γ) computed using the CFCMC technique [96][97][98][99][100]115 for various CO 2 -H 2 mixtures at 323.15 K and pressures of 5 and 50 MPa agreed within 3% of the REFPROP data. 108 Thu, for convenience and efficiency, thermodynamic factors for H 2 -CO 2 mixtures are computed by numerically differentiating the Margules activity coefficient model based on the excess Gibbs energies obtained from REFPROP, 108 as discussed in Section S5.To obtain accurate values of transport properties, an average value is calculated over 5 independent simulations at each thermodynamic state point. Uncrtainties in the computed transport properties are calculated as the standard deviation across these 5 independent simulations.
While viscosities are usually unaffected by finite-size effects, 116,117 self-, MS, and Fick diffusivities can depend on the system size, 116,118,119 particularly in dense systems. 117The finite-size correction to the self-diffusion coefficient proposed by Yeh and Hummer 118,120 is where D i self is the finite-size corrected self-diffusivity of the ith species, ξ is a dimensionless constant equal to 2.837298, and L is the box length of a cubic simulation box.Jamali et al. 116 showed that the MS diffusion coefficients can be corrected according to where D̵ MS is the finite-size corrected MS diffusivity, D YH is the Yeh−Hummer correction term in eq 1, 118 and Γ is the thermodynamic factor for diffusion.For ideal mixtures where cross-correlations between molecular displacements of different species are negligible, the MS diffusivities follow the Darken equation 11,116,121,122 Taking advantage of the relation between D̵ MS and D Fick , 12 Jamali et al. 116 showed that the correction for the Fick diffusivity follows 116,118

Monte Carlo Simulations.
VLEs of CO 2 -NaCl brine mixtures, H 2 -NaCl brine mixtures, and CO 2 -H 2 -NaCl brine mixtures are determined via MC simulations performed using the open-source Brick-CFCMC software. 98,100The focus is on thermodynamic conditions relevant to subsurface storage of H 2 : pressures between 5 and 50 MPa, temperatures between 323.15 and 423.15 K, and NaCl concentrations varying between 0 and 2 molal.The reaction equilibria of CO 2 -H 2 O systems are not considered as bicarbonate concentrations are low. 1236][97][98][99][100]115 Simulations are performed using the NPT version of the CFCGE technique, developed by Poursaeidesfahani et al., 115 to determine the VLEs of mixtures consisting of H 2 , CO 2 , and NaCl brine.A single fractional molecule per component i, where i can be either CO 2 , H 2 O, or H 2 , but not NaCl, is introduced in the simulations that is distinguishable from the whole molecules.115 Fractional molecules can be located in either of the two simulation boxes.The Lennard-Jones and electrostatic interactions of the fractional molecule with surrounding molecules are modulated by the coupling parameter λ i ∈ [0, 1].At λ i = 0, the interactions resemble those of an ideal gas, while at λ i = 1, interactions are fully scaled like those of a whole molecule.For the exact expressions of the scaled interactions the reader is referred to articles by Hens et al. 98 and Polat et al. 100 The details of the trial moves relevant for the CFCGE method: (1) molecule translations, (2) molecule rotations, (3) volume changes, and (4) hybrid insertions/deletions, are provided in the articles by Hens et al., 98 Polat et al., 100 and Rahbari et al. 99 To ensure that all values of λ i are sampled in both simulation boxes, a weight function or a biasing function W(λ i , j) is introduced for species i in box j, where j = 1 or 2. 98,100,115 The excess chemical potentials of each species i in box j in the NPT ensemble (μ i,j Ex ) equals 66,93,98,100,115 where p j (λ i ) is the probability that a fractional molecule of species i samples a value of λ i in box j, k B is the Boltzmann constant, and T is the absolute temperature.In addition to assessing equilibrium between systems, the excess chemical potential aids in determining fugacity coefficients, crucial for describing phase-coexistence in mixtures, 52,53 particularly used in equation of state-based thermodynamic modeling for phase Industrial & Engineering Chemistry Research equilibria.Fugacity coefficients of species i (ϕ i ) are directly related to its excess chemical potential in the NPT ensemble as 66,93 = i k j j j j j j y is the compressibility factor for a mixture in the gas phase whose average volume equals ⟨V⟩.For further details of the CFC method the reader is referred to the articles by Hens et al. 98 and Polat et al. 100 In the systems of H 2 -NaCl brine and CO 2 -NaCl brine, 200 molecules of H 2 (or CO 2 ) in the gas-rich phase and 400 molecules of H 2 O in the liquid-rich phase are initialized.For the ternary system of H 2 -CO 2 -NaCl, 100 molecules each of H 2 and CO 2 are initialized in the gas-rich phase, and 400 molecules of H 2 O is initialized in the liquid-rich phase.0, 7, and 15 ions of Na + and Cl − , corresponding to NaCl concentrations of 0, 1, and 2 molal, are distributed randomly inside the H 2 O-rich simulation box.A cutoff radius of 10 Å is applied for the LJ interactions, and analytic tail corrections are applied.Electrostatic energies are computed using the Ewald summation 124 with a real-space cutoff of 10 Å and a damping parameter α = 0.32 Å −1 for both simulation boxes.The number of k-vectors in Fourier space is k = 8 in the liquid-rich box.To ensure convergence of electrostatic energies in the expansive gas-phase boxes, the maximum k-vectors are varied between 9 and 19.−100 These settings ensure that the electrostatic energies are computed with a relative precision of 10 −6 .
Every cycle of a CFCGE simulation contains N Total Monte Carlo (MC) trial moves, where N Total is the total number of molecules.In the binary and ternary systems, molecule translations, rotations, and volume changes are selected with probabilities of 0.35, 0.24, and 0.01, respectively.Volume changes are computationally expensive since all configurational energies need recalculation based on the updated intermolecular distances after rescaling the simulation box.Consequently, volume changes are typically executed with a probability of 0.01. 98,100Changes of λ i and molecule transfer moves (swapping and identity changes) are selected with probabilities of 0.2 each.The swapping of a fractional molecule  S2 and S3 of the Supporting Information, respectively.
of species i between boxes is attempted if λ i < 0.3, and the identity change of the fractional molecules is attempted only if λ i > 0.7. 115The maximum displacements, rotations, volume changes, and λ changes are adjusted to obtain acceptance of ca.50%.Molecule transfers involving salt ions Na + and Cl − are excluded.For all simulations, an initialization phase lasting 10 4 cycles is simulated to remove overlaps between molecules.In an equilibration phase lasting 15 × 10 6 cycles, the weight functions W(λ i , j) are developed in each box, and a production phase lasting 10 6 cycles is used to calculate solubilities and fugacity coefficients (see eq 6).To obtain accurate values of solubilities and fugacity coefficients, 40 independent simulations are performed at each thermodynamic state point.The uncertainties (error bars) in the computation of solubilities and fugacity coefficients are estimated by dividing the 40 independent simulations into blocks of 5. 36 The solubilities of H 2 and CO 2 are determined as the ratio of the number of gas molecules to the total molecules in the liquid phase, with sodium and chloride ions considered separately.Liquid solubility in the gas phase is determined by the fraction of water molecules in the gas-rich phase.

Densities and Compressibilities of CO 2 -H 2
Mixtures.H 2 -CO 2 mixture densities (ρ) at 5, 25, and 50 MPa are plotted as a function of the mole fraction of H 2 (x Hd 2 ) in Figure 1a,c,e.The computed densities are compared with predictions from equations of state found in REFPROP. 108learly, densities decrease monotonically and nonlinearly with the hydrogen content in the mixture.In Figure 1, densities computed from MD simulations closely match predictions from REFPROP, 108 typically within 5% at all thermodynamic state points.The relative deviations between the computed and predicted densities are shown in Figures S5−S7 of the Supporting Information.
The ratio of the densities of pure CO 2 and H 2 can vary between 20 and 35, considering all pressures and temperatures.The large density contrast between the pure components also impacts the mixture densities, as seen by the positive curvature in Figure 1a,c,e, suggesting a strong influence of CO 2 content on the mixture densities.As expected, and as evident from Figure 1a,c,e, densities at a fixed pressure and mixture composition decreases with temperature.A ca. 31% increase in temperature (323.15 to 423.15 K), results in ca.30% and ca.40% reductions of densities of pure H 2 and CO 2 , respectively, at 5 and 50 MPa.Variations in H 2 -CO 2 mixture densities with temperature lie between the bounds set by the pure components.It can be concluded that density and temperature have a linear-like relationship.A closer look at the vertical axis in Figure 1a,c,,e, indicates that the density also increases with pressure.Barring the mixture at x Hd 2 = 0 and T = 323.15K, ρ(p) can be fitted to a second-degree polynomial with an error of up to ca. 10%.At x Hd 2 = 0 and T = 323.15K, densities sharply increase until p ≈ 10 MPa, followed by a more gradual rise.The proximity of this specific mixture to the critical point of CO 2 (T C = 304.15K, p C = 7.38 MPa) 125 could explain the transition from gas-like behavior (p ≲ 10 MPa) to fluid-like behavior (p ≳ 10 MPa).The reader is referred to Figures S1−S4 of the Supporting Information for additional plots of densities as functions of mole fraction of H 2 and pressure.
The compressibility factors (Z) of H 2 -CO 2 mixtures are investigated to understand intermolecular interactions.Compressibilities vary nonlinearly with x Hd 2 , showing values both above and below unity.Excellent agreement between simulations and REFPROP predictions is observed for various pressures and temperatures.In general, mixtures rich in CO 2 exhibit attractive forces (Z < 1), while H 2 -enriched mixtures display repulsive forces (Z > 1).Temperature affects compressibilities differently at various pressures.Compressibilities tend to approach unity at low pressures and high temperatures.The behavior of compressibilities becomes more complex at higher pressures.Plots of compressibilities as a function of pressure and mole fraction of H 2 are provided in Figures S8−S10 of the Supporting Information.
To summarize, CO 2 -H 2 mixture densities are dominated by the CO 2 content in the mixture.In mixtures with a high proportion of H 2 (x Hd 2 ≳ 0.7), repulsive interactions dominate, leading to densities lower than those of an ideal gas.Conversely, CO 2 -rich mixtures have higher densities due to attractive forces compared to ideal gas densities.Gas-and liquid-like behavior of H 2 -CO 2 mixtures is observed depending on the pressure, temperature, and mixture composition.The densities and compressibilities computed from MD simulations and REFPROP 108 for all mixtures are tabulated in Table S2 of the Supporting Information.
3.2.Viscosities of CO 2 -H 2 Mixtures.Viscosities (η) are usually independent of finite-size effects. 116,118The viscosities of various CO 2 -H 2 mixtures as a function of x Hd 2 for 3 different pressures are shown in Figure 1b,d,f, while the remaining plots are displayed in Figures S12 and S13 of the Supporting Information.Similar to the mixture densities, the viscosities of H 2 -CO 2 mixtures also show gas-and liquid-like behavior.Gaslike behavior is evident at 5 MPa, as depicted by viscosities within the range of a few μPas in Figure 1b. 100 Regardless of temperature, viscosities decrease monotonically with x Hd 2 .As the mixture becomes enriched with H 2 , the viscosities transition toward a more gas-like behavior, resulting in a gradual decrease of approximately 20% until x Hd 2 ≲ 0.8, followed by a steeper drop of about 40%.By increasing the pressure to 25 and 50 MPa, specifically below 373.15 K, the behavior of viscosities differs from that shown in Figure 1b, indicating a liquid-like behavior.For example, at 25 MPa and 323.15 K, the viscosity decreases sharply by ca.67% from about 75 μPas for pure CO 2 to around 25 μPas at x Hd 2 = 0.5.Subsequently, viscosities decrease nearly linearly with the mole fraction of H 2 .At 25 MPa and 423.15 K, η(x Hd 2 ) resembles gaslike behavior.Viscosities are significantly impacted by the H 2 concentration in the H 2 -CO 2 mixture, as indicated by the negative curvature of η(x Hd 2 ) in Figure 1d, . The variation of viscosities with temperature further exemplifies the gas-and liquid-like behavior of the H 2 -CO 2 mixtures.At 5 MPa and any fixed x Hd 2 , viscosities increase nearly linearly with temperature.An increase from 323.15 to 423.15 K results in a ca.20% rise in the viscosities.The increase in viscosities with temperature reflects a gas-like behavior wherein the kinetic component of the viscosity exceeds the virial contribution due to molecular interactions, signifying the escalating influence of thermal motion on momentum transport in molecules.In contrast to gases, liquids experience a drop in viscosity with temperature. 13At 25 MPa in Figure 1d, for x Hd 2 ≲ 0.4, viscosities decrease with temperature.While the same behavior is observed at 50 MPa, the liquid-like behavior is prominent for x Hd 2 ≲ 0.55 in Figure 1f.Above a crossover concentration, x Hd 2 ≈ 0.4 at 25 MPa and x Hd 2 ≈ 0.55 at 50 MPa, viscosities increase with temperature due to the gas-like behavior of H 2 -enriched mixtures.The crossover concentration signals the transition from a liquid-like (or dense-gas like) behavior in CO 2 -enriched mixtures to a gas-like behavior in H 2 -enriched mixtures.
The viscosities can change significantly with pressure.Pure CO 2 shows a 4 to 5 times increase in viscosity upon increasing the pressure between 5 and 50 MPa, typical of dense gas-like behavior. 11The addition of H 2 to the mixture substantially reduces viscosity.Pure H 2 maintains relatively low viscosities, with an insignificant dependence on pressure.The transition from gas-like to liquid-like behavior typically occurs within a pressure range of 30−40 MPa.These observations suggest that viscosity behavior is highly dependent on both pressure and composition, with clear distinctions between gas-like and liquid-like regimes.Plots for the viscosities as a function of pressure for various mole fractions of H 2 are shown in Figures S14 and S15 of the Supporting Information.
Computed viscosities agree (on average) within 5% with REFPROP, 108 barring a few systems at 5 and 10 MPa.Deviations up to 28% are observed, noticeable from the large uncertainties at 5 MPa in Figure 1b, for instance, the data point at 373.15 K.The sparse intermolecular interactions in the gas phase lead to large fluctuations in p 111 and require up to 40 ns to accurately compute viscosities.Another way to improve the accuracy of the computed viscosities is through more independent simulations.We refrained from conducting additional or longer simulations due to the satisfactory accuracy of the computed viscosities and the substantial computational costs associated with simulating the gas phase using MD.For instance, on a single core of an AMD7H12 2.6 GHz processor, simulating 1 ns of pure H 2 at 5 MPa and 423.15K requires ca. 15 CPU hours, while pure CO 2 at 50 MPa and 323.15K requires roughly 1 CPU hour.Computed and predicted values of viscosities for H 2 -CO 2 mixtures are tabulated in Table S3 of the Supporting Information.S3 of the Supporting Information.
be corrected for finite-size effects, 116−119 depending on the density of the system.Self-diffusivities are computed for systems containing 120, 250, 500, 1000, and 2000 molecules for 2 equimolar H 2 -CO 2 mixtures: (1) 5 MPa, 423.15K (33 kg/m 3 ), and (2) 50 MPa, 323.15K (389 kg/m 3 ).For a system of 120 molecules, finite-size corrections for self-diffusivities in system (2) can reach up to 18% for CO 2 and around 7% for H 2 , while in system (1), H 2 and CO 2 diffusivities remain unaffected by system size.The finite-size corrected selfdiffusivity (D self ) of pure CO 2 at 2 MPa and 323.15K (2.2 × 10 −8 m 2 /s) agreed within a percent of the value reported by Moultos et al. 126 Given the considerable variation in densities of H 2 -CO 2 mixtures spanning nearly 2 orders of magnitude, finite-size effects are expected to affect the self-diffusivities of CO 2 and H 2 only at specific densities.Consequently, the Yeh− Hummer correction is applied when the correction term in eq 1 exceeds 1% of the uncorrected self-diffusivities.The reader is referred to Figures S17 and S18 for the plots of the finite-size effects.
The range between 10 −8 and 10 −6 m 2 /s for self-diffusivities in Figure 2, implies dense-gas and gas-like behavior of H 2 -CO 2 mixtures. 11As clearly shown in Figure 2, the self-diffusivities of H 2 are always larger than those of CO 2 at any pressure, temperature, or mixture composition by a factor between ca.2.5 and 6.To explain why the self-diffusivities of H 2 exceed those of CO 2 , the Stokes−Einstein (SE) relation 11,127 is used.The SE equation links the self-diffusivity of a microscopic entity in liquid, viscosity, temperature, and effective radius (or a hydrodynamic radius, R eff ) as 11,127 Here, stick boundary conditions are assumed, 11,127 resulting in a factor of 6π in the SE equation, in contrast to the 4π factor used with slip boundary conditions. 11,127It is also implicitly assumed that R eff in eq 7 is unaffected by T. 11,45 Treating CO 2 and H 2 as spherical entities with an R eff equal to the bond length between carbon and oxygen (1.16 Å), and half the bond length between the two hydrogen atoms in H 2 (0.37 Å), the ratio of their self-diffusivities as per eq 7 amounts to ca. 3.14, a value larger than unity.According to the SE relation, by the virtue of geometry, H 2 diffuses faster than CO 2 by roughly 3 times.In contrast to the expected ratio of ca.3.14, according to the SE relation, the ratio of self-diffusivities for H 2 and CO 2 ranges from 2.5 to 6, indicating significant deviations.These deviations, as shown in Figure 3a, are particularly pronounced at 5 MPa, possibly due to the gas-like behavior of the mixtures, suggesting potential limitations in the applicability of the SE relation for gaseous systems.The agreement with the SE relation improves at larger pressures.For instance, Figure 3b indicates deviations less than 20% at 50 MPa, where the behavior of the mixtures resembles a liquid.The reader is referred to Figure S20a−d of the Supporting Information for similar plots between 10 and 40 MPa.While the SE equation provides a qualitative understanding of the self-diffusivities of CO 2 and H 2 , quantitative predictions may deviate due to (1) the incorrect assumption that R eff remains constant with temperature, as also noted by Saric and colleagues, 44 and/or, (2) the oversimplified treatment of linear molecules like CO 2 and H 2 as spherical entities.
The impact of mole fraction of H 2 and temperature on the self-diffusivities is now discussed.As anticipated and shown in Figure 2a−f, self-diffusivities of CO 2 and H 2 increase as the proportion of H 2 increases in the mixture.This is because molecules can diffuse more in a gaseous medium than in liquid.To evaluate the effect of mole fraction of H 2 on selfdiffusivities, the values of D self for CO 2 and H 2 in H 2 -CO 2 mixtures are compared with corresponding values in pure fluids, as shown for the diffusivities at 5 MPa in Figure 3c.Adding H 2 to the mixture can lead to a nearly exponential rise in the self-diffusivities of both CO 2 and H 2 , with a maximum increase of approximately 2 to 4 times.A similar exponential increase in the self-diffusivities is observed also at other pressures, as shown in Figures S21 and S22 of the Supporting Information.
As anticipated, increasing temperature leads to larger selfdiffusivities of H 2 and CO 2 , as seen in Figure 2a−f due to the more gas-like nature of systems at larger T. The ratio of selfdiffusivities of H 2 and CO 2 relative to their values at 323.15 K is plotted in Figure 3d at 20 MPa.The self-diffusivity of CO 2 and H 2 increases up to 300% at x Hd 2 = 0.1, while a modest increase up to 60% is observed for both species at x Hd 2 = 0.9.A  108 The dashed line corresponds to the ideal-mixture limit, where Γ equals 1. Lines are colored as per the legend in (a).The raw data of the thermodynamic factors of CO 2 and H 2 for all mixtures are provided in Table S3 of the Supporting Information.
Industrial & Engineering Chemistry Research similar trend is observed at all pressures, barring 5 MPa, where the self-diffusivity of both species increases by a maximum of ca.100% at x Hd 2 = 0.1, but at x Hd 2 = 0.9, the increase is between 50 and 60%.At 5 MPa, however, self-diffusivities of H 2 and CO 2 increase by ca.60% at all mixture compositions, exhibiting a linear-like response to temperature.It can be concluded that self-diffusivities of H 2 and CO 2 in mixtures rich in H 2 > x ( 0.7) respond similarly to temperature at all pressures.For mixtures rich in CO 2 the response of selfdiffusivities to temperature varies depending on the pressure, indicating a complex behavior.This complexity may be associated with the proximity of CO 2 to the supercritical phase (T C = 304.15K, p C = 7.38 MPa). 125he self-diffusivities of CO 2 and H 2 increase with the mole fraction of H 2 and temperature but decrease with pressure.Conversely, mixture densities decrease with the mole fraction of H 2 and temperature while increasing with pressure.The contrasting responses of self-diffusivities and densities imply a possible correlation between the two quantities.Figure 3e,f show the variation of self-diffusivities for CO 2 and H 2 with density at x Hd 2 = 0.2 and 0.8, respectively.A linear-like relationship, displayed on a log−log scale, between the self-diffusivities of both H 2 and CO 2 and density is apparent.A linear fit is calculated to best describe this relationship at all mixture compositions.The average value of the slopes of the linear fits for the case of CO 2 is −1.00 ± 0.05 and −1.1 ± 0.1 for H 2 .Minor deviations from the calculated slope of around −1 are found for ρ exceeding ca.400 kg/m 3 , in Figure 3e.Similar deviations in the slopes are also observed for pure CO 2 for ρ ≳ 500 kg/m 3 .Barring these exceptions, it can be concluded that the self-diffusivities for CO 2 and H 2 are inversely proportional to the mixture density.When comparing Figure 3e,f, a significant vertical shift of the data points is observed, along with a minor horizontal shift.The shifts imply dependencies of self-diffusivities not only on the density but also on x Hd 2 and T. It appears improbable to establish a master equation linking self-diffusivities solely to densities.Studies by Van Loef 128 and Harris and Trappeniers 129 studying pure liquids, found that self-diffusivities are functions of density and temperature.In a binary mixture, an extra thermodynamic degree of freedom, the mole fraction of a component, must be considered alongside density and temperature to fully describe self-diffusivities.Establishing a function correlating selfdiffusivities to density, temperature, and x Hd 2 will require  S3 of the Supporting Information.additional efforts.Currently, no experimental data are available to compare the computed self-diffusivities in H 2 -CO 2 mixtures.Additional plots of self-diffusivities as a function of ρ are shown in Figure S24a−d.The computed self-diffusivities in H 2 -CO 2 mixtures are tabulated in Table S3 of the Supporting Information.

Thermodynamic Factors for Diffusion for CO 2 -H 2 Mixtures.
Figure 4 shows the thermodynamic factors calculated from the Gibbs excess energies obtained from REFPROP. 108The values of Γ for the thermodynamic states considered in this work are positive and vary between 0.25 and 1. Positive values of Γ signify that the mixtures are thermodynamically stable and do not undergo phase separation. 10A notable deviation from unity in the values of Γ for the CO 2 -H 2 mixtures suggests highly nonideal behavior.Γ of a binary system must approach unity when approaching the pure-component limit. 10,12t a fixed mixture composition, the values of Γ tend to decrease with pressure.For instance, in Figure 4, for an equimolar mixture of CO 2 -H 2 the value of Γ for the mixture drops from ca. 0.9 at p = 5 MPa (Figure 4a) to ca. 0.3 at p = 50 MPa (Figure 4f).This holds true for the values of Γ at other mixture compositions.For a fixed mixture composition, the value of Γ approaches 1 as the temperature rises from 323.15 to 423.15 K, regardless of the pressure.For instance, at p = 50 MPa (Figure 4f), the value of Γ increases from ca. 0.3 at T = 323.15K to ca. 0.7.At a specific temperature and pressure, the values of Γ exhibit a slight asymmetry ca.x Hd 2 = 0.5, see Figure 4c (p = 20 MPa).The thermodynamic factors of all CO 2 and H 2 mixtures considered in our study are tabulated and displayed in Table S3 of the Supporting Information.

MS and Fick-Diffusion
Coefficients for CO 2 -H 2 Mixtures.−119 To investigate finite-size effects in MS diffusivity calculations and to confirm the validity of eq 2 for our systems, equimolar H 2 -CO 2 mixtures at 50 MPa and 323.15K (389 kg/m 3 ) consisting of 120, 250, 500, 1000, and 2000 molecules were chosen.Both methods, extrapolation of the line fitted to the MS diffusivities up to N tot → ∞ and the prediction from eq 2, agreed within 3%.The finite-size correction to the MS diffusivity for the system of 120 molecules was ca.20%, underscoring the significance of finite-size effects in the computation of MS diffusivities.Consequently, the correction suggested in eq 2 is applied if it exceeds 1% of the computed MS-diffusivities.The reader is referred to Figure S27 for the plots of the finite-size effects in MS diffusivities.
Figure 5a−f shows the finite-size corrected MS diffusivities (D̵ MS ) for various CO 2 -H 2 mixtures, compared to the Darken prediction from eq 3. Finite-size corrected self-diffusivities of CO 2 and H 2 are used to compute the Darken predictions.The Darken prediction underestimates the computed MS diffusivities for the mixtures considered in our study.The maximum deviations between eq 3 and the computed MS diffusivities can reach ca.30% around x Hd 2 ≈ 0.7, while a minimum deviation of about 5% occurs at x Hd 2 ≈ 0.1.
The significant deviations between the computed MS diffusivities and Darken predictions imply a significant nonideality of the H 2 -CO 2 mixtures.Since the thermodynamic factor (Γ) reflects the nonideality of the mixture, it is valuable to explore the potential correlation between D̵ MS − D Darken and Γ, as suggested by Wolff et al. 122 Figure 6a   4. Similar to the self-and MS diffusivities, finite-size corrections are disregarded if the corrections are less than 1% of the computed values, as shown in eq 4. Figure 7a−f shows the Fick diffusivities as a function of the mole fraction of hydrogen.The statistical uncertainties in the computed Fick diffusivities are noticeable, particularly below ca. 10 MPa, as shown in Figure 7a−c.This increased uncertainty, such as at 5 MPa, 323.15 K, and x Hd 2 = 0.5, is attributed to sparse intermolecular collisions in the gas phase, leading to slow convergence toward diffusive behavior of the cross-correlation of mean-squared displacements of H 2 and CO 2 .To enhance accuracy, more independent simulations could be performed, but we refrained from doing so due to the significant computational costs associated with gas-phase MD simulations.The Fick diffusivities computed for all pressures, temperatures, and mixture compositions vary between 10 −8 m 2 /s at 50 MPa and 10 −6 m 2 /s at 5 MPa, indicating dense-gas/ liquid-like behavior at high pressures and gas-like behavior at low pressures.In contrast to the smooth increase observed for self-diffusivities of H 2 and CO 2 in Figure 2a−f, the behavior of Fick diffusivities varies as the content of H 2 increases in the mixture depending on the pressure.At 5 MPa in Figure 7a, Fick diffusivities exhibit minimal dependence on mixture composition, whereas at 50 MPa in Figure 7f, they can increase by 3 to 4 times with higher H 2 enrichment in the mixture.The variation in Fick diffusivities contrasts with the consistent increase observed in the self-diffusivities of H 2 and CO 2 with the content of H 2 in the mixture.
In general, temperature elevation increases Fick diffusivities regardless of mixture composition or pressure.The extent of this increase, however, varies for different mixtures.Notably, for mixtures abundant in H 2 (x Hd 2 = 0.9), Fick diffusivities consistently double when the temperature rises from 323.15 to 423.15 K across all pressures.Conversely, for mixtures rich in CO 2 (x Hd 2 ≲ 0.2), the increase ranges between ca.100% (at 5 MPa) and ca.300% (at 20 MPa) depending on pressure, emphasizing the significant impact of temperature on Fick diffusivities, particularly for CO 2 -rich mixtures.The effect of temperature on Fick diffusivities aligns with trends seen in the self-diffusivities of both CO 2 and H 2 .

Prediction of Fick Diffusion Coefficients Using Kinetic Theory, Fuller Correlation, the Corresponding States
Principle, and Moggridge's Correlation.Due to the lack of experimental data or a specific theory to predict Fick diffusivities in CO 2 -H 2 mixtures at pressures larger than atmospheric pressure, the computed diffusivities are compared  S3 of the Supporting Information.
−134 Figure 8a shows the relative deviations of the computed Fick diffusivities from the kinetic theory expression of Chapman and Enskog.It can be expected that kinetic theory effectively predicts Fick diffusivities when mixtures behave ideally, typically observed at low pressures, high temperatures, and in mixtures rich in H 2 .In line with this expectation, a close agreement within 5% is evident at 423.15 K and 5 MPa for all mixtures.Below 423.15 K, predictions of kinetic theory are accurate within ca.10% at 373.15 K and ca.20% at 323.15 K.The computed Fick diffusivities of mixtures rich in H 2 (x Hd 2 ≳ 0.8) can be predicted within 15% regardless of the pressure.As anticipated, a severe loss of accuracy in kinetic theory (up to 75%) is observed for CO 2 -rich mixtures (x Hd 2 ≳ 0.5) at 25 and 50 MPa, regardless of the temperature.In conclusion, predictions of D Fick from kinetic theory are reliable up to ca. 20% at pressures of 5 MPa for CO 2 -H 2 mixtures.
−32 The relative deviations between the computed and predicted diffusivities of CO 2 -H 2 mixtures are shown in Figure 8b. Figure 8a,b exhibits striking similarity, indicating predictive capabilities of the FSG correlation comparable to the kinetic theory.This is not surprising, provided that Fuller's empirical relation is based on the Chapman and Enskog expressions from kinetic theory.It can be concluded that the kinetic theory and Fuller correlation can reliably estimate the Fick diffusivities within 25% in H 2 -CO 2 mixtures up to 5 MPa and temperatures between 323.15 and 423.15 K.For pressures below ca. 5 MPa and temperatures above ca.423.15 K, the accuracy of predicted Fick diffusivities from kinetic theory, and the Fuller correlation are expected to increase, consistent with findings in the literature. 11,13o test if accurate predictions of Fick diffusivities can be made at larger pressures, the performance of the corresponding state principle is examined.The ratio pD Fick /p 0 D 0 Fick is evaluated at all thermodynamic state points, where D 0 Fick is the Fick diffusivity evaluated at a reference pressure p 0 = 0.1 MPa and the corresponding temperature T. When the ratio equals unity, D Fick at any pressure and temperature can be predicted since the normalization factor is readily accessible.Figure 8c shows that the ratio is within 3% of unity at 5 MPa and 423.15 K. Like kinetic theory and FSG correlation, predictions are less accurate (deviations up to ca. 40%) at high pressures for mixtures rich in CO 2 .Conversely, mixtures rich in H 2 (x Hd 2 ≳ 0.8) agree within 15−20% with kinetic theory and Fuller correlation.Rather than using pressure p to assess the ratio pD Fick /p 0 D 0 Fick , an alternative approach is to use density ρ, as illustrated in Figure S31 of the Supporting Information.While the performance of the corresponding states principle remains qualitatively similar, there is a marginal improvement of ca.10% in the prediction of Fick diffusivities at 5 MPa, as shown in Figure S31 of the Supporting Information.In a nutshell, predicting Fick diffusivities is feasible through the kinetic theory, FSG correlation, or the corresponding states principle, specifically for mixtures at low pressures and high temperatures or enriched with hydrogen.
Due to the superior performance of Moggridge's correlation, specifically at high pressures and low temperatures, this empirical method could be advantageous for predicting Fick diffusivities at these conditions (see Figure 6a,b).By combining eqs 4 and 8, Moggridge's correlation can be extended to describe Fick diffusivities as As anticipated, in Figure 8d, the Fick diffusivities predicted using eq 9 closely match the computed values, staying within 5% at 25 and 50 MPa for temperatures of 323.15 and 348.15 K.For CO 2 -rich mixtures (x Hd 2 ≲ 0.2), the computed and predicted Fick diffusivities deviate by less than ca.5%, irrespective of pressure and temperature.In sharp contrast, for H 2 -rich mixtures (x Hd 2 ≳ 0.7) and mixtures at low pressures, deviations from the computed Fick diffusivities can reach approximately 35%.
In summary, Chapman and Enskog's kinetic theory, FSG correlation, and the corresponding states principle perform well for ideal-gas-like mixtures.Fick diffusivities for H 2 -CO 2 mixtures below 5 MPa can be accurately predicted within 5− 10%.These kinetic theory-based approaches also demonstrate good performance for mixtures rich in H 2 (x Hd 2 ≳ 0.8), with an accuracy of approximately 15−20%.Deviations can be as high as 75% when predicting Fick diffusivities of nonideal mixtures.For higher pressures and CO 2 -rich mixtures (x Hd 2 < 0.5), using Moggridge's correlation for MS diffusivities yields Fick diffusivities accurate to within 5−10%.Values of MS and Fick diffusivities for all CO 2 and H 2 mixtures considered in our study are provided in Table S3 of the Supporting Information.

Comparison of MD-Computed Fick Diffusivities with
Experiments.We aim to (indirectly) compare the Fick diffusivities of CO 2 -H 2 mixtures from MD simulations to existing experiments at p atm = 1 atm (see Table 1).Computing Fick diffusivities at atmospheric pressure using MD simulations is possible but impractical due to the severe computational demands accompanying gas phase simulations.
The principle of corresponding states can indirectly compute Fick diffusivity at an unknown pressure, in our case, the atmospheric pressure, by extrapolating known data from a reference pressure. 11,13For comparison with experiments, Fick diffusivities computed from MD simulations at a reference pressure of 5 MPa are used to estimate mutual diffusivities at atmospheric pressure.Comparison of the computed Fick diffusivities at 5 MPa with those estimated using kinetic theory at 0.1 MPa reveals that the product of pressure (or density) and mutual diffusivities remains nearly constant (see Figure 8e), with deviations of approximately 20%.This trend is consistent when analyzing densities instead of pressures, showing agreement within 10%, as shown in Figure S31 of the Supporting Information.The principle of corresponding states for the product of density and Fick diffusivity offers a practical approach to extrapolating computed Fick diffusivities to conditions of experimental relevance.The Chapman-Enskog expression for mutual diffusivity, denoted as D AB , of rigid spheres with unequal masses is given by 11,17,18  AB mix (11)   The functional form of eq 11 suggests that the product of density and Fick diffusivity can be equated between two thermodynamic states under the same temperature and mixture composition, as M mix depends on the mole fraction of species A. Although eq 10 assumes spherical molecules, the equation has been shown to be effective at low pressures, irrespective of the shape of the molecules. 11,13,17rom Figure 8e, we observe good agreement (within ca.15%) between the Fick diffusivities predicted from the principle of corresponding states using eq 11 (filled diamonds) for an equimolar mixture composition, and experimental data by Gavril et al., 135 Saxena and Mason, 136 and Weissman. 137he asymmetric error bars associated with each symbol in Figure 8e, representing the values of D Fick calculated using eq 11 with mixture compositions of x Hd 2 = 0.1 and 0.9, respectively, imply a weak dependence of Fick diffusivities on mixture concentration.Accounting for this variability in the Fick diffusivities by choosing different mixture compositions enhances agreement with the experiments (between ca. 1 and 5%).Leveraging the close agreement between MD computed densities and REFPROP 108 predictions (see Figures 1 and S5−S7 of the Supporting Information), Fick diffusivities can be evaluated at experimentally relevant conditions for mixture compositions of 0.001 and 0.999, yielding consistent results.Experimental measurements used in Figure 8e are not corrected for composition effects, known to affect diffusivities by ca.5%, 19 potentially influencing agreement between MD simulations and experimental measurements.Using the product of pressure and Fick diffusivity in eq 11 yields similar conclusions.In summary, MD simulations provide reliable estimates of Fick diffusivities at conditions relevant to atmospheric pressure and temperatures between 323.15 and 423.15 K. Solubilities of CO 2 in NaCl brine computed from simulations are compared to experimental data obtained from Weibe and Gaddy (WG), 137 Takenouchi and Kennedy (TK), 33 Duan and Sun (DS), 56 and Todheide and Frank (TF). 138In Figure 9a−c the computed solubilities at 323.15 K deviate from experimental measurements by between 20% at p ≲ 10 MPa and 30% at p ≳ 30 MPa.In contrast, an agreement between 5 and 10% is observed between simulations and most experiments at 423.15 K and for all salt concentrations.Previous simulation studies by Liu et al. 60,61 and Orozco et al. 62 report similar disagreements between experiments and MC simulations at low temperatures.Todheide and Frank (TF), 138 Takenouchi and Kennedy (TK), 170,171 and Duan and Sun (DS) 56,88 are considered for the comparison.(d,e) Reduction in the solubility of CO 2 as a function of NaCl concentration in brine.The logarithm of the ratio of CO 2 solubility in brine to its solubility in pure water is plotted as a function of NaCl concentration.Linear fits to the data are shown as dashed lines.The raw data of the solubilities are provided in Table S4 of the Supporting Information.
The salting-out effect occurs in solutions when certain salts (ionic compounds) are added, leading to a decrease in the solubility of nonionic solutes in the solvent, typically organic compounds.This is due to the presence of ions in the solution altering the interactions between solvent and solute molecules, resulting in reduced solubility. 139The salting-out effect has been well studied in the context of CO 2 solubility in aqueous solutions of many ionic salts, including NaCl. 55,56,88The influence of salt concentration on the solubility of CO 2 becomes evident when plotting the ratio of x COd 2 relative to its solubility in pure water under the same thermodynamic conditions, as a function of the salt concentration.This is the well-known Setschenow equation 140 = i k j j j j j j y where x c CO 0 2 is the reference solubility of CO 2 in H 2 O at the same pressure at a chosen temperature.In Figure 9d,e, the solubility ratios are plotted on the vertical axis as a function of the salt concentration, in a format consistent with the Setschenow equation (eq 12).A linear relationship is observed between the solubility ratios and salt concentration, consistent with eq 12, exhibiting a slope of ca.−0.21 regardless of temperature and pressure.The negative slope implies a reduced solubility of CO 2 in brine upon increasing the salt concentration.Specifically, the solubility of CO 2 reduces, at any chosen pressure and temperature, by ca.20% when a mole of NaCl is added to pure water results.A subsequent rise in salt concentration to 2 mol NaCl/kg H 2 O results in a further reduction of CO 2 by another ca.20%.
As an additional check of our simulations, the fugacity coefficients of CO 2 in the gas phase computed along the VLE of CO 2 -NaCl brine at a given composition (y COd 2 , y Hd 2 O ), pressure, and temperature are compared to the data from REFPROP. 108Fugacity coefficients are a direct reflection of the excess chemical and mixture density, as indicated by eq 6.In general, the values of φ COd 2 and H O 2 were below unity, signifying the nonideality of CO 2 -H 2 O mixtures in the gas phase.The fugacity coefficients were found to decrease with  79 Kling and Maurer (KM), 172 Wiebe and Gaddy (WG), 173 12e,f.Clearly, from Figure 12e,f, a substantial quantity of water is present in the gas phase at 423.15 K in comparison to 323.15 K.At 323.15 K, the water content in the gas phase is unaffected by pressure, while the water content decreases by ca. 3 times between 5 and 50 MPa at 423.15 K.The quantity of water in the gas phase is insensitive to the NaCl concentration in brine.
Additional plots pertaining to the effect of temperature, pressure, and salt concentration on the VLE CO 2 -NaCl brine systems is shown in Figures S35−S37 of the Supporting Information.
3.7.Phase Equilibria of H 2 -NaCl Brine Mixtures.Figure 10a shows the solubilities of H 2 in pure water at 323.15 and 423.15K for pressures between 5 and 50 MPa.In Figure 10a, the H 2 solubility in pure H 2 O shows a linear increase with pressure, indicating adherence to the Henry regime across the pressure range, with solubilities approaching 0 in the limit of p → 0. Increasing the pressure from 5 to 50 MPa leads to a near-10-fold increase in the solubility of H 2 .Interestingly, the impact of pressure on H 2 solubility remains consistent regardless of salt concentration or temperature.In contrast to CO 2 , increasing temperature consistently enhances the dissolution of H 2 in both H 2 O and NaCl brine.The impact of temperature however is less pronounced in NaCl brine, where a smaller increase in solubility is observed compared to pure H 2 O.For instance, an increase in temperature from 323. 15  The solubilities of H 2 in NaCl brine computed from MC simulations are compared to the available literature data in Figure 10a.Computed solubilities agree within ca.5−10%, of the experimental data from Tori ́n-Ollarves and Trusler 79 for the entire range of pressures at 323.15 K. Disagreements between simulations and experiments arise at 423.15 K, with the solubilities being overestimated by simulations by a maximum of ca.27% at p = 30 MPa.In Figure 10a, the disagreement between simulations and experiments at 423.15 K increases with pressure.The water content in the gas phase (  79 with deviations between 5 and 10% for various pressures. The salting-out effect in H 2 solubility on NaCl brine has been explored in experiments by Chabab et al. 57,59 Plotting the solubility ratios as a function of the salt concentration according to the Setschenow equation (eq 12) yields a linear relation, as shown in Figure 10c,d.Unlike CO 2 in Figure 9d,e, the fitted lines have different slopes, suggesting a minor sensitivity of the salting-out effect on the pressure.The slopes of the fitted lines in Figure 10c,d  in NaCl brine is examined as a function of fugacity rather than pressure, as fugacity serves as a measure of the chemical potential responsible for gas dissolution in water.The fugacity of component i (f i ), where i is either CO 2 or H 2 , is 52,53 = f p T y yp ( , , ) where ϕ i is the fugacity coefficient, as defined in eq 6, y i is the mole fraction of species i in the gas phase, and p is the pressure.
The fugacity coefficients and the fugacities of CO 2 and H 2 are shown in Figure 11a−d as a function of pressure between 323.15 and 423.15K along the VLE of CO 2 -H 2 -NaCl brine.In Figure 11a, the fugacity coefficients of CO 2 behave similarly to those in the CO 2 -NaCl brine system.Specifically, ϕ COd 2 decreases with pressure but increases with temperature.Computed fugacity coefficients derived from gas phase compositions along the VLE of CO 2 -H 2 -NaCl brine agree within ca.5−10% of the REFPROP predictions.In Figure 11c, fugacities of CO 2 computed using eq 13 vary linearly with pressure, as indicated by the fits at 323.15 and 423.15 K. Unlike CO 2 , the fugacity coefficients of H 2 shown in Figure 11b exhibit an increase with pressure but a decrease with  S6 of the Supporting Information.

Industrial & Engineering Chemistry Research
temperature, qualitatively reproducing its behavior in the binary mixture of H 2 -NaCl brine.Fugacities of H 2 computed using eq 13 vary linearly with pressure, as indicated by the fits at 323.15 and 423.15K in Figure 11d.To facilitate comparison, open symbols representing the fugacities of CO 2 and H 2 in the binary systems of CO 2 -NaCl and H 2 -NaCl are plotted alongside the ternary system in Figure 11c,d.The fugacities of CO 2 and H 2 in the ternary CO 2 -H 2 -NaCl systems are lower than those in their respective binary systems.For instance, in Figure 11c, the fugacity of CO 2 in the ternary system at 50 MPa and 423.15K is lower by ca.45% than f COd 2 in the binary system.Similarly, in the case of H 2 in Figure 11d, the fugacities in the ternary system are lower by ca.40% at 50 MPa and 423.15 K.
The solubilities of CO 2 and H 2 in NaCl brine solutions are shown in Figure 12a−d, comparing binary and ternary systems.The solubility of CO 2 varies quadratically with the fugacity of CO 2 in both binary and ternary systems, as shown in Figure 12a.In the low-fugacity regime ( f COd 2 ≲ 4 MPa), the solubilities of CO 2 in NaCl brine are similar for the binary and ternary systems.For intermediate and large fugacities (f COd 2 ≳ 4 MPa), the solubility of CO 2 in NaCl brine becomes substantially different in the binary and ternary systems.Specifically, CO 2 solubility in the ternary system is lower compared to the binary system at (approximately) the same fugacity.An even stronger suppression in the dissolution of CO 2 occurs at 423.15 K, as shown in Figure 12b.The solubilities of H 2 in NaCl brine shown in Figure 12c,d is analyzed next.In the low-fugacity regime (f Hd 2 ≲ 40 MPa), the solubilities of H 2 are nearly identical in the binary and ternary systems.Differences in the solubilities arise for f Hd 2 ≳ 40, wherein a greater H 2 dissolution is evident in the ternary case when compared to the H 2 -NaCl brine system.The enhanced dissolution occurs also at 423.15 K, as shown in Figure 12d.In contrast to CO 2 , the solubility of H 2 as a function of f Hd 2 is enhanced in ternary systems.
Figure 12e,f show the water content in the gas phase for both binary and ternary systems.The water content in the gas phase varies negligibly with the salt concentration in brine.In Figure 12e at 323.15 K, the mole fraction of H 2 O in the gas phase decreases with pressure for CO 2 -H 2 -H 2 O systems.As expected, the mole fraction of water in the ternary system falls between those of the CO 2 -H 2 O and H 2 -H 2 O binary systems.At 423.15 K, the water content in the gas phase is higher, yet it follows a trend similar to that at 323.15 K.In conclusion, an enhanced dissolution of water is observed compared to the H 2brine system but suppression compared to the CO 2 -brine system.

CONCLUSIONS
Understanding the thermodynamic behavior of mixtures containing H 2 and CO 2 is essential for addressing hydrogen losses in underground hydrogen storage within depleted gas fields.The loss of purity of H 2 resulting from molecular diffusion follows Fick's law, highlighting the importance of knowing the Fick diffusion coefficients for a variety of temperatures, pressures, and mixture compositions.The dissolution of H 2 in salt water, particularly in the presence of CO 2 , presents another avenue for hydrogen loss.Thus, a comprehensive knowledge of the Fick diffusion coefficients for H 2 -CO 2 mixtures and their phase equilibria with NaCl brine is crucial.The thermodynamic conditions relevant to subsurface hydrogen storage typically span temperatures of 323.15 to 423.15 K, pressures ranging from 5 to 50 MPa, and salt concentrations of 1 to 2 mol of NaCl/kg of H 2 O.
The Fick diffusivities of H 2 -CO 2 mixtures were investigated via MD simulations at storage-relevant conditions.At low pressures (5 Mpa), gas-like diffusivities of ca. 10 −6 m 2 /s remained composition-independent.The kinetic theory accurately predicted these diffusivities within 5−10% at 423.15 K, with deviations of up to 20% at 323.15 K.At high pressures (40−50 Mpa), kinetic theory predictions deviated up to 75% for CO 2 -rich mixtures and around 15% for hydrogenenriched ones.Approaches based on kinetic theory, such as the Fuller correlation 30−32 and the corresponding-states approaches, 11,13 yielded similar results.Kinetic theory-based methods can predict Fick diffusivities within 15% for H 2 -CO 2 mixtures exhibiting near-ideal gas behavior (1) at low pressures (ca. 5 Mpa) and high temperatures (423.15K) and ( 2) for H 2enriched mixtures (mole fraction of H 2 ≳ 0.8).Moggridge correlation 122,130−134 provides an alternative to accurately predicting diffusion coefficients within 5−10% at high pressures (40−50 MPa) and low temperatures (T ≲ 348.15 K) for CO 2 -rich mixtures, making it suitable for the highpressure regime.For intermediate pressures (10−35 MPa), molecular simulations are preferred.The efficacy of the Moggridge correlation stems from the fluid-like behavior of mixtures at high pressures.Mixtures at 323.15 K and pressures of 15−25 MPa exhibit oscillatory trends in Fick diffusivities due to the proximity of CO 2 to its critical point.
Phase equilibria of CO 2 -H 2 -NaCl mixtures are examined using CFCMC simulations showing close agreement with experimental results for CO 2 solubilities in NaCl brines.The Setschenow equation 140 used to describe the salting-out effect of CO 2 captures the exponentially decreasing solubilities of CO 2 in the liquid-rich phase with an increasing concentration of NaCl in brine.The water content in the gas phase decreased with pressure at 423.15 K but remained nearly constant at 323.15 K, regardless of salt concentration.H 2 solubilities in brine followed the Henry regime, varying linearly with pressure.The impact of salt concentration on H 2 solubility relative to pure water is quantified using the Setschenow equation.Similar to the CO 2 -NaCl brine system, water concentrations in the gas phase consistently decrease with pressure at both 323.15 and 423.15 K. Ternary systems of H 2 -CO 2 -NaCl brine mixtures revealed a suppressed CO 2 solubility and an enhanced H 2 solubility, compared to the solubilities in the respective binary mixtures.The water content in the gas phase is enhanced compared to H 2 -NaCl brine systems but suppressed compared to CO 2 -NaCl brine systems.
Mixture densities, mixture compressibility factors, mixture viscosities, self-diffusion coefficients, finite system-size effects on self-diffusion coefficients, selfdiffusivities, ratio of self-diffusion coefficients, ratio of self-diffusivities, effect of temperature on self-diffusion coefficients, density dependence of self-diffusion coefficients, Maxwell-Stefan diffusion coefficients from crosscorrelations of displacements, finite-size effects for the Industrial & Engineering Chemistry Research mutual diffusion coefficients, Fick diffusion coefficients, effect of mixture composition and temperature on Fick diffusion coefficients, principle of corresponding states for Fick diffusivities, effect of pressure, temperature, and NaCl concentration on the solubility, force field parameters, densities, compressibilities, and total energies from MD simulations, viscosities, thermodynamic factors, and self, Maxwell-Stefan, and Fick diffusion coefficients, phase equilibria, and fugacity coefficients (PDF) ■

Figure 1 .
Figure 1.H 2 -CO 2 mixture densities (ρ) and viscosities (η) at various pressures and temperatures.Computed properties (symbols) are compared to the REFPROP 108 data (solid lines).Uncertainties in simulation results are smaller than the symbol sizes.Lines are colored as per the legend in subfigure (a).Raw data on the densities and viscosities for all mixtures are provided in TablesS2 and S3of the Supporting Information, respectively.

3 . 3 .
Self-Diffusion Coefficients of CO 2 -H 2 Mixtures.Self-diffusivities computed from MD simulations might need to

Figure 2 .
Figure 2. Self-diffusion coefficients (D self ) of CO 2 (triangles) and H 2 (squares) in mixtures of CO 2 and H 2 , computed from MD simulations.The computed self-diffusion coefficients are plotted as a function of the mole fraction of H 2 in the mixture at various pressures.The uncertainties in the computed self-diffusion coefficients are smaller than the symbol sizes.The dashed lines act as guides to the eye.Lines are colored as per the legend in (a).Self-diffusivities of CO 2 and H 2 plotted on a linear scale are displayed in Figure S19a−f of the Supporting Information to aid comparison.The raw data of the self-diffusion coefficients of CO 2 and H 2 for all mixtures are provided in TableS3of the Supporting Information.

Figure 3 .
Figure 3. (a,b) Ratio of self-diffusivities of H 2 and CO 2 as a function of H 2 mole fraction at 5 and 50 MPa, respectively.The dashed line represents the ratio predicted by the Stokes−Einstein relation.(c) Comparison of the self-diffusivities of CO 2 (triangles) and H 2 (squares) in CO 2 -H 2 mixtures with their values in pure fluids at the same temperature at 5 MPa.(d) Comparison of the self-diffusivities of CO 2 (triangles) and H 2 (squares) in CO 2 -H 2 mixtures with their values at T 0 = 323.15K at 20 MPa.(e,f) Self-diffusivities of CO 2 (triangles) and H 2 (squares) as a function of mixture densities at x Hd 2 = 0.2 and 0.8, respectively, with symbols colored according to the legend in (c).The ratios of self-diffusivities in (c) are displayed on a linearly scaled vertical axis in Figure S22 of the Supporting Information along with other pressures.Similarly, Figure S25 displays the densities and self-diffusivities in (e,f) plotted on linearly scaled horizontal and vertical axis.Figures S20−S25 of the Supporting Information contain additional plots for (a,b,d) at various pressures.

Figure 4 .
Figure 4. Thermodynamic factors (Γ) of CO 2 -H 2 mixtures plotted as a function of the mole fraction of H 2 (x Hd 2 ) for various pressures.Values of Γ are obtained from the Gibbs excess energies of H 2 -CO 2 mixtures from REFPROP.108The dashed line corresponds to the ideal-mixture limit, where Γ equals 1. Lines are colored as per the legend in (a).The raw data of the thermodynamic factors of CO 2 and H 2 for all mixtures are provided in TableS3of the Supporting Information.

Figure 5 .
Figure 5. Finite-size-corrected Maxwell−Stefan diffusion coefficients as a function of the mole fraction of hydrogen (x Hd 2 ).The dashed-dotted lines are predictions from the Darken relation (eq 3, where the self-diffusion coefficients for CO 2 and H 2 are taken from Figure 2).Lines are colored as per the legend in (d).The raw data of the Maxwell−Stefan diffusion coefficients of CO 2 and H 2 mixtures are provided in TableS3of the Supporting Information.
shows the relative deviation of the MS diffusion coefficient from the Darken relation as a function of the thermodynamic factor, for all mixtures considered in our study.The Moggridge correlation122,130−134  when recast into a suitable form relative deviation of D̵ MS from D Darken and the thermodynamic factor, with the best agreement achieved at low temperatures.Figure6bhelps carefully assess the performance of the Moggridge correlation122,130−134  by showing the percentage deviations between the symbols and the dashed line in Figure6a.The deviations in Figure6bare within 20% when describing mixtures below 348.15K and above 40 MPa.In other words, the deviation between the computed MS diffusivities and the Darken relation is accurately captured in the liquid-like phase.The Moggridge correlation becomes less accurate when considering H 2 -CO 2 mixtures at high temperatures and low pressures.For instance, the black cross representing 5 MPa and 423.15K exhibits the maximum deviation.To conclude, Moggridge correlation can be advantageous in predicting the deviation between the computed MS diffusivities and the Darken relation for mixtures behaving as liquids.In summary, knowing the self-diffusivities of H 2 and CO 2 and the thermodynamic factor can facilitate the prediction of D̵ MS of the mixture with an accuracy of approximately 20% for mixtures between 40 and 348.15 MPa.Fick diffusivities for CO 2 -H 2 mixtures are the products of the MS diffusivities and the thermodynamic factors, as shown in eq

Figure 6 .
Figure 6.(a) Relative differences between the Darken relation and the computed MS diffusivities plotted as a function of the thermodynamic factor Γ for all pressure, temperature, and mixture compositions.The dashed line is the Moggridge correlation, 122,130−134 see eq 8.The symbols code for the pressure and the colors code for the temperature.For example, a blue star would correspond to 20 MPa and 348.15 K. (b) Relative deviations, in percentage, of the symbols in (a) from the Moggridge correlation (σ Mog ).The dashed line corresponds to a deviation of 20% to aid visual comparison.The percentage deviation extends across 5 orders of magnitude.

Figure 7 .
Figure 7. Finite-size-corrected Fick diffusion coefficients as a function of the mole fraction of hydrogen (x Hd 2 ).Lines are colored as per the legend in (d).The dashed lines act as guides to the eye.The raw data of the Fick diffusion coefficients of CO 2 and H 2 for all mixtures are provided in TableS3of the Supporting Information.
where D AB [m 2 s −1 ] is the mutual diffusivity of the A-B mixture, M AB is the reduced mass of A and B, k B is the Boltzmann constant, n [m −3 ] is the number density, σ AB [m] is the characteristic length of the intermolecular force law, and Ω(T) is the temperature-dependent collision integral for diffusion.Multiplying both sides of eq 10 by nM mix /N Av , whereM mix = M A x A + m B x Bis the molecular mass of the mixture in kg•mol −1 , and N Av [mol −1 ] is the Avogadro number, which allows us to express it as

3 . 6 .
Phase Equilibria of CO 2 -NaCl Brine Mixtures.The solubilities of CO 2 in pure water and NaCl brine at 323.15 and 423.15K between 5 and 50 MPa and for salt concentrations of 1 and 2 mol of NaCl/kg of H 2 O.At 323.15 K in Figure 9a−c, CO 2 undergoes a significant linear-like increase in solubility within the liquid-rich phase up to p ≲ 10 MPa.Beyond this pressure, the solubility continues to increase linearly, albeit more gradually.The distinct slopes in these two regions reflect a transition from a gas-like to a liquid-like behavior of CO 2 .At 423.15 K, x COd 2 varies almost quadratically with pressure.As p → 0, extrapolation of x COd 2 shows that the solubility of CO 2 also approaches 0, in agreement with the definition of the Henry regime.In general, insensitive to the concentration of NaCl in H 2 O, increasing pressure from 5 to 50 MPa results in a 2-fold increase in x COd 2 at 323.15 K and a 5-fold increase in x COd 2 at 423.15 K.The effect of temperature on the solubility of CO 2 is nontrivial.Beyond a certain threshold pressure p Th , the solubility of CO 2 in the NaCl brine increases with temperature.Below p Th , temperature has the opposite effect, leading to a decrease in CO 2 solubility.This threshold pressure appears to depend on the NaCl concentration in H 2 O.

Figure 9 .
Figure 9. (a−c) CO 2 solubilities in NaCl brine computed from MC simulations (symbols) compared to experiments (lines) for varying brine concentrations.Experimental studies of Weibe and Gaddy (WG),169 Todheide and Frank (TF),138 Takenouchi and Kennedy (TK),170,171 and Duan and Sun (DS)56,88 are considered for the comparison.(d,e) Reduction in the solubility of CO 2 as a function of NaCl concentration in brine.The logarithm of the ratio of CO 2 solubility in brine to its solubility in pure water is plotted as a function of NaCl concentration.Linear fits to the data are shown as dashed lines.The raw data of the solubilities are provided in TableS4of the Supporting Information.

Figure 10 .
Figure 10.(a) Comparison of simulations (symbols) with experimental (lines) predictions for H 2 solubilities in pure H 2 O. Experimental studies of Tori ́n Ollarves and Trusler (TT),79 Kling and Maurer (KM),172 Wiebe and Gaddy (WG),173 and Chabab et al. (SC) 57 are considered for the comparison.(b) Comparison of the water content in the gas phase between our simulations and the experiments from Tori ́n Ollarves and Trusler (TT).Predictions from the equation of state from Tori ́n Ollarves and Trusler (TT) are also plotted alongside.(c,d)Reduction in the solubility of H 2 as a function of NaCl concentration in brine.The logarithm of the ratio of H 2 solubility in brine to its solubility in pure water is plotted against NaCl concentration.Linear fits to the data are shown as dashed lines.The raw data of the solubilities are provided in TableS5of the Supporting Information.
Figure 10.(a) Comparison of simulations (symbols) with experimental (lines) predictions for H 2 solubilities in pure H 2 O. Experimental studies of Tori ́n Ollarves and Trusler (TT),79 Kling and Maurer (KM),172 Wiebe and Gaddy (WG),173 and Chabab et al. (SC) 57 are considered for the comparison.(b) Comparison of the water content in the gas phase between our simulations and the experiments from Tori ́n Ollarves and Trusler (TT).Predictions from the equation of state from Tori ́n Ollarves and Trusler (TT) are also plotted alongside.(c,d)Reduction in the solubility of H 2 as a function of NaCl concentration in brine.The logarithm of the ratio of H 2 solubility in brine to its solubility in pure water is plotted against NaCl concentration.Linear fits to the data are shown as dashed lines.The raw data of the solubilities are provided in TableS5of the Supporting Information.

2 )
shown in Figure 10b at 323.15 and 423.15K aligns with the experimental data from Tori ́n-Ollarves and Trusler,

Figure 12 .
Figure 12.Mole fractions of CO 2 (x COd 2 ) and H 2 (x Hd 2 ) in the liquid-rich phase are presented for the H 2 -CO 2 -NaCl system (closed symbols) and compared with the CO 2 -NaCl and H 2 -NaCl systems (open symbols).(a,b) x COd 2 as a function of the fugacity of CO 2 in the gas-rich phase at 323.15 and 423.15 K, respectively.Quadratic fits for x COd 2 in the H 2 -CO 2 -NaCl system are represented by dashed lines.(c,d) serve the same purpose for H 2 .(e,f) Comparison of the water content in the CO 2 -H 2 -H 2 O (ternary) system with the CO 2 -H 2 O and H 2 -H 2 O (binary) systems at 323.15 and 423.15 K.The raw data of the solubilities are provided in TableS6of the Supporting Information.

Table 1 .
List of Experimental Studies Reporting Diffusion Coefficients of Binary Mixtures of Hydrogen and Carbon Dioxide a O 2 in ternary mixtures of C 14 O 2 , CO 2 , and H 2 for temperatures of 308.15 and 373.15 K, and pressures of 0.6 and 25.33 MPa.Literature covering experimental studies of trace solutes in dense gases, such as supercritical CO 2 , has been summarized by Poling et al.

up to pressures of 40 MPa. Following Krichevsky and Kasarnovsky, 85 Tori ́n-Ollarves and Trusler 79 developed a correlation describing the solubility of H 2 in H 2 O and NaCl brines for temperatures between 273.15 and 423.15 K and for pressures up to 100 MPa. Although the solubilities of H 2 in salt-free water, as reported by Chabab et al., 57 conformed with
Table S5 of the Supporting Information.pressureand increase with temperature.For all thermodynamic conditions investigated in our study, the values of φ COd 2 and to 423.15 K results in a ca.80% increase in solubility in pure H 2 O, whereas the increase is around 40% in brine containing 1 mol NaCl/kg of H 2 O, and nearly 30% in brine containing 2 mol NaCl/kg of H 2 O.

. Phase Equilibria of CO 2 -H 2 -NaCl Brine Mixtures.
fall between [−0.11, −0.19].The solubility of H 2 reduces, at any chosen pressure and temperature, by ca.15% when a mole of NaCl is added to pure water results.A subsequent rise in salt concentration to 2 mol NaCl/kg H 2 O results in a further reduction of CO 2 by another ca.15%.The mole fraction of water in the gas phase (y Figure 11.Fugacity coefficients and fugacities of CO 2 and H 2 in the CO 2 -H 2 -NaCl brine system computed from MC simulations using eq 6. Fugacity coefficients of (a) CO 2 and (b) H 2 in the gas-rich phase compared to predictions from REFPROP 108 based on the gas-phase compositions from MC simulations depicted by dashed lines.(c,d) Fugacities of CO 2 and H 2 computed from eq 13 using fugacity coefficients from (a,b), respectively.The dashed lines are linear fits to the data points.The open symbols represent fugacities computed for CO 2 and H 2 in the respective binary systems of CO 2 -NaCl and H 2 -NaCl, respectively.The raw data of the fugacity coefficients are provided in Table S7 of the Supporting Information.The phase equilibria of H 2 -CO 2 -NaCl brine systems are investigated to compare the solubilities of CO 2 and H 2 in pure water and NaCl brine with their solubilities in binary systems: CO 2 -NaCl brine and H 2 -NaCl.The dissolution of CO 2 and H 2 2) at 323.15and 423.15K is plotted in Figure12e,f.Figure12e,f clearly shows that the water content in the gas phase at 423.15 K can be 20 to 40 times higher compared to 323.15 K. Similar to the