Hydrogen storage and diffusion in shale nanopores: molecular dynamics effects of pressure, temperature, confinement, and wall chemistry

1 Introduction

Hydrogen is increasingly considered an energy carrier for balancing variable renewable electricity and decarbonizing sectors that are difficult to electrify directly (Wang et al., 2018; Evangelopoulou et al., 2019; He et al., 2021; Shiraishi et al., 2024). Its large-scale deployment requires storage capacities and cycling durations beyond those of conventional surface vessels (Tarkowski, 2019; Ali et al., 2025). Underground hydrogen storage (UHS) in salt caverns, depleted hydrocarbon reservoirs, and saline aquifers can provide seasonal-scale capacity while reducing surface land use (Tarkowski, 2019; Ali et al., 2025). Salt caverns provide the clearest engineering precedent, but suitable salt formations are geographically restricted (Caglayan et al., 2020; Zheng et al., 2020). Porous formations may allow partial reuse of wells, gathering systems, and subsurface characterization developed during hydrocarbon production (Okoroafor et al., 2022; Ali et al., 2025).

Hydrogen presents distinctive containment and recovery challenges in porous media (Gomez Mendez et al., 2024; Ren et al., 2025; Shi et al., 2025). Its low molecular mass, low viscosity, and high diffusivity favor rapid migration, whereas adsorption, dissolution, capillary trapping, residual gases, water, geochemical reactions, and microbial conversion can alter both stored inventory and deliverability (Gomez Mendez et al., 2024; Baru et al., 2025; Shi et al., 2025). A formation therefore cannot be assessed from geometric pore volume alone. The spatial distribution of H₂ near pore walls and its diffusion through the connected pore system must be considered together, particularly in fine-grained rocks where much of the internal surface area occurs in nanopores and continuum descriptions may not resolve interfacial structuring or strong confinement.

Depleted shale gas formations are potential UHS targets because shale combines low matrix permeability with abundant organic and mineral nanopores (Zhang M. et al., 2024; Chen et al., 2026; Liang et al., 2026). Kerogen provides chemically heterogeneous carbonaceous surfaces (Craddock et al., 2018; Duan et al., 2018; Zhang et al., 2022), whereas silicates, clays, carbonates, and carbon-rich components provide different charge distributions, surface corrugations, and adsorption environments (Al-Harbi et al., 2023; Zheng et al., 2024; Zhao et al., 2025). Previous experimental and molecular-simulation studies have shown that H₂ storage and transport are sensitive to temperature, pressure, pore aperture, kerogen maturity, residual gases, water, and mineral identity (Kahzadvand et al., 2024; Zhang M. et al., 2024; Wang and Bai, 2026). However, these factors are commonly investigated separately, making it difficult to compare their relative influence on H₂ storage and diffusion within a unified molecular framework. In particular, the coupled response of pore-scale H₂ distribution and diffusion to thermodynamic conditions, nanoconfinement, organic-matter type, and pore-wall chemistry remains insufficiently resolved. A systematic multi-factor comparison is therefore required to compare the observed controls on H₂ retention and transport in shale nanopores.

This study compares H₂ storage and diffusion in shale slit-pore models using a unified molecular simulation and analysis framework. Its specific contribution is the joint evaluation of adsorption heat, H₂ population, pore-normal density profiles, and self-diffusion coefficients across variations in temperature, pressure, pore width, kerogen type, and pore-wall type. Rather than treating any single descriptor as a proxy for all storage behavior, the analysis examines whether energetic affinity, pore occupancy, and molecular diffusion respond consistently to changes in the investigated factors. Three questions are addressed: (i) what storage and diffusion characteristics H₂ exhibits in the reference quartz pore; (ii) how each investigated factor changes the complementary descriptors; and (iii) whether the rankings based on adsorption heat, H₂ population, and diffusion coincide. This comparison identifies the largest observed responses within the tested ranges and clarifies the distinction between energetic binding, pore occupancy, and molecular diffusion. The resulting molecular-scale insights can support the evaluation of shale pore structures and the assessment of retention–transport trade-offs in underground hydrogen storage.

2 Molecular models and simulation methods

2.1 Molecular model construction

2.1.1 Mineral crystal and pore-wall models

Albite (NaAlSi₃O₈) (Harlow, 1982), illite [KAl₂(Si₃Al)O₁₀(OH)₂] (Zhang X. et al., 2024), quartz (SiO₂) (Liu et al., 2018), calcite (CaCO₃) (Graf, 1961), and graphite (Kwiecińska and Petersen, 2004; Hou et al., 2019) were selected as representative inorganic pore-wall components of shale. Their crystal structures were constructed from crystallographic data in Materials Studio and geometrically optimized before replication. The optimized unit-cell structures are shown in Figure 1.

Each unit cell was replicated to form a mineral supercell with lateral dimensions of at least 40 Å and a thickness of at least 15 Å (Figure 2). The exposed surfaces were albite (010), illite (001), calcite (104), graphite (001) and quartz (001), as summarized in Table 1. During the production simulations, the mineral atoms were fixed to maintain the specified pore geometry. Related molecular studies of confined fluids and mineral interfaces provide the context for this idealized pore-wall representation (Greathouse et al., 2015; Mu et al., 2022).

TABLE 1. Molecular compositions and pore-wall characteristics used in the simulations.

2.1.2 Kerogen molecular and pore-wall models

Three kerogen models were constructed following the type-dependent molecular framework reported by Ungerer et al. (2015). Type I kerogen is aliphatic-rich and has the highest H/C ratio; type II contains both aliphatic and aromatic moieties; and type III has the highest aromatic fraction and oxygen content. Their molecular formulas, elemental ratios, and structural descriptors are listed in Table 1. Each kerogen molecule was first geometrically optimized and then packed into a periodic amorphous cell. Five annealing cycles between 300 and 600 K were performed to reduce configurational bias, followed by energy minimization to obtain stable kerogen wall structures. The molecular structures and equilibrated amorphous cells of the three kerogen types are shown in Figure 3.

2.1.3 Slit-pore model assembly

For each inorganic model, two identical mineral slabs were placed parallel to the xy plane, with their inner surfaces facing each other along the z direction. For the organic models, two copies of the equilibrated kerogen cell were arranged in the same manner. The separation between the opposing inner surfaces was adjusted to 2, 3, 4, or 5 nm to generate slit pores of different widths. Periodic boundary conditions were applied in all three directions, and the pore normal was aligned with z.

The quartz-based slit pore at 323 K, 5 MPa, and a width of 5 nm was used as the common reference system. Quartz was selected as a representative silicate component of shale, and retaining this single pore-wall type in temperature, pressure, and pore-width series allowed those effects to be compared without simultaneously changing the wall composition (Liu et al., 2018). The three kerogen models and the five mineral models were used for the respective composition comparisons. The models and their exposed surfaces are summarized in Table 1.

2.2 Hydrogen loading and simulation matrix

Hydrogen was represented as a rigid diatomic molecule with an H–H bond length of 0.074 nm. For each slit-pore volume, the bulk H₂ density at the target temperature and pressure was calculated with the Peng–Robinson equation of state (Peng and Robinson, 1976). The corresponding number of H₂ molecules was prescribed from that density and the accessible pore volume and inserted before the NVT trajectory. The pressure values reported for these comparisons therefore specify the bulk reference states used to initialize the loading; pressure was not independently controlled during the NVT production runs. The total H₂ number, pore volume, and temperature were fixed in each run. Consequently, the reported H₂ population depends partly on the prescribed loading and is not an independently predicted equilibrium adsorption capacity. The simulation matrix comprised temperatures of 313, 323, 333, 343, and 353 K; pressures of 5, 15, 25, and 35 MPa; pore widths of 2, 3, 4, and 5 nm; three kerogen types; and five inorganic pore-wall types.

The reference condition was 323 K, 5 MPa, and a 5 nm quartz-based slit pore. These values were adopted as a common comparison point shared by the different parameter series, rather than as an optimized storage condition. Temperature was varied at 5 MPa and a pore width of 5 nm; pressure was varied at 323 K and a pore width of 5 nm; and pore width was varied at 323 K and 5 MPa. The kerogen-type and pore-wall-type comparisons were also performed at 323 K, 5 MPa, and a pore width of 5 nm. This arrangement permits each reported comparison to be interpreted against the same reference settings (Table 2).

TABLE 2. Simulation matrix and common numerical settings.

2.3 Force field and molecular dynamics protocol

All models were built and simulated in BIOVIA Materials Studio 2023 using the COMPASS III force field. The COMPASS III framework combines bonded terms with Lennard–Jones and electrostatic interactions. Unlike-atom Lennard–Jones parameters were obtained with Lorentz–Berthelot mixing rules. A 1.25 nm cutoff was used for nonbonded interactions, and long-range electrostatics were evaluated by Ewald summation with a convergence criterion of 1 × 10⁻⁶ kcal mol^-1. Three-dimensional periodic boundary conditions were applied, with the slit-pore normal along z. COMPASS III was used consistently for the organic, mineral, and H₂ components to provide a unified force-field framework for comparison among the simulated systems (Zhang X. et al., 2024b; Dang et al., 2025).

Each H₂-filled slit-pore system was first minimized using the conjugate-gradient method, with a maximum-force criterion of 0.001 kcal mol^-1 Å^-1 and a maximum of 10,000 iterations. The minimized system was then equilibrated for 500 ps in the NVT ensemble at the target temperature using a Nosé–Hoover thermostat and a 1 fs time step. Equilibration was confirmed by the stabilization of temperature and total energy without systematic drift. A subsequent 100 ps NVT production simulation was performed under the same thermodynamic conditions (Jia et al., 2026). Only the stabilized 100 ps production trajectory was used to determine adsorption heat, H₂ population, density profiles, MSD, and diffusion coefficients. Three independently initialized simulations were conducted for each condition, and the mean values were used for comparison.

2.4 Data analysis

2.4.1 Adsorption heat, H₂ population, and density profiles

The adsorption heat, Q a d s, was evaluated using the interaction-energy difference for the H₂–pore-wall system. For each configuration sampled from the stabilized production trajectory, the mean binding energy per H₂ molecule was calculated using Equation 1:where E t o t a l is the total potential energy of the combined H₂–pore-wall system, E H 2 and E w a l l are the potential energies of the H₂ and wall subsystems used as the energy references, respectively, and N H 2 is the number of H₂ molecules used to normalize the energy difference for that configuration. The resulting Q a d s values were averaged over the stabilized production trajectory and are reported in kcal mol^-1 on a per-H₂ basis. An attractive interaction lowers E t o t a l relative to E H 2 + E w a l l and therefore gives a positive Q a d s; a larger positive value denotes stronger mean energetic binding under this sign convention. Here, “adsorption heat” denotes this interaction-energy-based descriptor, rather than an independently determined thermodynamic isosteric heat. Related energy-difference descriptions of fluid–solid adsorption have been used for shale nanopores (Liu et al., 2024; Zhao et al., 2025).

Q a d s = E H 2 + E w a l l − E t o t a l N H 2 (1)

H₂ population was defined as the mean number of H₂ molecules present within the slit-pore region during the stabilized production stage and is reported in molecules. It describes occupancy under the prescribed loading, rather than an adsorbed-only population or an independently predicted equilibrium adsorption capacity. The total H₂ number in each NVT trajectory remains fixed at its initialized value; the trajectory samples the distribution and motion of those molecules, rather than exchange with an external gas reservoir. Population comparisons among conditions therefore reflect the loading procedure in Section 2.2 together with the specified pore region. They should not be read as equilibrium adsorption isotherms or excess adsorption capacities. Pore-normal H₂ mass-density profiles were obtained by dividing the z direction into 0.08 nm slabs and averaging the molecular distribution over 1,000 production-trajectory records. The distance coordinate was referenced to the pore center. Together, Q a d s, H₂ population, and density profiles characterize mean energetic binding, pore occupancy, and the spatial distribution of H₂, respectively. Related pore-scale analyses have also used joint energetic and spatial characterization for shale fluids (Zhang X. et al., 2026b).

2.4.2 Mean squared displacement and diffusion coefficient

Mean squared displacement (MSD) was calculated from the center-of-mass trajectories of the H₂ molecules and is reported in nm². The total MSD contains the x, y, and z displacement components. Its definition and the scalar Einstein coefficient used for reporting the results are given in Equations 2, 3 (Ernst and Köhler, 2012; Nagai et al., 2020):

M S D t = 1 N ∑ i = 1 N r i t − r i 0 2 (2)

D = 1 6 lim t → ∞ d M S D t d t (3)

Here N is the total number of H₂ molecules, rᵢ(t) is the center-of-mass position vector of molecule i, and t is elapsed time in the production trajectory. The squared norm includes all three Cartesian components, so M S D x y z = M S D x + M S D y + M S D z. The linear portion of the total MSD from the stabilized 100 ps production trajectory was fitted with slope k in nm² ps^-1. The reported coefficient is D = k/6 in nm² ps^-1, or 100k/6 in Ų ps^-1 because 1 nm² = 100 Ų; equivalently, one Ų ps^-1 = 10⁻⁸ m² s^-1. The same analysis and unit convention were applied to every condition.

The slit pore is anisotropic, and the factor 1/6 is used to define a directionally averaged scalar, not to assume equal diffusion in all directions. If the directional slopes are evaluated over the same fitting interval, D = D x + D y + D z / 3, with each directional coefficient defined by one-half of the corresponding MSD slope. This trace-based average provides a consistent summary of the total molecular displacement for the present comparisons. In contrast, the in-plane coefficient is D ∥ = 1 4 d M S D x + M S D y dt. The pore-normal MSD is bounded at sufficiently long times by the opposing walls, so the reported D must not be identified with D ∥ or with an isotropic bulk diffusivity. Anisotropic H₂ diffusion in clay slits has also been reported by Liu et al. (2022). The values retained here are derived from the total three-dimensional MSD; no separately evaluated in-plane values are reported.

3 Results and discussion

3.1 H₂ storage and diffusion characteristics in a shale slit pore

At 323 K, 5 MPa, and a 5 nm quartz-based slit pore, H₂ occupies the full accessible pore (Figure 4a). The density profile is approximately symmetric, with a central value of about 0.0033 g cm^-3 and modest near-wall oscillations (Figure 4b). The adsorption heat is 0.67 kcal mol^-1 and the H₂ population is 94.96 molecules. These energetic, quantitative, and spatial quantities characterize the H₂ distribution under the reference loading; the population is not an independently calculated adsorption capacity.

The MSD increases approximately linearly over the stabilized 100 ps production trajectory (Figure 4c). Linear fitting gives the reported directionally averaged self-diffusion coefficient of 409.32 Ų ps^-1, equivalent to 4.0932 nm² ps^-1 (or 4.09 × 10⁻⁶ m² s^-1). Through the Einstein relation, this coefficient corresponds to an MSD slope of 24.5592 nm² ps^-1, making the conversion between the plotted MSD unit and the reported diffusion unit explicit. Together, the adsorption heat, H₂ population, density distribution, and diffusion coefficient provide a complementary basis for the subsequent factor comparisons.

The storage and diffusion characteristics obtained here are broadly consistent with previous molecular studies of H₂ in geological nanopores. Yang et al. reported weak physical adsorption of H₂ in quartz slit pores, with relatively low adsorption heats and no pronounced multilayer accumulation, consistent with the modest near-wall density modulation observed in the present quartz system (Yang et al., 2023). Previous molecular dynamics simulations further showed that H₂ self-diffusion in quartz and other geological slit pores is sensitive to pressure, pore confinement, and surface chemistry, with diffusivity generally decreasing with increasing pressure and increasing with pore aperture (Zheng et al., 2024). Similar transport trends have also been reported for kerogen and clay nanopores (Raza et al., 2022a; Oliver et al., 2024). These agreements in interfacial distribution and transport behavior support the physical consistency of the present molecular model.

3.2 Effects of controlling factors on H₂ storage and diffusion

The five factors were evaluated using four complementary descriptors: adsorption heat, H₂ population, pore-normal density distribution, and the MSD-derived, directionally averaged self-diffusion coefficient. The adsorption heat represents mean H₂–wall energetic binding, H₂ population and density profiles describe pore occupancy and its spatial distribution, and D quantifies molecular motion. Throughout the manuscript, “pore-wall type” identifies the investigated material variable, whereas “pore-wall chemistry” describes the interfacial properties used to interpret its effects. Each series varied one factor while keeping the remaining parameters fixed at their reference values: 323 K, 5 MPa, and a pore width of 5 nm, as applicable.

3.2.1 Temperature

At 5 MPa and a 5 nm quartz-based pore, adsorption heat decreases monotonically from 0.69 kcal mol^-1 at 313 K to 0.67, 0.65, 0.63, and 0.61 kcal mol^-1 at 323, 333, 343, and 353 K, respectively (Figure 5a). Over the same temperature range, the H₂ population decreases from 100.47 to 94.96, 94.22, 90.69, and 89.31 molecules (Figure 5b). The density profiles retain a broad pore-center distribution and exhibit temperature-dependent variations in local density and interfacial modulation (Figure 5c). The lower population is associated in part with the temperature-dependent initialization density. Together with the energy change, the results describe temperature-dependent occupancy and H₂–wall interactions under this loading protocol.

The self-diffusion coefficient exhibits the opposite trend, increasing from 352.88 Ų ps^-1 at 313 K to 409.32, 450.14, 481.96, and 500.30 Ų ps^-1 at 323, 333, 343, and 353 K, respectively (Figure 5d). The overall increase is 41.8%. Higher molecular kinetic energy accelerates translational motion and facilitates H₂ redistribution away from the pore walls, thereby increasing diffusion. The positive temperature dependence agrees with previous H₂ diffusion simulations in organic and mineral nanopores (Owusu et al., 2022; Yang and Jin, 2025; Zhao et al., 2025).

The coupled changes in adsorption heat, H₂ population, and D reveal a clear storage–deliverability trade-off. Lower temperature strengthens adsorption and favors H₂ retention, whereas higher temperature promotes molecular diffusion and facilitates injection–withdrawal redistribution. Temperature should therefore be considered together with fracture connectivity, water saturation, pressure gradients, and cycling frequency when selecting storage intervals and operating conditions.

3.2.2 Pressure

At 323 K and a 5 nm quartz-based pore, adsorption heat decreases from 0.67 kcal mol^-1 at 5 MPa to approximately 0.63 kcal mol^-1 at 15, 25, and 35 MPa (Figure 6a), indicating only a modest pressure dependence of the average H₂–wall interaction strength. Under the prescribed pressure-dependent loadings, the H₂ population rises from 94.96 molecules at 5 MPa to 285.28, 459.72, and 653.57 molecules at 15, 25, and 35 MPa, respectively (Figure 6b). The density profiles increase correspondingly throughout the pore across the prescribed loading states (Figure 6c). Because the initial molecular number increases according to the EOS-based protocol, this population increase is partly prescribed and must not be interpreted as an independently simulated increase in equilibrium adsorption capacity.

The self-diffusion coefficient decreases from 409.32 Ų ps^-1 at 5 MPa to 295.17 Ų ps^-1 at 15 MPa, 171.23 Ų ps^-1 at 25 MPa, and 136.05 Ų ps^-1 at 35 MPa (Figure 6d), giving an overall reduction of 66.8%. Higher H₂ density reduces intermolecular separation and increases collision frequency, thereby limiting the net displacement accumulated over a fixed time interval. The nearly stable adsorption heat above 15 MPa further indicates that pressure primarily regulates H₂ storage through pore filling and molecular crowding rather than through a substantial increase in interaction energy.

At the nanopore scale, the higher-pressure loading states contain more H₂ molecules and yield lower directionally averaged diffusion coefficients. The decrease in diffusivity is consistent with previous H₂ simulations in organic nanopores (Raza et al., 2022a), while the larger population must be interpreted in the context of the prescribed loading. These results describe an occupancy–diffusion trade-off for the simulated states; operating pressure in a storage project must additionally account for geomechanical constraints, caprock and well integrity, and withdrawal performance (Okoroafor et al., 2022).

3.2.3 Pore width

At 323 K and 5 MPa, adsorption heat increases from 0.62 kcal mol^-1 in the 2 nm pore to 0.63 kcal mol^-1 in the 3 and 4 nm pores and 0.67 kcal mol^-1 in the 5 nm pore (Figure 7a). Under the loading procedure in Section 2.2, the H₂ population increases with accessible pore volume, from 41.97 molecules at 2 nm to 59.51, 77.53, and 94.96 molecules at 3, 4, and 5 nm, respectively (Figure 7b). The density profiles broaden progressively with pore width, whereas the 2 nm pore contains no extensive pore-center region that is independent of both walls (Figure 7c). Thus, pore width controls storage through both accessible volume and the spatial extent of wall confinement.

The self-diffusion coefficient increases from 314.25 Ų ps^-1 at 2 nm to 389.61 Ų ps^-1 at 3 nm, 391.62 Ų ps^-1 at 4 nm, and 409.32 Ų ps^-1 at 5 nm (Figure 7d). The largest discrete increase occurs between 2 and 3 nm; above 3 nm, the diffusion coefficients change only slightly. This nonlinear response shows that diffusion becomes strongly restricted when the pore width approaches the combined influence range of the opposing pore walls.

Within the investigated quartz models, the largest discrete diffusion change occurs between 2 and 3 nm, consistent with stronger confinement in the narrower pore. This trend is supported by earlier molecular studies: Liu et al. (2022) reported increasing H₂ diffusivity with slit aperture in montmorillonite, while Oliver et al. (2024) found a similar pore-size dependence in graphene and kaolinite slit pores. Zheng et al. (2024) further showed that confinement depends on mineral–H₂ interactions in calcite, hematite, and quartz pores. These studies support the direction of the pore-size effect rather than a universal transition at exactly 2–3 nm. In the present models, narrower pores contain fewer H₂ molecules because of their smaller accessible volume and also restrict molecular motion. The observed contrast therefore illustrates the trade-off between available storage space and diffusion resistance; its relevance to connected shale pore networks requires consideration of the pore-size distribution and connectivity.

3.2.4 Kerogen type

At 323 K, 5 MPa, and a pore width of 5 nm, the adsorption heat is 0.63 kcal mol^-1 for type I and type II kerogen and increases slightly to 0.64 kcal mol^-1 for type III (Figure 8a). The corresponding mean H₂ populations are 53.80, 54.75, and 57.24 molecules (Figure 8b). The density profiles have similar overall shapes and show modest type-dependent differences (Figure 8c). In combination, the slightly higher mean adsorption heat and population in the type III model suggest a modest tendency toward stronger H₂ retention under the tested conditions.

The mean self-diffusion coefficients decrease from 417.75 Ų ps^-1 for type I kerogen to 407.23 Ų ps^-1 for type II and 398.74 Ų ps^-1 for type III (Figure 8d), a type I-to-type III difference of 4.6%. This spread is comparatively small relative to the pressure, pore-width, and pore-wall-type contrasts. The results are therefore described as a modest tendency toward enhanced retention in the type III model, rather than as evidence that this kerogen type is universally the most favorable for hydrogen storage.

This tendency may be associated with differences in aromatic structure, surface heterogeneity, and heteroatom content among the selected kerogen models, which can modify H₂–kerogen interactions and resistance to molecular translation. Previous studies have similarly related H₂ adsorption and transport to kerogen maturity, carbon content, microporosity, surface chemical heterogeneity, and residual methane (Raza et al., 2022b; Raza et al., 2023; Liang et al., 2026). These studies provide a mechanistic context, while the present comparison concerns kerogen type rather than a controlled maturity series. Kerogen type produced comparatively modest changes in the reported mean descriptors within the investigated models.

Organic–mineral heterogeneity provides an additional context for this comparison. Zhang L. et al. (2026) examined CO₂/CH₄ competition in type I–III kerogen–illite composite nanopores and found kerogen-dependent adsorption selectivity and different diffusion behavior for the two gases. This related result motivates consideration of composite organic–mineral interfaces when extending the present models. It concerns CO₂ and CH₄ rather than H₂, however, and is used here as mechanistic context rather than quantitative validation of the present H₂ ranking.

3.2.5 Pore-wall type

At 323 K, 5 MPa, and a 5 nm pore, the adsorption heat follows the order quartz (0.67 kcal mol^-1) > graphite (0.64 kcal mol^-1) > illite, albite, and calcite (approximately 0.63 kcal mol^-1) (Figure 9a). The H₂ population follows a different order: calcite (100.39 molecules) > albite (96.90) > quartz (94.96) > illite (87.33) > graphite (83.26) (Figure 9b). The pore-normal density profiles also display wall-dependent differences in local density and edge position (Figure 9c). These results show that adsorption heat and total H₂ population reflect complementary aspects of storage and do not necessarily follow the same mineralogical ranking.

The self-diffusion coefficients follow the order graphite (1,240.24 Ų ps^-1) > albite (637.10 Ų ps^-1) > illite (518.59 Ų ps^-1) > quartz (409.32 Ų ps^-1) > calcite (371.97 Ų ps^-1) (Figure 9d). H₂ therefore diffuses 3.3 times faster in the graphite pore than in the calcite pore. Calcite combines the highest H₂ population with the lowest D, whereas graphite combines the lowest population with the highest D. Quartz exhibits the highest adsorption heat but an intermediate population and diffusion coefficient. These population rankings refer to the individual loaded models and should not be taken as intrinsic mineral-specific adsorption capacities without a common volume or mass normalization.

The distinct rankings arise because wall-controlled storage and diffusion depend jointly on interaction energy, accessible adsorption sites, surface polarity, charge distribution, and structural corrugation. Previous molecular simulations have shown that mineral surface composition, interaction strength, and roughness can substantially modify H₂ adsorption and self-diffusion in geological slit pores. In particular, relatively smooth carbonaceous surfaces tend to permit faster molecular translation, whereas ionic and structurally heterogeneous mineral surfaces can impose stronger interfacial resistance to H₂ transport. The large wall-dependent contrast observed here therefore highlights mineralogical composition as an important control on H₂ retention and diffusion in shale nanopores (Al-Harbi et al., 2023; Zheng et al., 2024; Zhao et al., 2025).

3.3 Relative importance, engineering implications, and future perspectives

The integrated comparison shows that the largest response depends on the quantity considered (Table 3). Across the prescribed pressure-dependent loadings, H₂ population increases from 94.96 to 653.57 molecules and D decreases by 66.8%, while adsorption heat changes only modestly. Pore-wall type produces the largest diffusion contrast, with mean H₂ populations ranging from 83.26 to 100.39 molecules in the individual models. Heating lowers adsorption heat and H₂ population while increasing D, whereas pore width changes accessible space and confinement. Kerogen type produces smaller differences among the reported mean values. Within the investigated conditions, pressure-associated loading and pore-wall type therefore produced the largest observed variations in the coupled storage–diffusion quantities. The temperature, pressure, and pore-width intervals and the categorical material comparisons are not directly commensurate; these observations therefore describe the tested cases rather than a universal or normalized ranking of factor sensitivity.

TABLE 3. Observed responses of adsorption heat, H₂ population, and the directionally averaged self-diffusion coefficient to the five investigated factors. Ranges are taken from Figures 5–9; population is expressed in molecules and reflects the loading procedure in Section 2.2. The numerical and categorical ranges differ, so these descriptive comparisons do not constitute a normalized sensitivity ranking.

The different descriptor rankings are a central outcome of the comparison. Q a d s describes mean energetic binding per H₂ molecule, H₂ population describes occupancy of the specified pore region under the prescribed loading, and D describes their directionally averaged translational diffusion. Quartz has the highest Q a d s, whereas calcite combines the largest population with the lowest D; graphite exhibits the fastest diffusion. These observations show that strong energetic binding alone does not establish either the largest pore occupancy or the greatest diffusion resistance. Interfacial interactions must therefore be considered alongside accessible pore space and confinement when interpreting storage and transport. For shale evaluation, these model-specific relationships motivate joint consideration of pore-size distribution, mineral composition, kerogen type, and operating pressure, while retaining the distinction between molecular-scale retention and field-scale injection–withdrawal performance.

Formation water is a particularly important extension because hydrated mineral surfaces need not retain the gas-accessible interface of a dry pore. Zhang et al. (2027) studied CO₂-induced water-film reconstruction in hydrated shale-related nanoslits containing CO₂, CH₄, and He. Their results associate mineral-dependent hydration layers and asymmetric water films with differences in gas access to the walls. Although that system does not contain H₂, it provides a useful interfacial framework for examining whether hydration changes the pore-wall trends observed here. H₂–brine–mineral simulations also show sensitivity to salinity, mineral surface structure, and organic functional groups (Chen et al., 2023). Future hydrated H₂ models should therefore consider water-film distribution and gas-accessible pore space alongside mineral composition.

Future studies incorporating formation water, residual CH₄ and CO₂, surface roughness, flexible pore walls, and multicomponent interactions can further clarify how adsorption heat, pore population, and diffusion evolve during repeated injection–withdrawal cycles. Coupling these molecular descriptors with pore-network and reservoir-scale models will strengthen multiscale evaluation of shale formations for underground hydrogen storage. Pore geometry, connectivity, and geomechanical characterization are also required for reservoir-scale transport assessment (Wood and Hazra, 2017; Schwartz et al., 2019; Mu et al., 2025). At the UHS scale, Mu et al. (2025) further emphasize the joint assessment of lithology, water, and gas-flow performance when screening depleted reservoirs.

4 Conclusion

The molecular simulations lead to three main conclusions:

• At 323 K and 5 MPa in the 5 nm quartz slit pore, H₂ exhibits pore-center filling with modest near-wall density modulation. The adsorption heat is 0.67 kcal mol^-1, the mean H₂ population is 94.96 molecules, and the directionally averaged MSD-derived D is 409.32 Ų ps^-1. These complementary quantities describe energetic binding, pore occupancy, and molecular motion under a common reference loading.
• Heating increases diffusion while reducing adsorption heat and population under the prescribed initialization protocol; higher-pressure loading states contain more molecules and show lower diffusion. Narrowing the pore reduces the number of accommodated molecules and increases diffusion resistance, with the largest discrete change between 2 and 3 nm in the investigated quartz models. The population changes reflect the loading protocol as well as accessible pore space and do not constitute an independent prediction of equilibrium adsorption capacity.
• Pore-wall type produces nonidentical rankings of energetic binding, population, and diffusion: quartz has the highest adsorption heat, calcite the largest population and lowest D, and graphite the highest D. Differences among kerogen types are modest, with type III showing a tendency toward stronger retention. Joint interpretation of these descriptors provides a more informative basis for comparing shale pore models than any single adsorption or diffusion measure. Pressure-associated loading and pore-wall type produced the largest observed variations within the investigated conditions, rather than establishing a universal ordering of controlling factors.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Tao Hu, China University of Petroleum, Beijing, China

Reviewed by: Fuhua Shang, Inner Mongolia University of Technology, China

Hexin Huang, Chang’an University, China

Liyuan Zhang, Lomonosov Moscow State University, Russia