A Comparative Study on Transport and Interfacial Physics of H 2 /CO 2 /CH 4 1 Interacting with H 2 O and/or Silica by Molecular Dynamics Simulation 2

8 Underground H 2 storage (UHS), i.e., injecting H 2 into subsurface geological formation and its 9 withdrawal when needed, is identified as a promising solution for large-scale and long-term storage 10 of H 2 . In this study, molecular dynamics (MD) simulation was performed at a typical temperature 320 11 K with pressure up to 60 MPa to predict H 2 transport properties and H 2 -H 2 O-rock interfacial 12 properties, which are compared with those of CO 2 and CH 4 . The MD results show that the CH 4 13 profiles of property variations with pressure lie between those of H 2 and CO 2 and more comparable 14 to CO 2 . The interaction of H 2 with H 2 O/silica is much weaker than that of CH 4 and CO 2 . It is found 15 that the effect of H 2 pressure on altering the water contact angle and interfacial tension is negligible 16 under all conditions. Unlike the multi adsorption layers of the confined CO 2 and CH 4 , there is only 17 one adsorption layer of H 2 confined by silica nano-slit. The planar diffusion of H 2 in the confined 18 system is slower than that in the bulk system at pressures lower than 20 MPa. The data and findings 19 of this study will be useful for modelling the multiphase flow dynamics of UHS on reservoir scale, 20 optimizing UHS operation and assessing the performance of a cushion gas, e.g., CO 2 or CH 4 .


Introduction
Hydrogen is attracting enormous attention because it is identified as a zero-carbon energy carrier which will play a key role in achieving net-zero emissions.H2 has notable advantages such as a high energy density (141.86MJ/kg, 2-3 times higher than most hydrocarbon fuels [1]), and more importantly, clean products in combustion which can contribute significantly to decarbonize transportation, domestic heating and power generation.Although 95% of H2 in the market is still produced from fossil fuels such as the steam reforming of natural gas [2], green H2 produced by water electrolysis using renewable energy like solar/wind scales up rapidly in recent years.It is estimated that the production of low-emission H2 could reach 16-24 million tons (Mt) per year, with 9-14 Mt based on electrolysis by 2030 [3].
The storage of H2 remains a major barrier on the deployment of green H2 because the conventional surface storage facilities (e.g., high-pressure tanks) cannot meet the requirements of the capacity and time span on the scales of GWh-TWh and weeks/months, respectively, considering the fluctuation of H2 production due to the intrinsic intermittency of the renewable energy and the interseasonal discrepancies of demands (e.g., intensified demand for heat in winter) [4,5].Subsurface H2 storage in geological formations of salt caverns, saline aquifers, and depleted oil/gas fields, i.e., the so-called underground H2 storage (UHS), emerges as a promising method to achieve large-scale and long-term H2 storage.It has been proven that UHS in salt caverns is commercially feasible by an industrial-scale project (a mixture of 95% H2 and 3-4% CO2) implemented in Teesside, UK [6], whereas the progress is slow on UHS in depleted gas fields and aquifers.Salt caverns are impermeable while the other reservoirs feature porous media.
Injection of H2 into porous reservoirs displaces the original fluids, leading to complex multiphase flow processes controlled by rock properties (e.g., permeability and porosity), fluid properties (e.g., density, viscosity and diffusivity) and fluid-rock interfacial properties (e.g., wettability, interfacial tension and surface adsorption) [4,6].During injection/withdrawal cycles, cushion gas such as CH4 or CO2 is needed to work as a buffer to maintain the pressure (for H2 extraction at a desired pressure and flow rate), prevent water from entering the stored compartment, and optimize storage spaces [7,8].
In some aspects, UHS is comparable with other underground gas storage (UGS) of town gas and nature gas, as both involve compressed gases being injected, stored, and withdrawn cyclically in subsurface formations.It should be noted that UHS and UGS are different from CO2 geo-storage (CGS) in which there is no withdrawal or extraction process, and the dissolution and mineralization processes are regarded as the most secured CO2 trapping mechanisms.The dissolution and geochemical reactions of H2 would be clearly unfavorable in UHS scenarios, on the other hand.Previous extensive experiences in UGS and CGS can provide useful guides for UHS.It is, however, emphasized by Pan et al. [9,10] that direct extrapolation from these different scenarios should be avoided given that fundamental properties of H2 are significantly different from those of CO2 and CH4, as shown in Table 1.For example, H2 has a much higher diffusivity, lower solubility/density/viscosity than CH4 and CO2.H2/CO2/CH4 can switch between a gas phase and supercritical fluid at underground conditions.Meanwhile, CO2 is also likely to exist in liquid phase owing to its high critical temperature of 304.13 K.It is also found that UHS and CGS exhibit significant differences in gas saturation distribution due to their differences of gas-brine relative permeability hysteresis [10].Therefore, the transport and interfacial properties of H2/CO2/CH4 interacting with H2O and/or rock over a wide range of operating conditions should be comprehensively investigated and compared to better understand their differences in hydrodynamics and trapping mechanisms.Reliable property data is a precondition to reducing the uncertainties in site selection and cushion gas assessment, guiding future pilot and industrial-scale UHS projects.
As critically reviewed by Masoud et al. [11], a few experimental measurements were carried out recently to investigate the effects of temperature, pressure, fluid properties (gases, salinity, ion types, organic acids, etc.) and rock properties such as the surface chemistry and roughness on the contact angle (CA) and interfacial tension (IFT) for H2-brine-rock systems.CA is a quantitative measure of the wetting of a solid by a liquid, and IFT measures the force holding the surface of a particular phase together, which are both important in determining the physics of H2 flow in porous media and hydrodynamics such as rock-fluid and fluid-fluid interactions [12][13][14].Moreover, the injected gases would migrate upwards due to buoyancy until reaching the caprock, which is referred to as structural trapping [4].The interfacial interaction among the gas, water and rock also affects the storage safety and capacity as it determines the capillary entry pressure ( ! ) and the maximum column height (ℎ) of the gases sealed beneath the cap-rock [15]: where  () is the interfacial tension between the gas and water,  is the contact angle,  is the pore through radius, Δ is the gas-water density difference, and  is the gravitational constant.The data available for H2 is still sparse compared to CO2 and CH4.The wide variability and considerable inconsistencies for the measured contact angle in reported data of CO2 and CH4 [15][16][17][18][19][20] suggest that more elaborated experimental measurements (e.g., effective surface preparation, standard protocol for sample processing, distinguishing artificial contamination using advanced imaging techniques [18,19]) of H2-brine-rock systems are needed to reduce the uncertainty.Another property that is rather challenging to be experimentally measured but urgently needed is the surface adsorption and diffusivity of H2 in reservoir pores, particularly the nanopores of 5-15 nm [21], as the physical properties of gases may change dramatically due to the nanoconfinement.
Modelling methods at molecular scales such as molecular dynamics (MD) simulation provide a cost-effective and safe means of studying multiphase systems composed of hazardous substances under extreme conditions.Its capability of reducing experimental uncertainties and avoiding logistical challenges associated with flammable gases under high-temperature, high-pressure conditions makes it a valuable tool for UHS research and exploration.MD simulation has been used extensively to study CO2/CH4-brine-rock systems and predict interfacial properties of the wettability [23][24][25][26], IFT [27][28][29][30][31], surface adsorption and nanoconfined diffusion [32][33][34][35][36] for the applications of CGS or CO2 enhanced shale gas recovery.This study is motivated by the lack of H2 property data in realistic reservoir conditions and understanding of the differences in properties as H2/CO2/CH4 interacts with brine and/or rock.The novelty of the present study therefore lies in using molecular dynamics simulation to provide these important missing parameters under typical under-surface porous-reservoir conditions.
To this end, MD simulation was performed to predict and compare the properties of H2/CO2/CH4-H2O-silica systems at 320 K and pressure up to 60 MPa (the pressure can be higher than 30 MPa in the context of UHS [11]).The paper is organized as follows.Details of the molecular models, force fields, system configuration and MD setups are given in Section 2. Results of the effects of gas pressure on transport properties in the bulk phase, and interfacial properties of wettability, surface tension, surface adsorption and the diffusion coefficient of gas-water-silica multiphase systems are given in Section 3. Finally, conclusions are drawn in Section 4.

Molecular model, force field and system configuration
The molecular models and system configuration are shown in Fig 1 .Besides the molecular models of H2, CO2, CH4 and H2O, as shown in Fig 1 (a), two silica models i.e., Q2 (9.4 silanol groups/nm 2 ) and Q4 (no silanol group) developed by Heinz et al [37] are used to represent the solidphase component of two typical geological formations.A force field is needed to describe the atom interaction as a function of the position.The potential energy in this study includes the nonbonded and intramolecular energies of the bond stretch and angle bend.The dihedral and improper energies are not considered for silica [38].The potential energy is expressed as: where  4 is the dielectric permittivity of vacuum. / and  0 are the charge of particle  and  , respectively. /0 is the distance between particle  and . and  are the energy and size parameters for Lennard-Jones (LJ) potential, respectively. .and  & are the energy constants. is the bond length. 4 is the equilibrium bond length. is the angle between two bonds. 4 is the equilibrium angle.The interaction between unlike particles is modelled by the Lorentz-Berthelot combining rules:  /0 = ( // +  00 )/2 and  /0 = = //  00 > 2/" .
The Interface Force Field (IFF) developed by Wang et al. [39] and Heinz et al. [37] is used for H2 and silica; TraPPE (Transferable Potentials for Phase Equilibria Force Field) for CO2 [40]; OPLS-AA (Optimized Potentials for Liquid Simulations-All Atoms) for CH4 [41]; TIP4P for H2O molecules [42].The LJ parameters and charge of the components are listed in Table 2.These force fields, i.e., IFF, TIP4P, TraPPE and OPLS-AA, are compatible with each other and have been used in studying interfacial interactions in multiphase systems under reservoir conditions [23,24,[43][44][45].The validation of the force fields and comparison of our results with literature data will be presented in the following.

As shown in Fig 1 (b)
, a bulk gas system (S0) packed with 10 000 molecules is used to predict transport properties of gases.Three other multiphase systems are also built for studying the aforementioned interfacial properties of the contact angle in system 1 (S1, gas-liquid-solid), the interfacial tension in system 2 (S2, gas-liquid), and gas diffusion under nanoconfinement in system 3 (S3, gas-solid).There are 10 000 water molecules packed in S1 and S2.The

Simulation details
Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is used for all simulations [46], and visualization is performed in OVITO (the Open Visualization Tool) [47].The initial systems are built by PACKMOL [48].For S0 and S2, periodical boundary conditions are implemented in all directions.For S1 and S3 with solid silica slabs, periodical boundary conditions are implemented in x and y directions, while the nonperiodic fixed ppf boundary condition is used in the z direction.With ppf, a vacuum space, the length of which is three times the z-edge length, is virtually inserted in the z direction to avoid the artifacts generated by the interactions between the real and replicated systems [49,50].The cut-off distance for intermolecular interaction is 1.6 nm in real space.The particle-particle particle-mesh (PPPM) algorithm is used for long-range electrostatic interactions in the reciprocal k-space with an accuracy of 1×10 -5 .
All simulations are performed with a time step of 1 fs for two consecutive processes, i.e., equilibrium run for the system to reach stable conditions and production run for collecting data, as shown in Fig 1 (b).In S0, an isothermal-isobaric ensemble (NPT where N: number of molecules, P: pressure, T: temperature) is used to maintain the temperature and pressure during the equilibrium run, followed by a canonical ensemble (NVT where V: volume) in the production run to maintain a constant volume and the temperature.In NPT and NVT simulations, the damping factors for the Nosé-Hoover thermostat and barostat are 100 and 1000 times of the time step, respectively.For multiphase systems, the berendsen thermostat is used to maintain the temperature with a damping factor of 100.For S1 and S3, a 5-Å layer of the outmost silica molecules is kept fixed.respectively.Viscosity is calculated using the Green-Kubo method [51,52] in NVT ensemble by:

As shown in
where  is the volume,  @ is the Boltzmann constant,  is the temperature,  => represents the  components of the pressure tensor,  and  are any two of the x, y or z Cartesian coordinates, and á•••ñ indicates ensemble average.According to the work of Nie et al. [53], the correlation length is set to be 10 ps with a sampling interval of 1 fs for pressure tensor.The final value of the viscosity in In Fig 2 (b) and (c), both the density and viscosity of the gases agree well with the NIST data, indicating that the equilibrium molecular dynamics (EMD) simulation can be employed to predict the transport properties of binary/ternary mixtures over a wide range of reservoir conditions.H2 densities are underestimated slightly, and the maximum absolute relative deviation for H2, CO2 and CH4 is 5.1%, 2.4% and 4.2%, respectively.H2 and CH4 densities increase gradually with pressure passing across the critical pressure, whereas there is a dramatic change for CO2 from 5 to 20 MPa.The extremely low density of H2 makes it difficult to displace H2O in under-surface porous media.The substantial density difference between H2 and CO2/CH4 indicates that a strong gravity override and density segregation will occur during UHS operation, resulting in great upward movement of H2 towards the caprock than CH4 and CO2.Regarding the viscosity, the maximum deviation for H2 occurs at 5 MPa, which is 5.3% lower than NIST data, and the averaged absolute deviation of CO2 and CH4 viscosities are 12.1% and 8.0%, respectively.Compared with CO2 and CH4, the H2 viscosity is only slightly influenced by pressure.A low viscosity of H2 can facilitate its injection and extraction, but can also results in undesired viscous fingering [5,54].

Wettability of gas-water-silica S1 system
Wetting mechanics is the synergistic effects between the interactions of gas-water, gas-solid and water-solid phases.The effects of the gas temperature and pressure on wettability can be quantified by the alternation of the water contact angle.The droplet morphology and the contact angles are demonstrated in Fig 3 .Firstly, the time evolution of the center of mass (COM) of the water droplet on Q2 silica in the z direction is evaluated in Fig 3 (a).The decrease of the COM-z of the droplet is because the silanol groups on Q2 surface can form the hydrogen bonding with water, which leads to the gradually spreading of water molecules over the silica surface.The water droplet spread much more slowly when enveloped by CO2 than by CH4, while H2 can hardly alter the spreading process.CO2 has high affinity with silica because of the strong electrostatic interaction as CO2 has quadrupole polarizability.There is only vdW intermolecular interactions between H2 and silica.The morphology of the water droplet is stable as the COM would not change with time after 4 ns, and the snapshots of the systems at 20 MPa are shown in Fig 3 (b).The water droplet on Q2 is flattened, while the droplet can almost be detached from the Q4 surface when enveloped by CO2.As shown in Fig 3 (c), the contact angle is computed by the 2D density profile of the water droplet, which is averaged over 4-6 ns using bins of 0.5´0.5 Å.A circular profile is fitted at the gasliquid interface where the iso-density is in the range of 0.2-0.4g/ml [26,55].The contact angle is determined by: where  and  are the coordinates of the center of the fitted circle in x and z directions, respectively. is the radius. 4 is the height of the contact plane, which is the average position of the outmost silica atoms [56], i.e., 19.3 Å for Q2 and 26.5 Å for Q4.
The contact angles of H2/CO2/CH4-H2O-silica systems are summarized in Fig 4 .For comparison, representative experimental data of a gas-water/brine-quartz system is collected and plotted.The tilted plate technique was used in experiments, and two contact angles were generated, i.e., advancing contact angle ( C ) and receding contact angle ( D ).The equilibrium contact angle ( E ) can be obtained by [11]:  E = arccos X The contact angle of water without an ambient gas (i.e., 0 MPa) on Q2 silica at 320 K is 27.3°, which indicates the silica is strongly hydrophilic; while the water contact angle on Q4 silica is 106.2°,indicating the silica is strongly hydrophobic.At the interface, a CA smaller than 90° corresponds to high wettability or hydrophilicity, whereas a contact angle larger than 90° signifies low wettability or hydrophobicity.The results are consistent with the earlier reported MD simulation results of 103.8° for a H2O (SPC/Fw model)-Q4 system at 300 K [50].The contact angle in the CO2-H2O-Q4 system is 141.9° at 10 MPa, which is close to 146.7° at 10.5 MPa and 318 K reported by Chen et al. [23], who used the flexible SPC water and EPM2 CO2 models.
In Fig 4 , CA follows the order of  JK + >  JL , >  L + in both experiments and MD simulations under the same thermodynamic condition and with the same silica subtract.This implies that CO2 has higher wettability and therefore is more favorable than CH4 as the cushion gas for withdrawing H2.At 323 K, the measured CA is increased by 16°, 21°, and 16° for H2, CO2 and CH4, respectively, when the pressure is increased from 5 MPa to 20 MPa.The experimental results also reveal that CA increases with the temperature for CO2/H2-water/brine-quartz systems, but not for CH4 where the CA has the same value at 300 K and 320 K.The effects of the temperature on wettability depend on minerals, and the mechanism is not yet fully understood [20,60,61].Moreover, different measurement procedures and sample cleaning processes in experiments could result in the uncertainties of CA and the inconsistency of its temperature dependence [24].It has been widely accepted that CA increases with the CO2 pressure for various minerals, which is attributed to an increased gas molecular density and enhanced intermolecular interaction between the gas and rock surface [60,62,63].CO2-deionised water-quartz [58]; Reproduced with permission from J. Chem.Thermodyn.93, 416 (2016).Copyright 2016 Elsevier.CH4-brine (1.5 wt.% NaCl)-quartz [59].Reproduced with permission from Energy Fuels 33, 788 (2019).Copyright 2019 American Chemical Society.
As for H2, a linear correlation between cos and gas density () is well fitted for H2-H2Omica/clay/basalt systems, with  " (goodness of fit) higher than 90% at 330 K and 0-20 MPa [11].But it should be mentioned that some experiments found that CA is not affected by the temperature and pressure in H2/CH4-water/brine-sandstone systems when the captive bubble method was used [64,65].In MD simulations using a sessile water droplet, CA is increased by 21° and 50° when the CO2 pressure is increased to 10 MPa and 20 MPa in Q2 and Q4 silica systems, respectively.The CA is only increased by 4° and 7° in Q2 and Q4 silica systems, respectively, when the CH4 pressure is increased to 30 MPa, which is much less noticeable than in experiments.Effects of the H2 pressure on CA show a different trend between MD and experimental data, and there is almost no change to CA in the MD simulations for both Q2 and Q4 silica in conditions up to 60 MPa.It can be seen from Fig 4 that for hydrophobic Q4 silica, H2 virtually does not change the contact angle of water at 0 MPa, i.e., under no impact of any gas components, even when the pressure is increased to 60 MPa.For hydrophilic Q2 silica, H2 tends to increase the water contact angle in vacuum, but only slightly compared to CH4 and CO2.And the slight increase is also independent of the pressure, similar to what is observed for Q4 silica.Due to a much weaker intermolecular interaction with the rock surface together with a low molecular mass and thus volume, H2's capability of adsorbing on to either the Q2 or Q4 silica surface is low, even in supercritical fluid state as the system pressure far exceeds its critical one (the system temperature is much higher than the critical one of all the three gases).
According to [60], cos  + 1 ∝  M,L + K A2 , i.e., the contact angle is largely determined by the surface tension force between water and the gas in such a system.As shown in Fig 5 (c), the H2/H2O surface tension is independent of pressure, which is also corroborated by experimental data shown in the figure.Therefore, the effects of H2 on wettability is low, and the water contact angle is independent of the H2 pressure, either with a hydrophilic or hydrophobic silica substrate.Moreover, although considerable discrepancy or inconsistency of experimental CA data for CO2-water/brine-rock/shale [16,18,66] and H2-water/brine-rock systems have been seen, our MD results are consistent with literatures [67,68].The CA at above 30 MPa is almost constant for all gases in MD simulation, and a further increase of the pressure has negligible effects on wettability.At  > 30 MPa, the adsorbed CH4 and CO2 on silica surface would not change notably with the pressure, indicated by the peak densities shown in Fig 6 .In addition, the contact angle is largely determined by and inversely proportional to the surface tension force between water and the gas in such a system, as stated above.As shown in Fig 5 (c), the surface tension of the three gases, both experimental and MD data, shows little dependence on the pressure.Therefore, the CA has little dependence on pressures higher than 30 MPa.
For CGS applications, most previous MD modelling studies of CO2-water/brine-rock systems is at pressures below 25 MPa [23,24,26,54,55,69,70].In a recent study by Zhou et al. [49], it is found that contact angles keep increasing gradually with a reduced growth rate in a CO2-H2O-kerogen system, changing from strongly H2O-wet (CA=60.4°at 0 MPa) to strongly CO2-wet (CA=180° at 44 MPa), which implies that the effects of a higher pressure on H2 wettability over an organic subtract should be investigated, considering coal seams are also a feasible reservoir type for UHS [11,71].

Interfacial properties of gas-water S2 system
The interfacial properties of the gas-water system S2 are evaluated by the gas/water density, surface excess () and IFT (), as demonstrated in Fig 5 .The Gibbs dividing surface (GDS) is used to identify the position of the interface where the surface excess of water is zero from the density profile of water [30,72], as shown in Fig 5 (a).The gases have a positive surface activity and can accumulate at the gas-water interface, where the gas density is higher than that away from the interfacial region.The positions of the peak values of CH4 and H2 are closer to GDS than CO2.
The surface excess of gas  relative to water ( / ) ) is used to quantify the gas adsorption tendency to water at the interface, expressed as [30,72]: where  is the total number of gas molecules,  is the number density of the gas, I and II denote, respectively, the gas-and water-rich bulk phases distinguished by the GDS,  is the volume of phase I or II, and  is the area of interfaces.
The surface tension is computed by [30,72]: where  == is a diagonal component of the pressure tensor, and the pre-factor of 1/2 considers the existence of two interfaces in the simulation box.
The surface excess follows the order of ) under all conditions, as shown in Fig 5 (b). JK + ) increases rapidly to the peak value of 0.156 mole/Å 2 at 10 MPa, then drops to 0.05 mole/Å 2 at 20 MPa, followed by a slow decrease with the pressure. JL , ) increases gradually with the pressure to the peak value at 20 MPa, and then decreases. L + ) , on the other hand, increases monotonously with the pressure.The IFT results of MD simulation agree well experimental data for both CO2 and CH4, as shown in Fig 5 (c).The H2-water IFT is underestimated by ~11% on average as compared to the experimental values.Contrary to the surface excess, IFT follows the order of  JK + <  JL , <  L + under all conditions.The CH4-water IFT decreases gradually with the pressure with a reduced rate.There is a clear turning point for CO2-water IFT at ~10 MPa where the IFT decreases remarkably by 42%, followed by a slight decrease with the pressure.In the experiment, the H2-water IFT decreases slightly with the pressure at a reduction rate of -0.01 mN•m -1 /MPa on average for temperatures ranging from 298 K to 423 K [73], while there is no pressure dependence observed in MD simulation at 320 K. H2 can hardly alter the IFT because the IFT-pressure correlation is associated with the density difference between the gas and water [11].Moreover, the intermolecular vdW interaction contributes scarcely to the gas-water interaction at GDS from the energy decomposition perspective [79].

Nanoconfined adsorption and diffusion behavior of gas-silica S3 system
The gas adsorption in nanoslits can be characterized by the gas density profile, as demonstrated in Fig 6 .The simulation time is 2 ns for all cases.The density profile is computed by: where D is the bin size 0.2 Å,  is the molecule weight, and  is the number of gas molecules in the bin.The gas density varies with the distance away from the outmost atoms of the silica subtract.The silica nano-slit can be divided into three regions according to the density profile, i.e., adsorption zone, transition zone, and bulk zone [80].Gas molecules accumulate in the region close to the silica surface forming adsorption layers, featuring peak values in the density profile.The density in the adsorption zone decreases further away from the silica surface owing to a diminishing gas-silica interaction.Further away from the transition zone and approaching the slit center, the bulk-fluid region would be formed where surface effects are negligible, molecules can move freely, and the density does not change with position anymore.The transition zone of CO2 at 5 MPa and 10 MPa is 0.5-1.0nm, which is wider than that of CH4 and H2, both less than 0.3 nm.The thickness of the first adsorption layer ( \] ) of H2 in both Q2 and Q4 systems is 0.4 nm, which is consistent with the theoretical prediction 0.3-0.5 nm [81].
The density profile depends strongly on the pressure and gas type, but only slightly on the silica type of Q2/Q4.For H2 in both Q2 and Q4 systems, there is only one adsorption layer at pressures below 50 MPa.The second layer can be discerned at 60 MPa, which is only 1.03 times denser than in the bulk region.The second adsorption layers are formed clearly at pressures above 10 MPa and 20 MPa for CO2 and CH4, respectively, in both Q2 and Q4 systems.There is even a third adsorption layer of CO2 and CH4, which can be seen at pressures above 30 MPa and 50 MPa, respectively.The density of H2 in the adsorption zone increases linearly with the pressure in both Q2 and Q4 systems, with  " of the linear fitting equal to 99.1% and 98.6%, respectively.The peak values of the first adsorption density in Q2 is 1.44 times higher than that of the bulk density obtained from the linear fitting of  ^ECU ~  RSTU , and 1.32 for H2-Q4 system.For CH4 and CO2 in both Q2 and Q4 systems, the density increases with a declined growth rate for both the first and second adsorption layers.
The effects of nanoconfinement on mobility properties of the gases are compared by the mean square displacement (MSD) and diffusion coefficient (), as shown in MSD is a measure of the deviation of the particle position with respect to its initial position over time, expressed as: where  is the number of particles, and vector  (/) () is the position of the  `a particle at time .
As shown in Fig 7 (a), MSD's of the gases in x, y and z directions overlap with each other due to the isotropy intrinsic of the bulk system.MSD correlates linearly with time as gas molecules undergo Brownian motion, featuring normal diffusion [82].The diffusion coefficient can be extracted from the slope of the MSD using 〈 " ()〉 = 2, where  is the number of dimensions.According to Zhong et al. [83], a system of more than 500 gas molecules is large enough to avoid the finite-size effect on self-diffusion.
The MSD of H2 has a much larger slope than that of CH4 and CO2 at 20 MPa.For gases in the nanoconfined system S3, the MSD's in the planar x and y directions also overlap with each other and increase linearly with time.However, in the confined z direction, gas molecules undergo sub-diffusive motion, featuring anomalous diffusion [82,84].As shown in Fig 7 (b), the MSD in the z direction increases nonlinearly with a reduced growth rate until it reaches the plateau value at ~15-17 nm 2 , and the H2 case reaches a plateau much faster than CH4 and CO2.This is consistent with the theoretical value of 16.67 nm 2 , as analyzed by Wang et al. [85] that MSD approaches the value of  8 " /6 after the time is longer than a characteristic time scale.The diffusion in the normal direction would affect the swap frequency of the gases between the adsorption and bulk regions.The planar motion is the primary reason causing the leakage of gases in underground porous formation [6].The effect of space confinement is studied by comparing the difference between the planar diffusion coefficient ( || ) in the nanoconfined system and the diffusion coefficient in the bulk system ( RSTU ), as shown in Fig 8 .The diffusion coefficient follows the order of:  L + >  JL , >  JK + in both the bulk and confined systems under all conditions.The diffusion coefficient of a gas decreases with the pressure due to the increased intermolecular friction according to gas-kinetic theory.The  RSTU of H2/CO2/CH4 decreases drastically from 5 MPa to 20 MPa by a factor of 3.94, 10.43 and 4.36, respectively, followed by gradual decreasing with the pressure.The  || of CO2 and CH4 in the confined system has a similar trend to  RSTU , while  || of H2 decreases linearly with the pressure in the confined system.The silica type has negligible effects on planar diffusion of the gases.At 5, 10 and 20 MPa, the confinement effect on H2 diffusion is evident as the ratios of  || / RSTU are ~0.23,0.46 and 0.70, respectively.At higher pressures of 30-60 MPa, the confinement effect weakens as the ratio ranges in 0.85-0.99 for all gases.It should be mentioned that the pore size also plays a crucial role in gas adsorption and diffusion.For example, the bulk zone would disappear when the slit is narrower than the critical size, e.g., 2 nm for CH4 in Montmorillonite slit pore at 323 K [86].

Conclusions
The reliable data of transport and interfacial physics of H2/CO2/CH4 interacting with H2O and/or rock is essential for understanding H2 injection, storage, and withdrawal mechanisms, particularly when CO2 or CH4 is used as the cushion gas.Molecular dynamics simulation was performed at 320 K to study the effects of H2 pressure (up to 60 MPa) on transport properties (density, viscosity, and diffusion coefficient) in the bulk system, and interfacial properties (wettability, interfacial tension, surface adsorption, and confined diffusion) in contact with water and/or silica.The H2 properties in the three systems are compared with those of CO2 and CH4.Two silica molecular models are built, i.e., hydrophilic Q2 and hydrophobic Q4.The conclusions are summarized as follows: The density and viscosity predicted by equilibrium molecular dynamics agree well with experimental data.The density and viscosity of H2/CH4 increase gradually with the pressure, while there is a drastic change for CO2 as the pressure crosses the critical value.
The contact angles of water on both the Q2 and Q4 silica increase with the gas pressure of CH4 and CO2 until reaching a plateau at 10-30 MPa.H2 has a much lower wettability than CO2 and can barely alter the water contact angle even at 60 MPa.
Similar to CO2 and CH4, H2 molecules can also accumulate at the interface with water.The surface excess of H2 increases monotonically with the pressure, but at a much lower rate than that of CO2.Effects of the H2 pressure on the interfacial tension are negligible, while CH4 can slightly reduce the IFT.The interfacial tension of CO2-water decreases significantly when the pressure is increased to 10 MPa but can remain almost unchanged at pressures above 20 MPa.
Confined by the nano-slit of 10 nm, the first adsorption layer of H2 has a thickness of 0.4 nm.The second adsorption layer of H2, CO2 and CH4 emerges at 50, 10 and 20 MPa, respectively.The density of H2 in the adsorption zone increases linearly with the pressure.Gas molecules undergo subdiffusive motion in the confined direction, and the planar motion is still in the normal diffusion regime.H2 planar diffusion is slowed down evidently due to the nanoconfinement effects at pressures lower than 20 MPa.Silanol functional groups on the silica surface have insignificant effects on gas diffusion.These fundamental data and findings can be used as input parameters for predicting the capillary pressure, gas column height and the IFT of the porous media in the context of UHS.Moreover, the MD systems configured in this study can be extended to predict properties of the binary gas mixtures as there is a mixing zone in UHS between H2 and the cushion gas.
Fig 1.(a) Molecular models of gases (H2, CO2, and CH4), H2O, and the silica unit cells (Q2 and Q4).The size of the Q2 and Q4 unit cells in x, y and z is (34.7,34.3, 17.8) Å and (33.4,34.8, 23.9) Å, respectively; (b) Initial configurations and equilibrated H2 systems at 320 K and 20 MPa after 4 ns simulation:  6 and  7 are the box edge lengths in x and y directions, respectively;  8 of S1 and S3 is the width of the silica nano-slit.The thickness of the H2O film in S2 is 8.4 nm.

Fig 2 ,
the densities and viscosities obtained from bulk MD simulations are compared with the data from NIST Chemistry Webbook [22].The converging process of the viscosity is shown in Fig 2 (a), and the final values of the density and viscosity are shown in Fig 2 (b) and (c),

Fig 2 (
c) is the averaged results of the last 2 ns data of Fig 2 (a).

Fig 2 .
Fig 2. (a) Convergence of gas viscosity using GK method; (b) Effects of pressure on density; (c) Effects of pressure on viscosity.

CH 4 Fig 3 .
Fig 3. (a) Time evolution of H2O droplet COM-z at 20 MPa; (b) System snapshots at 20 MPa after 6 ns; (c) Density contour plots of H2O droplet in H2-H2O-silica S1 systems.The black dashed lines show the fitted circles.

Fig 5 .
Fig 5. (a) Density profiles of water and reduced density profiles (/ RSTU ) of gases at 20 MPa, where the black dotted line indicates the GDS of H2-H2O located at  = 41.84Å; (b) Effects of pressure on surface excess adsorption at 1 MPa and 5-60 MPa; (c) Effects of pressure on IFT.

Fig 6 .
Fig 6.Gas density profile along the direction normal to the silica surface at different pressures.(a) Gas confined by Q2 silica; (b) Gas confined by Q4 silica.The dashed lines indicate the thickness of the first adsorption layer.

Fig 7 .
Fig 7. (a) MSD of gases in x, y and z directions in the bulk-gas system S0 at 20 MPa; (b) MSD of gases in Q2 nanoslit at 20 MPa.The dashed lines in (b) indicate the boundaries of the plateaus.

Table 2 .
LJ parameters and charge q of silica, H2, CO2, CH4, and H2O.M is the additional virtual 158 particle of the TIP4P H2O model.159

Table 3 .
Number of gas molecules in multiphase systems.