Experimental Measurements

Adsorption isotherms at 273 K and 300 K were measured on a Micromeritics ASAP 2020 surface area analyzer (Norcross, GA, USA). Clinoptilolite, > 97% purity, was purchased from KMI Zeolite, Inc. (Amargosa Valley, NV, USA) and originated in Nevada, USA. The sample was powdered to an average grain size of 76 μm. A water circulation bath was used to maintain the constant temperature of 300 K for the sample tube. A mixture of ice and water was used to achieve 273 K in the analytical bath. Before each measurement, the powders were degassed and vacuumed for a minimum of 4 h at room temperature. The measurement consisted of sequential quantification of adsorbed gas as a function of pressure up to 1 bar. Adsorption equilibrium was assumed when the differences from consecutive pressure readings were within 5% or 5 torr, whichever was smaller. Typically, the time between successive pressure readings were set at 10 s. An additional measurement for Kr at 300 K with successive readings set at 1000 s was made to determine if adsorption equilibrium had been achieved.

Bulk Zeolites

Classical grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulations were performed on bulk zeolite models for hydrated clinoptilolite and mordenite. Model systems were prepared from the published structures (Simoncic & Armbruster, Reference Simoncic and Armbruster2004; Smyth et al., Reference Smyth, Spaid and Bish1990) corresponding to three hydration states: fully hydrated, partially hydrated (half of the waters removed), and dry (all waters removed). Cation and water (oxygen) positions from the published structures were used in the initial configurations of the simulation models.

The clinoptilolite model contains the same number of Na⁺, K⁺, and Ca²⁺, corresponding to the chemical formula (Ca,Na,K)₆Al₆Si₃₀O₇₀⋅24H₂O (Smyth et al., Reference Smyth, Spaid and Bish1990). The C2/m unit cell was orthogonalized with the crystallographic b and c axes aligned with z and x, respectively. A 4 × 2 × 2 supercell was created with xyz dimensions 29.60 Å × 31.59 Å × 35.88 Å. Water contents for fully hydrated and partially hydrated clinoptilolite models were 15.5 wt. % and 7.6 wt. %, respectively. Cross-section models of the a and c axes and the corresponding pore size distributions are shown in Fig. 1. The a-pore is very irregular with several local minima at diameters <3.0 Å, corresponding to windows connecting adjacent cages. This implies a high diffusion barrier for all gases in this direction. There are two pores along the c-direction, each with a minimum pore size of ~3.3 Å. This also suggests a high diffusion barrier, even for Ar which has a slightly larger kinetic diameter of 3.4 Å (Li et al., Reference Li, Kuppler and Zhou2009).

Fig. 1 Models of the silicate framework for clinoptilolite (Si = yellow, O = red) showing pores parallel to the (A) a and (B) c axes. Pore size profile from the HOLE program (Smart et al., Reference Smart, Neduvelil, Wang, Wallace and Sansom1996) along the (C) a and (D) c axes. Extra-framework cations have been removed for clarity

The mordenite unit cell is orthorhombic with cation sites occupied by Na⁺ corresponding to the synthetic sample with chemical formula Na₆Al₆Si₄₂O₉₆⋅28H₂O (Simoncic & Armbruster, Reference Simoncic and Armbruster2004). A 2 × 1 × 3 supercell was created with xyz dimensions 36.02 Å × 20.51 Å × 22.57 Å. Water contents for the fully hydrated and partially hydrated mordenite models were 11.1 wt. % and wt. 5.7%, respectively. The pore structure in mordenite (Fig. 2) indicates that diffusion in the b-direction requires a zig-zag path between cages, with a minimum pore diameter of 3.6 Å at the cage window. Two pores exist along the c-direction, but only the larger pore is accessible to the gases. In fact, gas transport through mordenite is thought to occur through the large c-pore (Viswanadham & Kumar, Reference Viswanadham and Kumar2006), which has a consistent pore diameter of 6.5 Å.

Fig. 2 Stick models of the silicate framework for mordenite (atoms colored as in Fig. 1) showing pores parallel to the (A) b and (B) c axes in the bc and ab planes. The pink spheres in A show the diffusion path through two pores parallel to the b axis. The gray lines in B represent the unit cell. Pore size profile from the HOLE program (Smart et al., Reference Smart, Neduvelil, Wang, Wallace and Sansom1996) along the (C) b and (D) c axes. Extra-framework cations have been removed for clarity

All simulations were performed using the LAMMPS code (Thompson et al., Reference Thompson, Aktulga, Berger, Bolintineanu, Brown, Crozier, in’t Veld, Kohlmeyer, Moore, Nguyen, Shan, Stevens, Tranchida, Trott and Plimpton2022) with three-dimensional (3D) periodic boundary conditions at a temperature of 300 K. Forces and potential energies were computed from pairwise atomic interactions based on van der Waals (vdW) and electrostatic terms with a short-range cutoff distance of 10.0 Å. Zeolite framework atoms (Si, Al, O) remained fixed at their crystallographic coordinates throughout the GCMC and bulk MD simulations, but framework flexibility was allowed in the slab MD simulations. Electrostatic interactions between all atoms are based on partial charges with the long-range component computed using a particle-mesh Ewald algorithm (Plimpton et al., Reference Plimpton, Pollock and Stevens1997). Rather than assigning specific sites for tetrahedral Si and Al atoms, a single tetrahedral atom type was used with the atomic charge corresponding to the net negative charge for each zeolite. This approach has been used by Boutin and co-workers to model several charged zeolites (Buttefey et al., Reference Buttefey, Boutin, Mellot-Draznieks and Fuchs2001; Coudert et al., Reference Coudert, Cailliez, Vuilleumier, Fuchs and Boutin2009; Di Lella et al., Reference Di Lella, Desbiens, Boutin, Demachy, Ungerer, Bellat and Fuchs2006; Jeffroy et al., Reference Jeffroy, Boutin and Fuchs2011, Reference Jeffroy, Nieto-Draghi and Boutin2014). There were no vdW interactions between pore species and zeolite Si atoms. The Buckingham potential was used for cation-O_zeolite vdW interactions with parameters taken from previous work (Jeffroy et al., Reference Jeffroy, Nieto-Draghi and Boutin2014, Reference Jeffroy, Nieto-Draghi and Boutin2017). All other vdW interactions were based on the Lennard–Jones potential, with parameters taken from previous work: cation-cation (Jeffroy et al., Reference Jeffroy, Boutin and Fuchs2011, Reference Jeffroy, Nieto-Draghi and Boutin2014), flexible SPC water (Teleman et al., Reference Teleman, Jonsson and Engstrom1987), cation-water (Aqvist, Reference Aqvist1990; Dang, Reference Dang1995), and noble gases (Pellenq & Levitz, Reference Pellenq and Levitz2002; van Loef, Reference van Loef1981). Lorentz-Berthelot combination rules were used to calculate vdW parameters for interactions between unlike atoms. Intramolecular interactions included bond stretch and angle bending terms for SPC water (Teleman et al., Reference Teleman, Jonsson and Engstrom1987) and the zeolites (Jeffroy et al., Reference Jeffroy, Nieto-Draghi and Boutin2014). Silanol (SiOH) groups at zeolite slab terminations were described by the SPC bond stretch term and a Si–O-H angle bending term previously used for silica surfaces (Harvey & Thompson, Reference Harvey and Thompson2015). All potential parameters are included in the Supplementary Information (SI).

Equilibrium configurations of dry and hydrated zeolites were obtained from 1.5-ns constant volume MD simulations with a timestep of 0.5 fs for short-range interactions and 1.0 fs for long-range electrostatic interactions. Noble gas adsorption isotherms were obtained from GCMC simulations, in which gas atoms are randomly inserted into the bulk zeolite model based on acceptance criteria. The initial zeolite configuration was taken from the final MD frame of the corresponding pure dry or hydrated zeolite. Chemical potential values for the noble gases corresponded to pressures ranging from 0.1 atm to ~100 atm as determined from separate pure gas simulations (Perry et al., Reference Perry, Teich-McGoldrick, Meek, Greathouse, Haranczyk and Allendorf2014). Each simulation consisted of 1.5 × 10⁷ steps, where each step involved 20 gas insertion/deletion moves and 10 gas translation moves. Average adsorption loadings were obtained from the final 5 × 10⁶ steps.

Gas mobility in dry and partially wet zeolites was investigated from 10-ns constant-volume MD simulations at fixed gas loadings ranging from 0.4 to 2.3 atoms per unit cell. The actual number of gas atoms in the supercells ranged from 5 to 14, which provided adequate statistics for analysis. Initial configurations were taken from GCMC snapshots with some gas atoms removed to achieve the desired loading. Also, some Ar atoms in the dry mordenite model were moved from the small c-pores to the large c-pores (discussed below), as Ar mobility in the former is severely limited. Based on our analysis of the pore-size profiles (Fig. 2), it is unlikely that any Ar atoms would be able to enter the small c-pores. Specific gas loadings and approximate pressures are given in Tables S1 and S2 of the SI. Framework zeolite atoms (Si and O) remained fixed, but extra-framework cations and gas species were mobile with velocities controlled by a Nosé-Hoover thermostat with a relaxation time of 100 fs.

Zeolite Slabs

Unlike the GCMC simulations described above in which gas species are directly inserted into the bulk zeolite structure, in the slab models the gas species enter the zeolite pores through an external surface. Zeolite surfaces were created by cleaving the bulk crystal from the corresponding pure zeolite MD simulation and converting undercoordinated surface sites to isolated silanol (SiOH) groups. As seen in Fig. 3, the MD models are periodic in the two directions parallel to the surfaces, and nonperiodic in the perpendicular direction. The diffusion pathways (nonperiodic direction) correspond to the c-axis of both clinoptilolite and mordenite, as these channels have the largest pore-limiting diameters in these minerals (Fig. 1 and Fig. 2). Simulations from models oriented along other channels resulted in little or no gas uptake. Expanded supercells were used so that each lateral dimension of the external surface was a least 35 Å. Surface dimensions in the ab planes for slab models were 31.59 Å × 35.88 Å for clinoptilolite (3 × 2 supercell) and 36.02 Å × 41.01 Å for mordenite (2 × 2 supercell). Models were created for dry zeolite with and without extra-framework cations, and partially hydrated zeolite (7.6 wt.% and 5.7 wt.% for clinoptilolite and mordenite, respectively). For zeolite models without extra-framework cation removed, Si and O charges were adjusted to 1.4 e and − 0.7 e, respectively (e is the elementary charge), so that the zeolite framework was charge neutral (i.e., siliceous zeolite).

Fig. 3 MD snapshots of krypton gas and A 4.4-nm thick clinoptilolite slab and B 3.8-nm thick mordenite slab with the zeolite c-pores oriented in the horizontal direction. Black lines denote the simulation cell boundaries. Zeolite slabs are represented as stick models with Si and O atoms colored as in Fig. 1. Additional atom colors are as follows: Na = blue, Ca = green, K = purple, H = white, Kr = light blue, and C = gray. Insets show surface SiOH groups

Initial configurations for the two-phase simulations consisted of zeolite slabs with a 50 Å vacuum gap on either side of the zeolite surface. Three slab models were created for clinoptilolite (43.7 Å, 58.6 Å, and 73.8 Å which are equivalent to 6, 8, and 10 unit cells in the c-direction, respectively) to determine the effect of slab thickness on gas adsorption kinetics. A single slab model was created for mordenite (37.9 Å equivalent to 8 unit cells in the c-direction) because all gases rapidly adsorb and diffuse in this zeolite. For single gas simulations, gas atoms were randomly inserted in each vacuum region at low concentration (50 atoms on each side) or high concentration (100 atoms on each side). For mixture gas simulations (Ar/Xe and Kr/Xe), 50 atoms of each gas were randomly inserted on each side. Based on the pure gas properties using the same force field parameters (Perry et al., Reference Perry, Teich-McGoldrick, Meek, Greathouse, Haranczyk and Allendorf2014), these initial gas concentrations correspond to pressures of 15–20 bar (50 atoms on each side) to 30–40 bar (100 atoms on each side). Corresponding gas pressures at the end of each simulation can be obtained by comparing steady-state gas loadings in zeolite slabs with loadings from adsorption isotherms (discussed below). A wall consisting of a single layer of graphite C atoms (Trucano & Chen, Reference Trucano and Chen1975) prevented gas atoms in the two vacuum regions from comingling. The C atoms remained fixed during the simulations, but they interacted with gas atoms via vdW interactions as effective hard spheres (0.0001 kcal⋅mol⁻¹).

Constant volume MD simulations were performed at temperatures of 250 K, 300 K, and 350 K. Separate Nosé-Hoover thermostats were used for the zeolite and gas velocities, each with a relaxation time of 100 fs. Separate thermostats have been used in previous simulation studies of two-phase systems (Inzoli et al., Reference Inzoli, Simon and Kjelstrup2009; Simonnin et al., Reference Simonnin, Marry, Noetinger, Nieto-Draghi and Rotenberg2018), and it was confirmed that the total temperature matched the target temperature. To avoid drift, the center of mass of the zeolite framework was constrained using the ‘fix recenter’ command in LAMMPs. The length of the simulations varied from 40 ns up to 1000 ns depending on how rapidly steady state loading was achieved. For each model system, results were averaged from ten independent simulations with different initial configurations.

Kinetic Models of Diffusion

To estimate both the quantity of the adsorbed gas and the diffusivity in the zeolite, a kinetic model was adapted from earlier works (Barrer, Reference Barrer1949; Crank, Reference Crank1975; Wilson, Reference Wilson1948). Whereas Karger and Ruthven refer to this as the diffusion-immobilization model (Karger & Ruthven, Reference Karger and Ruthven2016), the mathematical solution does not require immobilization but only that the transport diffusivity, $D_t$ , is constant and uniform throughout the slab. Here, the transport diffusivity is the product of the Maxwell–Stefan diffusivity and the thermodynamic factor that relates chemical potential to concentration (Krishna, Reference Krishna2012). Let $n$ be the number density or volumetric mass concentration in the pores that includes no distinction between free and adsorbed particles. Let $x$ be the distance from the centerline of the slab in the $c$ -coordinate and $t$ be the time elapsed in the simulation. Due to the periodicity of the system, the concentration was assumed to be uniform in the $a$ and $b$ directions. In this idealized system, the mass balance within the slab is a restatement of Fick’s Second Law in one dimension.

(1) $\frac{\partial n}{\partial t}=D_t\frac{{\partial }^{2}n}{\partial x^2}$

(2) $n\left(x,,,0\right)=0$

(3) ${\frac{\partial n}{\partial x}}_{x=0}=0$

(4) $n\left(w,,,t\right)=Kf\left(t\right)$

The zeolite slab is initially devoid of the gas; thusly the initial condition is zero concentration throughout the slab (Eq. 2). The first boundary condition (Eq. 3) is the no-flux condition at the centerline $x=0$ due to the symmetry of the model. The second boundary condition (Eq. 4) requires the concentration at exterior edges of the slab, that is at distance $x=w$ , to be a function $f$ of the free gas concentration. Due to the lower chemical potential in the slab, the concentration within the slab at equilibrium will be larger than the free gas concentration. To a first approximation, the concentration within the slab pores is related to the concentration in the void by the equilibrium constant $K$ .

Rather than solving the preceding equation directly for the concentration in the slab, instead this work followed the method proposed by Wilson (Reference Wilson1948) and seeks the total mass of the species within a subset of the slab starting at centerline. Let the linear mass loading function be described as $m\left(x,,,t\right)$ .

(5) $m\left(x,,,t\right)=\frac{{\int }_0^z{\int }_0^y{\int }_0^{x}n\left(x^{\prime },,,t\right)d{x}^{\prime }d{y}^{\prime }d{z}^{\prime }}{{\int }_0^z{\int }_0^{y}d{y}^{\prime }d{z}^{\prime }}=\varphi {\int }_0^{x}n(x^{\prime },t)d{x}^{\prime }$

where $y$ ’ and $z$ ’ represent integrals in the a and b coordinates (Fig. 3). Because only a fraction of the slab region consists of the accessible pore volume, but the tallies were conducted over the entire volume, the mass loading function is function of both the pore concentration and the areal porosity $\varphi$ .

Performing this integral transformation on the earlier differential equations (Eq. 1) results in the similar differential equation but altered boundary conditions:

(6) $\frac{\partial m}{\partial t}=D_t\frac{{\partial }^{2}m}{\partial x^2}$

(7) $m\left(x,,,0\right)=0$

(8) $m\left(0,,,t\right)=0$

(9) ${\frac{\partial m}{\partial x}}_{x=w}=\varphi n\left(w,,,t\right)=\varphi Kf\left(t\right)$

The second boundary condition (Eq. 9) now relates the $x$ -derivative of the mass-loading function $m$ to the free gas concentration $f$ . An expression for the free gas region concentration was then constructed using the overall mass balance on the system.

In these simulations, the mass loading inside the slab was twice the mass in a singular side; that is $2m\left(a,,,t\right)$ . Defining the linear extent of the void region as $l$ , the mass of free gas outside of the slab is twice the free gas concentration and linear extent; that is $2fl=\frac{2l}{\varphi K}{\frac{\partial m}{\partial x}}_{x=a}$ . As recognized by Barrer (Reference Barrer1949), a non-negligible quantity of gas may adsorb as a film on the external surface of the slab. In this work, it was assumed that this surface is saturated and hence constant and independent of the gas concentration in the other regions. Thus, the total mass within the slab and void regions was the total mass of the gas $M_{\mathrm{total}}$ minus the mass in the external film $M_{\mathrm{film}}.$ In doing so, the external boundary condition (Eq. 9) can be replaced with the following:

(10) $m\left(w,,,t\right)+\frac{l}{\varphi K}{\frac{\partial m}{\partial x}}_{x=w}=\frac{M_{\mathrm{total}}-M_{\mathrm{film}}}{2}.$

The solution to this partial differential equation can be concisely written using two parameters to non-dimensionalize the system. Firstly, the capacity ratio $\alpha$ is the ratio of the capacity of the void versus the adsorptive capacity of the slab:

(11) $\alpha =\frac{l}{\varphi Kw}$

Secondly, let the kinetic parameter $\beta$ be defined as the ratio of the transport diffusivity to the width of the slab squared.

(12) $\beta =\frac{D_t}{w^2}$

To solve this heterogeneous but constant boundary condition, the solution to the mass-loading function $m$ was modeled as the summation of the steady-state and transient solutions:

(13) $m\left(x,,,t\right)=s\left(x\right)+u\left(x,,,t\right).$

where the steady-state solution is a linear function of the number density at equilibrium $n_{\infty }$ :

(14) $s\left(x\right)=n_{\infty }x=\frac{M_{\mathrm{total}}-M_{\mathrm{film}}}{2\left(1,+,\alpha \right)w}x$

and the transient solution is composed of the general solution, where $A_i$ and $B_i$ are undetermined coefficients and $q_n$ are eigenvalues the satisfy the second mixed boundary condition:

(15) $u_i\left(x,,,t\right)=e^{-\beta q_i^{2}t}\left(A_i,\sin ,\left(\frac{q_{n}x}{w}\right),+,B_i,\cos ,\left(\frac{q_{n}x}{w}\right)\right)$

Applying the two boundary conditions requires all $B_n$ coefficients to be zero and that the eigenvalues $q_n$ satisfy the transcendental equation.

(16) $\tan q_i+\alpha q_i=0$

The associated coefficients $A_i$ were then found by eigenfunction expansion on the initial condition after subtracting the steady state solution.

(17) $A_i=\frac{M_{\mathrm{total}}-M_{\mathrm{film}}}{1+\alpha +{\alpha }^2{q}_i^2}\frac{1}{q_i\cos q_i}$

Thus, the complete solution for the mass-loading function is:

(18) $m\left(x,,,t\right)=\frac{M_{\mathrm{total}}-M_{\mathrm{film}}}{2\left(1,+,\alpha \right)}\frac{x}{w}+{\sum }_{i=1}^{\infty }\frac{M_{\mathrm{total}}-M_{\mathrm{film}}}{1+\alpha +{\alpha }^2{q}_i^2}\frac{\sin \left(q_i,\frac{x}{a}\right)}{q_i\cos q_i}{e}^{-\beta q_i^{2}t}$

To determine the values of the parameters $\alpha$ , $\beta$ , and $M_{\mathrm{ext}}$ , it is prudent to fit the total mass within the slab, which maximizes the number of particles counted, reducing variance due to discrete fluctuations in the tallies. Having a mass-loading function for the half-slab, the quantity within the tallied region of the simulation is simply twice the mass in the half-slab plus the mass with the thin film on the interface. Let the mass loading function for the tallied region be $M_{\mathrm{tally}}\left(t\right)$ .

(19) $M_{\mathrm{tally}}\left(t\right)=M_{\mathrm{film}}+\frac{M_{\mathrm{total}}-M_{\mathrm{film}}}{1+\alpha }\left(1,-,{\sum }_{i=1}^{\infty },\frac{2\alpha \left(1,+,\alpha \right)}{1+\alpha +{\alpha }^2{q}_i^2},e^{-\beta q_i^{2}t}\right)$

To fit this function, a subroutine in MATLAB (v R2021b 9.11.0.1769968) was programmed both to evaluate the first 1000 eigenvalues $q_n$ for a given $\alpha$ , and calculate the time-dependent mass loading for the specified $\beta$ and $M_{\mathrm{ext}}$ . The non-linear least squares parameter estimation routine in MATLAB was utilized to perform the parameter fits. All fitted parameters and their corresponding 95% confidence intervals are reported in the SI.

Series truncation error is small when the non-dimensional time, that is $\beta t$ , takes on modest values of 0.01 or greater (Crank, Reference Crank1948). However, a problem can arise when the preponderance of adsorption occurs at small values of time. Wilson (Reference Wilson1948) showed that, when the infinite series summation becomes impractical, an alternative representation is more effective. In this small-value time approximation, the system was dominated by a different characteristic timescale $\gamma =\frac{\beta }{{\alpha }^2}$ . That is, early in the process, the quantity entering the slab is independent of slab thickness as few if any particles reach the centerline within the span of the simulation. Consequently, it is not possible to make independent estimates of $\alpha$ and $\beta$ in this regime. The mass loading function using this small-value time approximation is as follows:

(20) $\underset{\beta t\to 0}{\lim }{M}_{\mathrm{tally}}\left(t\right)=M_{\mathrm{film}}+{(M}_{\mathrm{total}}-M_{\mathrm{film}})\left(1,-,e^{-\gamma t},\mathrm{erfc},\sqrt{\gamma t}\right).$

In the special case of an infinite reservoir, and hence infinite $\alpha ,$ this function appears as the more commonly referenced $\sqrt{t}$ ‘law’. However, for the limited extent of the void region here, this formulation is more necessary. Where the small-value time approximation is necessary, it is not possible to isolate the capacity ratio $\alpha$ from the kinetic parameter $\beta$ , however, the characteristic timescale can, conversely, be calculated from independent estimates from the large-value time formulation, that is Eq. 17. In doing so, results can be compared across all simulations, even where the small-value time approximation is necessary.