Adsorption of hydrogen, methane, CO₂ and their binary mixtures in silicalite-1: role of pore characteristics revealed by molecular simulations

Abstract

Physisorption in nanoporous materials is an important alternative for storing hydrogen fuel. Use of methane and CO₂ as cushion gases can aid in keeping hydrogen stored in geological repositories. For this, an understanding of the coadsorption of these gases and the underlying mechanisms that are determined by the pore characteristics of the repository rocks is essential. Here we use grand canonical Monte Carlo (GCMC) simulations to understand coadsorption of hydrogen, methane, CO₂ and their binary mixtures in microporous silicalite-1. To understand the effects of pore characteristics like surface area, connectivity, and tortuosity, we progressively block some channel-like pores of the adsorbent by filling it with immobile organic material (represented by methane). At 100 atm in unmodified silicalite the adsorption amounts of hydrogen, methane and CO₂ in pure state are respectively (1.08 ± 0.11) mmol g⁻¹; (2.86 ± 0.06) mmol g⁻¹ and (3.52 ± 0.03) mmol g⁻¹. Both CO₂ and methane are found to suppress hydrogen adsorption from an equimolar mixture, while presence of hydrogen in the mixture has no discernible impact on the adsorption of these gases. Suppression of hydrogen adsorption by CO₂ is stronger by an order of magnitude. This is because the adsorption energy of CO₂ in silicalite is about four times that for hydrogen and 1.5 times that for methane. Adsorption of all gases in pure state as well as the carbon fluids in the mixtures exhibits clear dependence on the ratio of surface area to volume of the pores. Information obtained in this study can guide the design of a future study that recreates the UHS scenario of a cushion gas being pushed on top of hydrogen present in a reservoir.

---

1. Introduction

Growing population and industry across the world demand a large supply of energy necessitating the search for alternative sources of energy.^1–3 Hydrogen, with a high energy density, is one such renewable source of energy.^4,5 While the thermodynamical properties of hydrogen make it an attractive source of energy, storing it in large quantities for use at industrial scales remains a challenge.⁶

Underground hydrogen storage (UHS) is a promising alternative that uses sub-surface salt caverns or depleted oil and gas reservoirs to store hydrogen.^7,8 When hydrogen is retrieved from these reservoirs a substitute gas is needed to maintain the reservoir pressure. In addition, such a gas also prevents water intrusion and provides a permanent cushion to maintain the pressure and is therefore called a cushion gas.⁹ Non-reactive gases like methane or CO₂ can be used as cushion gases.^9,10 The storage capacity of the UHS sites is characterized by the porous structure of the formation rock. The important parameters characterizing this structure are permeability, pore size distribution, geometry and connectivity of the pores.¹¹ Understanding the effects of pore characteristics on the mixture of a potential cushion gas with hydrogen is therefore important.

While the UHS formations exhibit a wide range of pore sizes, some UHS sites, especially shale gas reservoirs can have small micropores (<2 nm).¹² In some previous works we have explored the effects of pore characteristics like shape, connectivity, surface area and volume of the pores on the adsorption and dynamics of fluids in sub-nanometer pores of siliceous zeolites.^13–17 These studies revealed important trends observed in the properties of confined fluids and their mixtures. In particular, the use of silicalite – an all-silica zeolite of MFI type pore topology – provided an important combination of different pore geometries in the same system. Further, selectively blocking some pores in this material helped us systematically vary the pore surface area, volume and the degree of connectivity.

While the adsorption of hydrogen, methane and CO₂ has been reported in porous media including kerogen¹⁸ and organic-rich rocks,¹⁹ to our knowledge, mixture coadsorption of these fluids in sub-nanometer pores remains unexplored. Further, no systematic investigation addressing the effects of pore characteristics like tortuosity, connectivity, surface area and volume of the pores has been reported. To address this knowledge gap, in the current work we discuss these effects on mixtures of hydrogen with two potential cushion gases – methane and CO₂. This serves as an exploratory study with the objective of understanding the behavior of such mixtures in UHS formations. The focus is on the interplay of the effects of the pore characteristics and the presence of the cushion gases on the adsorption and distribution of hydrogen in the sub-nanometer pores of silica. In the next section, we provide details on the making of adsorbent models; the force-fields used to represent different entities; pore characterization and the simulation methodology. After this, salient results are presented in Section 3 and discussed in Section 4 in the context of current literature and potential implications. Finally, we summarize the main conclusions in Section 5.

2. Methods

2.1 Adsorbents

All adsorbent models were based on silicalite with MFI geometry characterized by a network of channel-like pores ∼0.55 nm in diameter.²⁰ Straight channels with an elliptical cross-section run along the Cartesian direction Y and are intersected by tortuous channels having a zigzag pore axis lying in the Cartesian XZ plane. The intersections provide a slightly larger pore space in an ellipsoidal shape. The atomic coordinates of silicalite provided by van Koningsveld et al. are used in this work.²¹ A unit cell was replicated 2 × 2 × 3 times using VESTA²² to result in a cuboidal simulation cell of sides ∼4 nm. This constituted the unmodified silicalite with all the pores (straight as well as zigzag) open and is henceforth referred to as S4Z4 (see Fig. 1). To vary the pore characteristics, some of the pores were selectively blocked by immobile methane molecules. To do this, a grand canonical Monte Carlo simulation was carried out for adsorbing methane in the pores of S4Z4 at 200 K and 200 atm pressure. This resulted in saturation loading of all the pores filled with methane to maximum capacity. After this, some methane molecules were removed in such a way as to selectively free some of the straight or zigzag channels. S4Z4 has a total of 12 each of straight and zigzag channels. Out of these 12, some were freed in steps of three each. Thus, adsorbents with 12, 9, 6, 3 or 0 of straight or zigzag channels were left free. This resulted in 12 different adsorbent models with different combinations of straight/zigzag channels open. In general, the adsorbents are named SnZm (n, m = 0, 1, 2, 3 or 4) with n and m denoting, respectively, the fraction (out of 4) of straight and zigzag channels open. For example, S4Z0 has all straight channels open, all zigzag channels blocked, S0Z4 has all straight channels blocked and all zigzag channels open, and S2Z2 has half of each type of straight and zigzag channels open.

Fig. 1 Simulated adsorbents (a and b) and absorbates (c). (a) Unmodified silicalite (S4Z4) in X–Z plane. Red and blue spheres represent oxygen and silicon atoms respectively. Elliptical cross-sections of straight channel-like pores running along the Y direction (perpendicular to the plane of the figure) can be seen. (b) S3Z4 with three out of 12 straight channels blocked with methane molecules (grey spheres), leaving three fourths of the straight channels open. (c) Models of the three adsorbates – hydrogen and methane in the united atom formalism represented by a structureless particle; and CO₂ in all-atom model with two oxygen atoms connected to a central carbon atom.

2.2 Force-fields

This work used a combination of TraPPE^23,24 to represent the carbon bearing fluids; ClayFF²⁵ to represent the adsorbent and Buch force-field²⁶ to represent hydrogen molecule. Methane belonging to the adsorbent as an immobile pore-blocking entity as well as the adsorbate fluid methane were both represented in the same united atom model. In this model all hydrogen atoms and carbon atom are joined together and are represented by a structureless particle interacting with other atoms solely with the van der Waals interaction.²³ Similarly, the hydrogen molecule is represented by a single structureless point particle like entity in the Buch force-field that interacts with other atoms with van der Waals interactions alone.²⁶ In a previous publication,²⁷ it was found that this force-field provides good computational economy compared to the force-field proposed by Darkrim and Levesque²⁸ that incorporates the quadrupole moment of the hydrogen molecule. Besides, as we shall see in the results section, this force-field provides good agreement with experimental data. For CO₂ an all atom TraPPE force-field was used.²⁴ The force-field parameters used in this work for all atoms and molecules are detailed in the SI.

2.3 Simulations

Grand canonical Monte Carlo simulations were carried out using DL-Monte.²⁹ The starting configuration had one molecule of the adsorbate fluid placed in an open channel of the adsorbent. As the simulation progressed, adsorbate molecules were inserted/deleted, translated or rotated with respective probabilities of 0.5, 0.25 and 0.25 in case of CO₂. For the case of hydrogen or methane as the adsorbate, the molecules were inserted/deleted, or translated with respective probabilities of 0.5 and 0.5. In all simulations, all silicalite atoms including the blocker methane molecules (not to be confused with adsorbate methane molecules) were kept rigid with no changes involved with their positions throughout the simulation. For comparison with the experimental data available for unmodified silicalite (S4Z4) simulations were carried out with the temperature and pressure conditions reported for the experimental data.^29,30,31 All other simulations were carried out at a temperature of 298 K and partial pressure (henceforth referred as pressure for brevity) of the adsorbate fluid ranging between 0.05 and 100 atm. For each of the 12 adsorbents, simulations with these environmental conditions were carried out for adsorption of either hydrogen, methane or CO₂ in pure state, or a mixture of hydrogen with either methane or CO₂. For pure fluid adsorption, a total of 2 million Monte Carlo steps were simulated, whereas for the mixture adsorption, 4 million steps were simulated. Simulations for all adsorbent–adsorbate combinations were started with the lowest pressure of 0.05 atm. All subsequent simulations at progressively higher pressures were carried out using the final configuration obtained in the simulation of the preceding pressure value. This ascertained a quick achievement of equilibration. For all simulations, it was observed that the number of adsorbed molecules started increasing and reached a stable plateau signifying equilibrium before 500 000 steps for the pure adsorbate and 2 million steps for the mixture adsorbates. Data corresponding to the initial 500 000 and 2 million steps were therefore discarded from these simulations and averages for different quantities calculated from the subsequent 1 500 000 and 2 million steps respectively for pure adsorbate and mixture adsorbate cases. Mixture adsorption in all 12 adsorbents was carried out with an equimolar composition. In addition, 13 simulations of the mixture adsorption with the feed mixture composition varying between 10 : 90 and 98 : 2 (hydrogen : carbon fluid) were also carried out for each of the three adsorbents – S4Z4, S4Z0 and S0Z4 at a total gas pressure of 100 atm.

2.4 Pore characterization

Pore networks in the 12 adsorbent models were characterized in terms of the accessible surface area and pore volume using Zeo++0.3.³² A small spherical test particle of radius 1.2 Å (smaller than the kinetic radius of all the three adsorbate molecules) was used to efficiently probe the surface area and pore volume with reasonable accuracy. The 12 adsorbents exhibit a range of pore surface area (S), volume (V) and their ratio, with S ranging between 456 and 784 m² g⁻¹; V ranging between 0.03 and 0.06 cm³ g⁻¹ and the ratio S/V ranging between 1.4 and 1.7 (×10¹⁰ m⁻¹).¹⁷ Further, the number of pore connections range between 0 (S0Z4 and S4Z0) and 48 (S4Z4) while the number of tortuous (zigzag) pores as a fraction of total range between 0 and 100% with 11 distinct values.

3. Results

3.1 Validation of the simulation data

Fig. 2 shows a comparison of the simulated adsorption isotherms with the low-pressure experimental data reported in literature. As the experimental data are reported at different temperatures (hydrogen at 305 K,³⁰ methane and CO₂ at 308 K (ref. 30 and 31)), we carried out simulations at these temperatures specifically for a fair comparison. Two different sets of experimental data are available for comparison of CO₂ and CH₄ simulation data obtained here. As can be seen the comparison is good for all the three fluids, validating the accuracy provided by the force-fields used in the simulations.

Fig. 2 Comparison between the simulated pure gas adsorption isotherms of hydrogen (left); and CO₂ and CH₄ (right) in silicalite; and the corresponding experimental data reported in the literature. The experimental data are obtained from Golden and Sircar³⁰ and Sun et al.³¹

3.2 Comparative adsorption isotherms

Fig. 3 shows the adsorption isotherms of hydrogen (top) and the carbon fluids CO₂ and CH₄ (bottom) in pure state and binary mixture with hydrogen for the three representative adsorbents – S4Z4 (unmodified silicalite), S4Z0 (all straight channels open; zigzag channels blocked) and S0Z4 (all straight channels blocked; zigzag channels open). The adsorption isotherms of the remaining 9 adsorbents can be found in the SI, Fig. S1 and S2. Hydrogen adsorption in silicalite is linear up to about 30 atm, beyond which the slope of n_ads vs. P begins to decrease. This low-pressure regime where n_ads vs. P is linear extends only for pressures below 1 atm for the carbon bearing fluids (note: the linear regime for carbon fluids is not visible in Fig. 3 because of the logarithmic scale used for the pressure axis in these cases). At all pressures, the amount of fluids adsorbed in all three adsorbents follows the order CO₂ > CH₄ > H₂. While at higher pressures the amount of hydrogen adsorbed is lower but comparable with that of carbon fluids, at lower pressures this difference gets amplified. Comparing the adsorption amounts of a fluid in the pure state vs. the mixture it can be seen that while co-adsorption with carbon fluids significantly impacts hydrogen adsorption, the presence of hydrogen makes no discernible difference to the adsorption amounts of the carbon fluids. For example, in unmodified silicalite (S4Z4) the adsorption amounts of hydrogen, methane and CO₂ in pure state at 100 atm are respectively (1.08 ± 0.11) mmol g⁻¹; (2.86 ± 0.06) mmol g⁻¹ and (3.52 ± 0.03) mmol g⁻¹. In the mixture of hydrogen with methane the amount of hydrogen adsorbed at 100 atm in S4Z4 reduces to (0.12 ± 0.04) mmol g⁻¹; while that in the mixture with CO₂ is reduced to (0.02 ± 0.01) mmol g⁻¹. In contrast the adsorption amounts of methane and CO₂ in mixture with hydrogen at 100 atm in S4Z4 remains almost unchanged at (2.76 ± 0.05) mmol g⁻¹ and (3.63 ± 0.03) mmol g⁻¹.

Fig. 3 Simulated adsorption isotherms of hydrogen (top panels) and CO₂ and CH₄ (bottom panels) in three representative adsorbents based on silicalite. Data corresponding to both pure fluid adsorption and one component in the mixture adsorption is shown for comparison. Note that the pressure axis in the top panels is in linear scale but in the bottom panels it is in logarithmic scale.

Selectively blocking pores of silicalite results in a lowering of the adsorption amounts for all fluids. In particular, blocking of all straight channels (S0Z4) results in a higher reduction in adsorption amounts compared to when all zigzag channels are blocked (S4Z0). A similar reduction in adsorption amounts by blocking of pores was also seen earlier for ethane and CO₂.^14,16 At 100 atm the adsorption amounts of hydrogen, methane and CO₂ in pure state in S0Z4 are respectively (0.67 ± 0.08) mmol g⁻¹; (1.73 ± 0.04) mmol g⁻¹ and (2.17 ± 0.03) mmol g⁻¹. The corresponding numbers for S4Z0 are (0.76 ± 0.08) mmol g⁻¹; (2.01 ± 0.05) mmol g⁻¹ and (2.49 ± 0.05) mmol g⁻¹.

While the Buch force-field used to model hydrogen in this work is validated against experimental data as shown in Fig. 2, it does not incorporate the quadrupole moment of the hydrogen molecule. Could ignoring the quadrupole moment of hydrogen result in underestimating its coadsorption with the carbon fluids as seen in Fig. 3? To investigate this, we carried out a simulation of hydrogen coadsorption with the carbon fluids in silicalite using the force-field proposed by Darkrim and Levesque (DL).²⁸ This force-field models the hydrogen molecule as a three-site structure with partial charges attached on two hydrogen atoms and the center of mass (details in the SI). As shown in Fig. S3 of the SI, the adsorption amounts of hydrogen and the carbon fluids in the mixture using DL force-field are very similar to that obtained using Buch force-field.

To understand the strong suppression of hydrogen adsorption by the carbon fluids, a comparison of the adsorption energies could be useful. For this, Clausius–Clapeyron relation in the following form is used.^29,33

Here Q_ads is the isosteric heat of adsorption, R is the ideal gas constant (8.314 J mol⁻¹ K⁻¹), P is the pressure that corresponds to a given adsorption amount and T is the temperature. The subscript w on the bracket on the RHS indicates that the bracketed quantity is to be evaluated at a constant coverage. This calculation assumes that the adsorbate can be considered an ideal gas and the volume of the adsorbed phase is negligible compared to the bulk gas phase.³³ These assumptions can be valid at low adsorptions. It can be identified that the bracketed quantity is the slope of an isostere (i.e. a plot of ln P vs. 1/T). For this, several isotherms at low pressures were simulated for the three fluids in S4Z4 at seven closely spaced temperatures. The pressure corresponding to the adsorption of one fluid molecule was obtained by interpolation of these isotherms. The corresponding ln P was then plotted as a function of 1/T (Fig. 4) and the Q_ads determined as in eqn (1).

Fig. 4 Adsorption isosteres corresponding to the adsorption of one fluid molecule in S4Z4. Circles represent the simulation data while the straight lines are linear fits to the data. The isosteric heat of adsorption for the three fluids calculated from the slope using eqn (1) is displayed.

The isosteric heat of adsorption of CO₂ (29.19 kJ mol⁻¹) in silicalite is about four times that of hydrogen (7.25 kJ mol⁻¹) and 1.5 times that of methane (19.28 kJ mol⁻¹). This strong difference in the adsorption heats is consistent with the high selectivity for CO₂ and methane for adsorption in silicalite from a mixture with hydrogen.

While the coadsorption isotherms displayed in Fig. 3 correspond to 298 K, in real UHS applications, higher and more variable temperatures can be envisaged. To explore the temperature dependence of coadsorption, we simulated the adsorption of the equimolar fluid mixtures in S4Z4 at a total pressure of 200 atm at several temperatures with a range that can be expected in UHS applications. Shown in Fig. 5 is the variation of adsorption amounts with temperature. The adsorption amounts of the carbon fluids follow the expected trend, i.e. adsorption amounts decrease with an increase in temperature. For hydrogen, the adsorption amounts seem to increase with temperature. However, this increase is barely beyond uncertainty because of the unusually small amounts of adsorption. At higher temperatures, as the carbon fluids occupy smaller number of adsorption sites, these sites become available for hydrogen. Nevertheless, the amount of hydrogen adsorbed remains extremely low.

Fig. 5 Amounts of fluids adsorbed in S4Z4 from an equimolar mixture with hydrogen at a total pressure of 200 atm, as a function of temperature.

3.3 Effects of mixture composition

Fig. 6 shows the effects of mixture composition on the adsorption amounts. As the hydrogen content in the feed mixture increases, the number of adsorbed hydrogen molecules increase while that of the carbon fluids decrease, as expected. However, while for methane–hydrogen mixture the number of the adsorbed molecules of the two species become comparable beyond a 95% molar composition of hydrogen in the mixture, for CO₂–H₂ mixture, CO₂ adsorption remains dominant even up to 98% hydrogen content in the feed. No significant qualitative difference can be seen in the variation of adsorption amounts as a function of hydrogen content between different adsorbents.

Fig. 6 Amounts of fluids adsorbed as a function of hydrogen content in the feed mixture for three adsorbent models. H₂–CO₂ and H₂–CH₄ respectively denote the data for hydrogen in a mixture with CO₂ and CH₄.

3.4 Effects of pore characteristics

Varying the pore space accessible for adsorption facilitates changing the pore characteristics like S and V of silicalite. However, the variation in S and V on selectively blocking the pores are correlated¹⁷ and therefore, it is advisable to study the effects of the combined variable S/V. Further, as the topologies of the straight and zigzag pores are different, a fair study on the effects of S/V would require keeping either all straight or all zigzag pores open. Thus, we assessed the effects of S/V for either S4Zm (m = 0, 1, 2, 3 or 4) adsorbents (labelled S-major because all straight channels are open in them) or SnZ4 (n = 0, 1, 2, 3 or 4) adsorbents (labelled Z-major). The effects of S/V in these two classes of adsorbents on the adsorption amounts of the three fluids in pure state as well as in the mixture at a pressure of 100 atm are shown in Fig. 7. For S-major adsorbents, a consistent decrease of the adsorption amounts with S/V can be seen for all fluids. For hydrogen in the mixture the variation is apparently absent because of the very small number of molecules adsorbed. For Z-major adsorbents, while the variation is not as systematic and smooth as for S-major adsorbents, in general higher S/V can still be seen to result in a lower adsorption amount for all fluids. The variation of adsorption amounts with S/V at a lower pressure of 1 atm shows a similar trend as at 100 atm as shown in the SI.

Fig. 7 Variation of adsorption amounts of hydrogen (left) and carbon fluids (right) at 100 atm with S/V of the adsorbents.

The adsorbent models studied here exhibit a wide range of pore connectivities ranging from 0 to 48. These different adsorbents however, have different pore volumes and therefore, for a fair comparison of the adsorption amounts between these adsorbents, the adsorption amounts are normalized to the available pore volume. These normalized adsorption amounts for all fluids in the pure as well as mixed state at a pressure of 100 atm are shown in Fig. 8 (see the SI for the corresponding data at 1 atm; Fig. S5). The variation of adsorption amounts of hydrogen with number of pore connections is within uncertainty. However, for the carbon fluids, an unambiguous decrease of adsorption amounts with the number of pore connections can be seen. This is consistent with the previous findings at higher temperature.¹⁴ When there is no connection between the pores (S4Z0 and S0Z4), the adsorption is noticeably higher.

Fig. 8 Variation of adsorption amounts of hydrogen (left) and carbon fluids (right) at 100 atm with the number of pore connections in the adsorbents.

Another parameter that characterizes the pore network is tortuosity. While tortuosity can be parametrized in terms of the diffusion coefficient of a fluid confined in the pore,³⁴ here we use the fraction of available pores that are tortuous as the relevant parameter, with S4Z0 and S0Z4 representing the extremes. As observed earlier, these extreme cases exhibit noticeably larger amounts of adsorption. Apart from these extremes, the adsorption amounts exhibit no systematic variation with tortuosity (see SI).

3.5 Structure and distribution of the confined fluids

Distribution functions for the pair of atoms belonging to the guest molecules and the adsorbent S4Z4 calculated from the simulations at 100 atm are shown in Fig. 9. The left panel shows the pair distribution of atoms of the guest molecules with the oxygen belonging to the adsorbent. As the silicon atoms in the adsorbent lie behind the oxygen atoms, the corresponding distributions (not shown) are shifted to larger distances with no other noticeable feature. Amongst all, oxygen atoms belonging to the CO₂ molecules are closest to the adsorbent oxygen which form the pore surface; followed by hydrogen. The first peak for the carbon belonging to CO₂ and methane are coincident. It is important to note here that both methane and hydrogen are represented in a united atom model and thus, the corresponding distribution represents the center of mass of the molecule instead of the individual atoms at the extremeties. The center of mass in case of CO₂ which is represented in fully atomistic model is coincident with the carbon atom. Thus, comparing only the molecular center of mass, it can be seen that hydrogen is located closest to the adsorbent pore surface. This is partly because a smaller kinetic diameter (d_k(H₂) = 2.89 Å;³⁵ d_k(CO₂) = 3.30 Å;³⁶ d_k(CH₄) = 3.80 Å)³⁷ lets hydrogen approach closer to the pore surface.

Fig. 9 Pair distribution functions calculated from the simulation of fluid adsorption in S4Z4 at 100 atm. The left panel shows the guest-adsorbent distributions with the adsorbent represented by the oxygen atom of silicalite. The middle panel shows the self atom distributions. The right panel shows the distribution of hydrogen with the carbon fluids from mixture adsorption simulation.

The like atom distributions (e.g. C–C) for all the guest atoms are shown in the middle panel of Fig. 9. The strong delta peak at 2.32 Å in the distribution of oxygen atoms belonging to the CO₂ molecules is an intramolecular peak resulting from the fixed distance of the two oxygen atoms belonging to the same molecule. Here too, two hydrogen molecules can be seen to come closer to each other than two CH₄ or CO₂ molecules.

The distributions shown in the left and middle panels are for pure gas adsorption, however, the correspondent data for the mixture adsorption is coincident and is therefore not shown. In case of mixture adsorption, it is instructive to investigate the distribution of atoms of the carbon fluid with hydrogen. These distributions are shown in the right panel of Fig. 9. For reference, the self (H–H) distribution is shown from the simulation of pure hydrogen as the corresponding distribution in the mixture simulation has poor statistics because of a smaller number of hydrogen molecules adsorbed in the mixture. Disregarding the oxygen atom, and considering only the center of mass, the positioning of the first peaks follows the order H₂ < CO₂ < CH₄. This is consistent with the kinetic diameters of the three molecules (d_k(H₂) < (d_k(H₂) + d_k(CO₂))/2 < (d_k(H₂) + d_k(CH₄))/2).

While pair distribution functions provide a one-dimensional picture of the interactions between atoms, a two dimensional map of atomic distribution in the simulation cell can provide additional insights. In Fig. 10, we show the distribution of the center of mass of the three fluid molecules in representative straight (left vertical panels) and zigzag (right horizontal panels) channels of S4Z4 from the corresponding simulations of pure gas adsorption at 100 atm. The intensity in these maps represent the probability of finding the center of mass of the fluid molecule at a given location. White regions represent zero probability with the zeolite atoms occupying the corresponding locations. The instensity varies logarithmically between low probability represented by bluish hues and higher probabilities represented by yellower hues. Several differences between the distributions of the molecules can be noted. Firstly, the distribution of hydrogen is the most homogenous, with relatively less variation in the intensity, followed by CH₄. CO₂ on the other hand shows strongest heterogeneity with several high intensity regions separated by zero intensity. These disparate regions represent points of strong adsorption. In particular, while hydrogen and methane occupy the intersections almost homogenously, for CO₂ distinct regions can be seen separated by a region of almost zero probability. This is equivalent to adsorption layers on opposite pore surfaces separated from each other with a bulk like low intensity region. As shown in a previous publication, this has conseuqences for the S/V dependence of the diffusion coefficient.¹⁷

Fig. 10 Probability distribution maps of the three fluid molecules in the straight (left vertical panels) and zigzag (right horizontal panels) channels of S4Z4 at 100 atm. The intensity represents the probability of finding a molecule at a given location and varies logarithmically between lower probabilities represented by bluer hues and higher probabilities represented by yellower hues. White regions represent zero probability.

Another notable feature is the width of the finite probability region, which is largest for hydrogen. This is because hydrogen can approach the pore surface to closer distances as demonstrated in Fig. 9 and therefore has an effectively wider pore volume available for adsorption.

Finally, the heterogeneity in the distribution of all fluid molecules shown in Fig. 10 is seen to be smoother at lower pressure of 1 atm (see SI, Fig. S8 for the corresponding plot of the probability distribution). A similar reduction in the heterogeneity at lower pressure was earlier observed for both CO₂ and ethane at 308 K and was attributed to the dominance of fluid–fluid interactions over the fluid–host interactions at higher pressures.¹⁶ This suggests a coverage dependent enhancement in adsorption at higher pressures as seen for n-hexane on a silica surface.³⁸

4. Discussion

As noted in the introduction, although hydrogen has a high energy density, storing it for the purpose of energy usage is challenging because of its thermophysical properties. Being a light molecule, at room temperature it has a kinetic energy which dominates over the adsorption energies at most surfaces. This means that at room temperature, it is difficult to adsorb it on a surface while methane and CO₂ are more readily adsorbed. This is also the reason for the relatively more homogenous distribution of hydrogen in the adsorbent compared to the carbon fluids. Because of a smaller kinetic diameter, hydrogen has relatively larger number of adsorption sites available. However, because of a higher kinetic energy, not all of these sites are occupied, resulting in an inefficient adsorption and overall lower adsorption capacities of hydrogen. This is also reflected in the low isosteric heat of adsorption of hydrogen in silicalite compared to methane and CO₂. Further, in presence of the two carbon fluids, the relative uptake of hydrogen is reduced even further. This is in spite of a significantly smaller kinetic diameter of hydrogen compared to the carbon fluids. The effects of CO₂ on the uptake of hydrogen are significantly stronger than that of methane because of (i) a smaller kinetic diameter of CO₂ compared to methane; and (ii) a stronger interaction of CO₂ with silicalite atoms compared to methane, which in turn is probably because of the quadrupole moment of CO₂. This combined effect is also seen in the smaller uptake of methane compared to CO₂ in the pure gas adsorption and is reflected in the isosteric heat of adsorption of methane in silicalite being 1.5 times lower than CO₂.

For UHS applications, the host geological formations can have a diverse mineralogical and chemical composition. SiO₂ which forms silicalite studied here is expected to be a major component of these formations. In addition, clay mineral and organic content can be significant components. The isosteric heats of adsorption calculated for the three fluids in this work can be expected to decrease with higher adsorption amounts and thus provides an upper bound estimate for a single molecule adsorption. For hydrogen adsorption in clay-rich geological media, Wang et al.³⁹ report isosteric heat of adsorption values ranging between 7.4 and 25.1 kJ mol⁻¹, which are higher than that obtained here for hydrogen in silicalite. Ramirez-Vidal et al.⁴⁰ also report values ranging between 8 and 10.1 kJ mol⁻¹ for adsorption of hydrogen in activated carbon materials at 253 K. For hydrogen in low-maturity shales, Zheng et al.⁴¹ report adsorption heats in the range 11.81–12.29 kJ mol⁻¹. Alanazi et al.⁴² also report similar range (6.43–11.8 kJ mol⁻¹) for hydrogen in organic-rich shales. Thus, hydrogen adsorption in SiO₂ seems to be marginally weaker than other minerals expected to constitute UHS formations. Also notable are the adsorption energies of methane in shales reported by Li et al.⁴³ to fall in the range 17–25 kJ mol⁻¹, which is similar to the adsorption energy of methane in silicalite obtained here. Further, Chang et al.⁴⁴ report very high adsorption energies for CO₂ in shale (30–51 kJ mol⁻¹), roughly four times that of methane (9–13 kJ mol⁻¹) in the same material. The comparative analysis of the adsorption energies mentioned above suggests that the order of adsorption selectivity (CO₂ > CH₄ > H₂) is similar in different adsorbents and should thus be applicable to most UHS formations.

The temperature variation of adsorption amounts of hydrogen in the mixture with the carbon fluids in unmodified silica at 200 atm total pressure suggests that at temperatures higher than the room temperature, the adsorption amounts of the carbon fluids decreases, whereas that for hydrogen remains very low (Fig. 5). In pure adsorption, while methane and CO₂ seem to have reached saturation loadings at 100 atm, hydrogen adsorption remains unsaturated (Fig. 3). This suggests that at total pressures higher than 100 atm, hydrogen might start to populate the pores even as the carbon fluids are saturated. Alternatively, higher adsorption amounts can be obtained for hydrogen at very low temperatures which are practically difficult to achieve for UHS applications.

Selectively blocking some pores of silicalite in this study helped constrain the effects of pore characteristics on the adsorption of the three fluids and their mixtures. However, due to the lower amounts of hydrogen adsorption, the effect of pore connectivity is not strong enough beyond uncertainty, while for the two carbon fluids, the connectivity of pores can be seen to unambiguously result in a suppression of adsorption amounts. A similar effect was also seen for ethane and CO₂ at 308 K and was explained in terms of the blocking matter providing additional adsorption sites in low connectivity adsorbents.¹⁴

While the effect of pore connectivity on the adsorption of hydrogen is not clear, it is remarkable that the S/V of the adsorbents has a clear effect on the adsorption amounts of all the fluids in spite of the much inefficient and sparse adsorption of hydrogen. A larger S/V is in general seen to suppress the adsorption. A higher S/V suggests a combined effect of a larger S and a smaller V. S of an adsorbent can be expected to favor adsorption by providing a larger number of adsorption sites. In mesoporous systems, an increase in V will favor adsorption only at higher pressures when the adsorption is dominated by the available volume. At lower pressures, up to a complete surface coverage, the volume of the mesopore might not affect the adsorption. For narrower pores with diameters comparable to the kinetic diameter of the adsorbate however, the pore volume becomes important even at lower pressures as it allows access to the adsorption sites. Thus, adsorption as a function of S/V is a trade-off between higher S that favors adsorption and a lower V that suppresses adsorption. In the present study, it appears that the effect of V is stronger than that of S to result in a lowering of adsorption amounts at higher S/V. This is especially true for inefficiently adsorbed hydrogen, as an increase in S may result in an increase in adsorption sites that remain largely unoccupied. On the other hand, with a smaller kinetic diameter, and closer approach to the pore surface as seen in Fig. 9 and 10, hydrogen has an effectively wider volume of the pore available and is therefore relatively more sensitive to small changes in V. For larger pores where the changes in V do not impact the adsorption except at very high pressures, the relation might be different. Amongst the porous media relevant for UHS applications, only shale formations exhibit a preponderance of pore diameters below 2 nm.¹² For sandstones and other media, the pores are wider and hence the S/V lower. The S/V values investigated in this work are closer to the upper limit for UHS applications.

This study focused on co-adsorption, providing useful information on how the pore characteristics and potential cushion gases like methane and CO₂ might impact the adsorption of hydrogen. Also, it has been found that for hydrogen, the pore shape and connectivity might not be important determinants of adsorption. This is because the room temperature adsorption is inefficient with most of the adsorption sites remaining unoccupied. For UHS applications, a study that involves pushing the carbon fluids in the pores that already have hydrogen can provide further important insights. Such a study should account for transport properties of the fluids including dynamic displacement and breakthrough curves and might require multiple tools. The results obtained here might inform such a study in the future that could provide insights on reservoir selection.

5. Conclusions

We carried out a grand canonical Monte Carlo simulation study on the room temperature adsorption of hydrogen and its co-adsorption with two potential cushion gases – methane and CO₂ – in models of silicalite with varying degrees of pore characteristics like connectivity, tortuosity and S/V ratios. At room temperature, the adsorption of hydrogen is inefficient leading to low amounts of adsorption in all adsorbent models. In unmodified silicalite the adsorption amounts of hydrogen, methane and CO₂ in pure state at 100 atm are respectively (1.08 ± 0.11) mmol g⁻¹; (2.86 ± 0.06) mmol g⁻¹ and (3.52 ± 0.03) mmol g⁻¹. This difference in the adsorption amounts is reflected in the isosteric heat of adsorption, which for CO₂ is four times stronger compared to hydrogen and 1.5 times higher compared to methane. In the presence of methane and CO₂ the adsorption amounts of hydrogen are reduced further, with CO₂ presence exhibiting a stronger reduction due to its stronger interaction with the pore surface and smaller kinetic diameter. Hydrogen adsorption is found to be relatively more homogenous with different adsorption sites exhibiting similar adsorption strength. This is in contrast to CO₂ that adsorbs more heterogeneously with strong adsorption sites separated in space. The S/V of the adsorbents is found to unambiguously suppress the adsorption amounts of all fluids. The information obtained here can be useful in designing a future study that recreates the UHS scenario where the cushion gases are pushed into a reservoir that already has hydrogen gas.

Author contributions

Conceptualization, methodology, validation, formal analysis, investigation, visualization, project administration and data curation, S. G.; writing—original draft preparation, S. G.; writing—review and editing, S. G. and D. C.; funding acquisition, D. C. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Data availability

All data related to this study are available in this article and the associated supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5ra06520h.

Acknowledgements

This research was funded by the State of Ohio through the Third Frontier Ohio Research Scholar Program. We would like to acknowledge STFC's Daresbury Laboratory for providing the package DL-Monte,²⁹ which was used to obtain the quantities used in this work. Figures in this manuscript were made using the freely available plotting software Gnuplot.⁴⁵ We would also like to thank the reviewers for their constructive comments.

References

1. J. Symons, in Energy Security in the Era of Climate Change: The Asia-Pacific Experience, ed. L. Anceschi and J. Symons, Palgrave Macmillan, London UK, 2012, pp. 1–10.
2. W. Ma and W. Wang, Heliyon, 2024, 10(8), e29712.
3. C. Greig, Engineering, 2017, 3(4), 436–438.
4. L. Schlapbach and A. Züttel, Nature, 2001, 414(6861), 353–358.
5. M. Yue, H. Lambert, E. Pahon, R. Roche, S. Jemei and D. Hissel, Renewable Sustainable Energy Rev., 2021, 146, 111180.
6. M. D. Allendorf, V. Stavila, J. L. Snider, M. Witman, M. E. Bowden, K. Brooks, B. L. Tran and T. Autrey, Nat. Chem., 2022, 14(11), 1214–1223.
7. R. Tarkowski, Renewable Sustainable Energy Rev., 2019, 105, 86–94.
8. G. O. Taiwo, O. S. Tomomewo and B. A. Oni, J. Energy Storage, 2024, 90, 111844.
9. M. Zamehrian and B. Sedaee, J. Petrol. Sci. Eng., 2022, 212, 110304.
10. M. Kanaani, B. Sedaee and M. Asadian-Pakfar, J. Energy Storage, 2022, 45, 103783.
11. X. Qiu, H. Liu, M. Liu, H. Mao, D. Wang, Q. Ying and S. Ban, Energies, 2023, 16(5), 2096.
12. F. N. Yang and B. P. Wang, Acta Pet. Sin., 2013, 34(2), 301.
13. M. M. Kummali, D. Cole and S. Gautam, Membranes, 2021, 11(2), 113.
14. S. Gautam and D. R. Cole, ChemEngineering, 2021, 5(3), 55.
15. S. Gautam and D. R. Cole, Phys. Chem. Chem. Phys., 2022, 24(19), 11836–11847.
16. S. Gautam and D. R. Cole, C, 2023, 9(4), 116.
17. S. Gautam and D. R. Cole, C, 2025, 11(3), 57.
18. A. Raza, M. Mahmoud, S. Alafnan, M. Arif and G. Glatz, Int. J. Mol. Sci., 2022, 23(21), 12767.
19. A. Alanazi, H. R. Abid, M. Usman, M. Ali, A. Keshavarz, V. Vahrenkamp, S. Iglauer and H. Hoteit, Fuel, 2023, 346, 128362.
20. I. Díaz, E. Kokkoli, O. Terasaki and M. Tsapatsis, Chem. Mater., 2004, 16(25), 5226–5232.
21. H. Van Koningsveld, H. Van Bekkum and J. C. Jansen, Struct. Sci., 1987, 43(2), 127–132.
22. K. Momma and F. Izumi, Appl. Crystallogr., 2011, 44(6), 1272–1276.
23. M. G. Martin and J. I. Siepmann, J. Phys. Chem. B, 1998, 102(14), 2569–2577.
24. J. J. Potoff and J. I. Siepmann, AIChE J., 2001, 47(7), 1676–1682.
25. R. T. Cygan, J. J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108(4), 1255–1266.
26. V. Buch, J. Chem. Phys., 1994, 100(10), 7610–7629.
27. S. Gautam, D. R. Cole, Z. I. Dudás and I. Dhiman, ChemPhysChem, 2024, 25(17), e202400360.
28. F. Darkrim and D. Levesque, J. Chem. Phys., 1998, 109(12), 4981–4984.
29. J. A. Purton, J. C. Crabtree and S. C. Parker, Mol. Simul., 2013, 39(14–15), 1240–1252.
30. T. C. Golden and S. Sircar, J. Colloid Interface Sci., 1994, 162(1), 182–188.
31. M. S. Sun, D. B. Shah, H. H. Xu and O. Talu, J. Phys. Chem. B, 1998, 102(8), 1466–1473.
32. T. F. Willems, C. H. Rycroft, M. Kazi, J. C. Meza and M. Haranczyk, Microporous Mesoporous Mater., 2012, 149(1), 134–141.
33. S. Poyet and S. Charles, Cem. Concr. Res., 2009, 39(11), 1060–1067.
34. N. Epstein, Chem. Eng. Sci., 1989, 44(3), 777–779.
35. M. Kanezashi and T. Tsuru, in Membrane Science and Technology, ed. S. Ted Oyama and S. M. Stagg-Williams, Elsevier, Oxford, UK, 2011, vol. 14, pp. 117–136.
36. J. Yang, Q. Zhao, H. Xu, L. Li, J. Dong and J. Li, J. Chem. Eng. Data, 2012, 57(12), 3701–3709.
37. S. Inthawong, A. Wongkoblap, W. Intomya and C. Tangsathitkulchai, Molecules, 2023, 28(14), 5433.
38. J. I. Chevallier-Boutell, R. H. Acosta, M. B. Franzoni and J. A. Olmos-Asar, Microporous Mesoporous Mater., 2024, 375, 113174.
39. L. Wang, J. Cheng, Z. Jin, Q. Sun, R. Zou, Q. Meng, K. Liu, Y. Su and Q. Zhang, Fuel, 2023, 344, 127919.
40. P. Ramirez-Vidal, R. L. Canevesi, G. Sdanghi, S. Schaefer, G. Maranzana, A. Celzard and V. Fierro, ACS Appl. Mater. Interfaces, 2021, 13(10), 12562–12574.
41. R. Zheng, K. Hu, Y. Wu, X. Li and G. Li, Fuel, 2025, 401, 135805.
42. A. Alanazi, H. R. Abid, I. S. Abu-Mahfouz, S. A. Bawazeer, T. Matamba, A. Keshavarz, S. Iglauer and H. Hoteit, Fuel, 2025, 379, 132964.
43. J. Li, Z. Chen, K. Wu, K. Wang, J. Luo, D. Feng, S. Qu and X. Li, Chem. Eng. J., 2018, 349, 438–455.
44. X. Chang, S. Lin, C. Yang, Y. Guo and X. Dai, Phys. Fluids, 2024, 36(11), 116614.
45. T. Williams, C. Kelley, R. Lang, D. Kotz and J. Campbell, Gnuplot 4.6. 6: An Interactive Plotting Program, 2014. http://gnuplot.info, accessed on August 28, 2025.

---