Densified HKUST-1 Monoliths as a Route to High Volumetric and Gravimetric Hydrogen Storage Capacity

We are currently witnessing the dawn of hydrogen (H2) economy, where H2 will soon become a primary fuel for heating, transportation, and long-distance and long-term energy storage. Among diverse possibilities, H2 can be stored as a pressurized gas, a cryogenic liquid, or a solid fuel via adsorption onto porous materials. Metal–organic frameworks (MOFs) have emerged as adsorbent materials with the highest theoretical H2 storage densities on both a volumetric and gravimetric basis. However, a critical bottleneck for the use of H2 as a transportation fuel has been the lack of densification methods capable of shaping MOFs into practical formulations while maintaining their adsorptive performance. Here, we report a high-throughput screening and deep analysis of a database of MOFs to find optimal materials, followed by the synthesis, characterization, and performance evaluation of an optimal monolithic MOF (monoMOF) for H2 storage. After densification, this monoMOF stores 46 g L–1 H2 at 50 bar and 77 K and delivers 41 and 42 g L–1 H2 at operating pressures of 25 and 50 bar, respectively, when deployed in a combined temperature–pressure (25–50 bar/77 K → 5 bar/160 K) swing gas delivery system. This performance represents up to an 80% reduction in the operating pressure requirements for delivering H2 gas when compared with benchmark materials and an 83% reduction compared to compressed H2 gas. Our findings represent a substantial step forward in the application of high-density materials for volumetric H2 storage applications.

Adsorption experiments (up to 1 bar) for different pure gases were performed on Micromeritics 3Flex surface area and pore size analyzer ( Fig. S2-S3). About 200 mg of activated samples were used for the measurements. Liquid nitrogen was used to maintain constant temperature in the bath through the duration of the experiment. Samples were degassed on a Smart VacPrep prior to the analysis.

Molecular Simulations
The crystallographic information files (cif) of the 2932 MOFs studied in this work are taken from the DDEC database 2 . We also included 10 MOFs with benchmark performance for H2 storage in our high-throughput screening study -HKUST-1, MOF-5, NU-100/PCN-610, NU-1501-Al, NU-1500-Al, Ni(dobdc), MIL-101, IRMOF-10, UMCM-9, and IRMOF- 20. More details about the cif of these benchmark materials are given in Table S1. The low pressure adsorption of N2 and absolute adsorption capacity for H2 at 77 K is simulated using the Grand Canonical Monte Carlo method as implemented in the RASPA simulation package 3 where is the depth of the potential well, is the distance at which the intermolecular potential between the two particles is zero, and is the distance between atom and atom . Lorentz Berthelot mixing rules are used for all the cross-interaction terms: where and are the LJ parameters (as described above) and the subscripts indicate atom and atom . Electrostatic interactions are modeled using the coulombic potential (3) and are computed using the Ewald summation method with the precision set to 10 -6 : where, is the charge, 0 is the permittivity, and is as defined above. N2 is modeled using parameters taken from the TraPPE force field 4 . Table S2 shows the parameters used for N2 in this study. H2 is modeled as a rigid three-site model with parameters taken from the work of Bucior et al. 5 . Table S3 shows the parameters used for H2 in this study. To account for the quantum effects, which are important at cryogenic temperatures 6 , the LJ interactions between H2-H2 and H2-framework are modified to include the Feynman-Hibbs correction 7 where ℏ is the planck's constant, is the Boltzmann constant, = /( + ) is the reduced mass, and is the temperature. The framework atoms are modeled using LJ parameters taken from both the DREIDING force field 8 and the Universal force field 9 . Table S4 shows the LJ parameters used to describe the atoms in the framework. The framework is modeled as a rigid structure with periodic boundary conditions applied in all directions. The number of unit cells in the simulation box is computed individually for each of the frameworks to make sure the simulation box is large enough to ensure that a distance of at least twice the cutoff radius is maintained between the periodic images.

Principal Component Analysis (PCA)
We performed Principal Component Analysis (PCA) using the data generated from the highthroughput screening of 2940 MOFs at the following conditions: adsorption at 25/50/100 bar and 77 K with desorption at 5 bar and 160 K. The raw data used for the PCA can be found at: https://aam.ceb.cam.ac.uk/mofexplorer.html -please refer to Part 1 of the dataset. PCA was performed using the Principal Component Analysis Visualization Tool which can be found at 16 : https://hydrogen-storage-pca.herokuapp.com, considering the following properties: pressure (bar), density (g cm -3 ), gravimetric deliverable capacity (GDC, g g -1 ), void fraction, accessible surface area per mass (m 2 g -1 ), accessible surface area per volume (m 2 cm -3 ), largest cavity diameter (LCD, Å), volumetric deliverable capacity (VDC, g cm -3 ), and the pore limiting diameter (PLD, Å). In this work, PCA was based on a correlation matrix. Figure S8 shows the correlation plot for the selected geometric properties and pressure. We can see that the LCD, VDC, Acc. SA per mass, Acc. SA per volume, void fraction, and GDC are highly correlated among themselves. This is not surprising as these geometric properties are just different ways of characterizing the free space within the material. The pressure is not correlated to any of the geometric properties considered. However, the density is inversely correlated to the LCD, VDC, Acc.
SA per mass, Acc. SA per volume, void fraction, and GDC. This is not surprising either, as a material with larger pores, will have a lower density.

S16
PCA is commonly used for dimensionality reduction, helping us to choose the minimum number of variables needed to explain the maximum amount of variance in the dataset. Figure S9 shows the cumulative scree plot, with the region in red depicting that a minimum of 4 principal components are needed to explain ~ 93% of the variance seen in the dataset. Using 1 principal component would help us to explain only ~70 % of the variance in the dataset. Hence, for further analysis we pick 4 principal components. The relationship between the properties and the principal components can be quantified/measured through loadings and contributions.  Figure S10 shows the % contribution of each feature to the 4 different principal components selected (PC1, PC2, PC3, and PC4). We can see that the GDC, void fraction, Acc. SA per mass, Acc. SA per volume, LCD, VDC, and PLD all exhibit strong positive loadings (Table S7), contributing in the range of 10-14 % each to PC1 ( Figure   S10). This means that an increase in any of these geometric properties would most likely lead to an increase in performance. The density on the other hand exhibits strong negative loadings, and contributes ~10% to PC1. This means that any further increase in density will most likely lead to a decrease in performance. The pressure exhibits very strong negative loadings, contributing ~ 98% to PC2. This implies that any change in PC2 scores for a material is primarily due to the influence of pressure. Figure S9. Cumulative scree plot of the explained variance as a function of the number of principal components. The region shaded in red depicts that 93% of the variance in the data can be explained using 4 principal components.   scores. This implies that as the pressure increases from 50 to 100 bar, the pressure starts to negatively influence the performance of the materials. This tells us that there is an optimum pressure at which the performance of the materials peaks. This is most likely in the range of 50-55 bar, i.e. when most materials in the dataset have PC2 scores close to 0. We also analyzed the top 1 % of materials and found that the 50 bar dataset was best represented in that 90% of the structures present in the top 1 % of the 50 bar dataset were also present in either/both the top 1 % of the 25 bar and 100 bar datasets, with the top 5 materials being same for both the 50 bar and 100 bar datasets. We hope these insights will lead to the rational selection of process conditions for H2 storage in benchmark MOF materials.

NREL high-pressure gas adsorption measurements
The H2 gas adsorption measurement at -197 °C at higher pressures (~5 to 105 bar Figure S12) was collected on a volumetric adsorption apparatus at the National Renewable Energy Laboratory (NREL). (The temperature of liquid N2 at the elevation in Golden, CO is 75.6 K). The sample was packed in the custom-made sample holder in a dry box and transported to NREL in argon atmosphere.
The sample holder includes a VCR diaphragm isolation valve, 1/8" tubing, and a 1/4" Swagelok VCR union and cap where the sample is located. There is a 2 micron filter gasket to ensure the sample remains in the sample holder. The sample was degassed to 120 °C at a heating rate of 2 °C/min, held at 120 °C for 12 hours, and cooled to ambient temperature at 2 °C/min. The H2 adsorption isotherm was measured on a custom modified PCTPro 2000. The mass of the sample after measuring the isotherm was 248.5 mg. This system has been previously described in detail 17 . Modifications to the system include water-cooled copper jacket fittings to maintain extended temperature control of the gas dosing arm and the sample holder. For measurements at -197 °C, the bottom portion of the sample 2 holder containing the sample is placed in a liquid nitrogen bath. It is essential that the level of the liquid nitrogen bath remain constant throughout the measurement, as to not change the temperature profile of the sample holder as this would change the headspace calibration for the instrument. The same experimental setup was used to measure the isotherm at 303 K (1 to 140 bar, Figure S13) with the exception that a cryostat was used to control the temperature of the sample and the sample mass was 248.5 mg. It is important to note that the NREL experimentally measured values are excess amounts adsorbed (Nexc), which were transformed into total uptakes (Ntot) by using equation (5): where ρ is the density of the gas at the given adsorption pressure and temperature, obtained from the National Institute of Standards and Technology (NIST) 18 , and Venvelope is the envelope volume and Vskeletal is the He skeletal density of the material 19 .
Open symbols represent excess uptake data, while closed symbols represent the calculated total uptake. Total adsorption values were calculated using the NREL method (Equation 5). An envelope density (ρenvelope) of 1.07 g ml -1 was used to calculate the volumetric H2 uptake of the monoHKUST-1 material, while a packing density of 0.2979 g cm -3 was used to calculate the volumetric H2 uptake of the powdered HKUST-1 material. Figure S13. Excess (Nexc) and total (Ntot) H2 adsorption isotherms collected by NREL for monoHKUST-1 and powdered HKUST-1 at 303 K. Open circles represent excess uptake data, while closed circles represent the calculated total uptake. Total adsorption values were calculated using the NREL method (Equation 5). An envelope density (ρenvelope) of 1.07 g ml -1 was used to calculate the volumetric H2 uptake of the monoHKUST-1 material, while a packing density of 0.2979 g cm -3 was used to calculate the volumetric H2 uptake of the powdered HKUST-1 material. analyses, the samples were activated overnight at 120 °C (vacuum) before measuring the mass, and then degassed in situ thoroughly before the gas adsorption. It is important to note that the experimentally measured values are excess amounts adsorbed (Nexc), which are transformed into absolute uptakes (Nabs) by using equation (6): where ρ is the density of the gas at the given adsorption pressure and temperature, obtained from the National Institute of Standards and Technology (NIST) 18 , and Vpore is the pore volume of the adsorbent 20 . Figure S14. Absolute (Nabs) H2 adsorption isotherms collected at the University of Cambridge for monoHKUST-1 at 77, 160, 233 and 273 K. Absolute adsorption values were calculated using Equation 6. An envelope density (ρenvelope) of 1.064 g ml -1 and total pore volume (Vpore) of 0.634 cm 3 g -1 was used to calculate the volumetric H2 uptake of the monoHKUST-1 material. Figure S15. Excess (Nexc) and absolute (Nabs) H2 adsorption isotherms collected at the University of Alicante for monoHKUST-1 at 77, 195 and 298 K. Open circles represent excess uptake while closed circles represent absolute uptake. Absolute adsorption values were calculated using Equation 6. An envelope density (ρenvelope) of 1.064 g ml -1 and total pore volume (Vpore) of 0.634 cm 3 g -1 was used to calculate the volumetric H2 uptake of the monoHKUST-1 material.

Dual-Process Langmuir Model
The single-gas DPL model describes the adsorption of gas i on a heterogeneous adsorbent that is composed of two homogeneous but energetically different sites. Assuming that the adsorbate-adsorbent free energy on each patch is constant, the amount adsorbed mi for component i is given by where m1 and b are respectively the saturation capacity and affinity parameter on site 1, m2 and d are respectively the saturation capacity and affinity parameter on site 2, and P is the absolute pressure. All of the assumptions of the Langmuir model apply to each patch, and the two patches do not interact with each other.
The Henry's law constant is simply the sum of those on each patch or site (i.e., [(m1b)site 1 + (m2d)site 2]). In this formulation, the saturation capacity for each component on each site is allowed to be different. The free energy or affinity parameter for two different sites is expressed as: where the subscript b and d represents the free-energy level of site 1 or 2, Q1 and Q2 are the adsorption energy of components on site 1 or 2, and bo and do are the pre-exponential factor or adsorption entropy of component of site 1 or 2. In these calculations, site 1 always denotes the higher adsorbate-adsorbent free energy and site 2 always denotes the lower adsorbate-adsorbent free energy. For single-gas adsorption, this always makes the free energy of site 1 higher than that of site 2. The fitting parameters used can be found in Table S16. Absolute adsorption values were calculated using excess adsorption collected at the University of Alicante and Equation 6. A total pore volume (Vpore) of 0.634 cm 3 g -1 was used to calculate the H2 uptake of the monoHKUST-1 material. Dual-Process Langmuir fitting values are presented in Table S16. Figure S25. Dual-Process Langmuir (Equation 7) calculated H2 isotherms and corresponding experimental absolute (Nabs) isotherms calculated using data collected by NREL at 303 K for monoHKUST-1. Absolute adsorption values were calculated using excess adsorption collected at NREL and Equation 6. A total pore volume (Vpore) of 0.634 cm 3 g -1 was used to calculate the H2 uptake of the monoHKUST-1 material. Dual-Process Langmuir fitting values are presented in Table S16.

Mercury porosimetry density measurement
Mercury porosimetry was obtained up to a final pressure of 2,000 bar using an AutoPore IV 9500 instrument from Micromeritics. This technique was used to estimate the envelope density of the monoHKUST-1 at atmospheric pressure. Prior to the analysis, all samples were activated overnight at 120 °C (vacuum) before measuring the mass, and then degassed in situ thoroughly before the mercury porosimetry. As a part of the NREL calculated measurements, the envelope density was confirmed to be 1.07 g/cm 3 by mercury porosity measurement performed by Particle Authority 21 .

Solid State Nuclear Magnetic Resonance (NMR)
HKUST-1 samples were activated at 120 °C in flowing N2 for 3-4 h. Samples were packed into 3.2 mm rotors inside a N2 filled glovebag. 13 C NMR experiments were then performed at 16.4 T with a sample magic angle spinning frequency of 20 kHz. In each case a spin echo pulse sequence was used with a single rotor period as the interpulse delay time, and with a recyle delay of 100 ms. Owing to the large range of 13 C chemical shifts for these paramagnetic samples, two separate spectra were acquired with different transmitter frequencies. 13 C NMR spectra are referenced relative to the tertiary (eft-hand) carbon resonance in adamantane at 38.5 ppm.

X-ray Total Scattering Studies
Synchrotron X-ray total scattering data was collected at beamline 11-ID-B of the Advanced Photon Source (APS), Argonne National Laboratory (ANL), IL, USA 23 , using a Perkin Elmer area detector and X-ray wavelength = 0.2115 Å. To allow 2D mapping, a monolithic sample of HKUST-1 was sliced into four parallel sections, each approximately 1 mm thick, and secured in the centre of a nylon washer between two sheets of Kapton™ tape ( Figure S27, Figure 4(a)). Data suitable for Bragg diffraction analysis were collected at a sample-to-detector distance of ca. 947 mm, while data suitable for pair distribution function studies were collected at a sample-to-detector distance of ca. 197 mm.
Using a 500 μm 2 X-ray beam, the cross sections were scanned both horizontally and vertically in steps of 500 μm to afford a map consisting of 500 μm 2 pixels. Data was collected across the entire sample holder, and the edges of the monolith sections were estimated based on the presence of (i) diffraction rings or (ii) PDF atom-atom distances consistent with Cu. A CeO2 diffraction standard was used to calibrate experimental geometry in both cases and all scattering images were integrated using GSAS-2 24 . Pair distribution functions, G(r), were extracted from the total scattering data within xPDFsuite 25 , with X-ray scattering measured for an empty sample holder used for background subtraction. G(r) S38 was calculated using scattering data in the range 0.1 Å -1 ≤ Q ≤ 22.7 Å -1 . Inspection of the diffraction patterns reveals a number of spurious peaks present around the edge of the mapped monoliths and are emphasised upon LeBail fitting of HKUST-1 lattice parameters to the data in Topas (Figure S28).To compare the relative presence of impurities, the normalized integral intensity of the spurious diffraction peak at Q = 0.57 Å -1 and the (222) peak (Q = 0.83 Å -1 ) of HKUST-1 were fitted using the cumulative trapezoid method as implemented in the Python package scipy.integrate. (Figure 4(e-h)).
To further probe the monolithic HKUST-1 PDF, we analyzed the data using previously described non-negative matrix factorization (NMF) techniques 26 . Two components were used to describe the data (Figure 4 (d, i-l) & Figure S29). Comparison of these maps ( Figure S30) Where Δρ = contrast, F(Q,r) = scattering form factor, V(r) is the particle volume, N is the total number of scattering particles, and P(r) is the probability of occurrence of a scatterer size r, and were modelled with a log-normal distribution using a spheroidal form factor. The Guinier regions of these materials, the region which describes a particles radius of gyration (Rg), were outside the range probed and could S42 not be fit. The power law slope of I(Q) at low scattering angle for each sample was therefore fitted using the simplified relationship: Where, intensity (I) can be fitted by a prefactor (B) across the scattering range (Q) and a power law slope (P) and indicated surface fractals for both powder and monolithic samples (ca. 3.6). Size distribution models were fitted with a log normal size distribution and spheroidal form factors.
Powder samples consistently indicated two broad size distributions (mean diameters = ca. 24 nm and ca. 92 nm), whereas monolithic samples consisted of a single size distribution (diameters = ca. 20 nm) (Figure S32).

Raman microscopy experimental
Raman microscopy was performed using a Horiba LabRAM HR Raman microscope equipped with an automated xyz stage (Märzhäuser). To simultaneously scan a range of Raman shifts, a 600 lines mm -1 rotatable diffraction grating along a path length of 800 mm was employed. Spectra were acquired using a Synapse CCD detector (1024 pixels) thermoelectrically cooled to −60 °C. Before spectra collection, the instrument was calibrated using the zero-order line and a standard Si (100) reference band at 520.7 cm -1 .
For single point measurements, spectra were acquired using a 532 nm laser (at a power of 0.2 mW), a 100× objective, and a confocal pinhole of 200 µm, over the range 50-4000 cm -1 (3 spectral windows) with an acquisition time of 60 seconds and 2 accumulations to improve the signal to noise ratio and remove the spikes due to cosmic rays. The spectral resolution in this configuration is better than 1.9 cm -1 . Spectra were averaged by collection from a minimum of two locations with similar colouration from optical images.
An interfacial spectroscopic map was acquired using a 532 nm laser (at a power of 0.       Figure S37. Raman spectra intensities of peak intensities associated with dark blue (805-845 cm -1 ) and light blue (765-805 cm -1 ) over the mapped section for monoHKUST-1.

Isosteric heat of adsorption (Qst) calculations.
A virial-type expression of the form below was used to fit the combined isotherm data for all the compounds at 273 and 298 K, where P is the pressure described in Pa, N is the adsorbed amount in mmol/g, T is the temperature in K, ai and bi are virial coefficients, and m and n are the number of coefficients used to describe the isotherms. Qst is the coverage-dependent enthalpy of adsorption and R is the universal gas constant. All the related fitting curves are shown in Figure   ) Pressure (Bar) Figure S41. Absolute (Nabs) H2 adsorption isotherms at 77 K for monoHKUST-1 compared to existing benchmark MOFs for high pressure H2 adsorption 29,30 . Note: Volumetric performance for IRMOF-20, DUT-23-Co, MOF-177, SNU-70, UMCM-9, NU-100 and MOF-5 are based on a packing density of 0.2 g cm - 3 29 . A packing density of 0.447 g cm -3 was used to estimate the volumetric H2 uptake for MIL-101 30,31 . HKUST-1 powder and monoHKUST-1 volumetric performances were based upon an experimental packing density (0.2979 g cm -3 ) and envelope density (1.07 g cm -3 ), respectively.     HKUST-1 powder and monoHKUST-1 volumetric performances were based upon an experimental packing density (0.2979 g cm -3 ) and envelope density (1.07 g cm -3 ), respectively. Figure S45. Comparison of the gravimetric total pore volume vs. gravimetric BET area of monoHKUST-1 compared to existing benchmark MOFs 12,27-31 . Figure S46. Comparison of the volumetric BET area vs. volumetric H2 adsorption capacity at 100 bar and 77 K for monoHKUST-1 compared to existing benchmark MOFs. Note: Volumetric performance for IRMOF-20, DUT-23-Co, MOF-177, SNU-70, UMCM-9, NU-100 and MOF-5 are based on a packing density of 0.2 g cm - 3 29 . A packing density of 0.447 g cm -3 was used to estimate the volumetric H2 uptake for MIL-101 30,31 . Ni(dobdc) performance was calculated using the previously reported packing density of 0.366 g cm -3 28 . HKUST-1 powder and monoHKUST-1 volumetric performances were based upon an experimental packing density (0.2979 g cm -3 ) and envelope density (1.07 g cm -3 ), respectively. Figure S47. Comparison of the mechanical degradation of MOF vs. compaction density for monoHKUST-1 compared to existing benchmark MOFs. The y-axis corresponds to the ratio between the maximum excess adsorption at 77 K for a MOF compacted to a specific density divided by the value for the initial value measured for the powder. The x-axis corresponds to the density of the compacted MOF divided by its crystal density. Note: Values for MOF-5, MOF-177, SNU-70 and NU-100 have been adapted from Ref. 29 . nanoMOF-5 has been adapted from Ref. 41 . monoHKUST-1 volumetric performances were based upon an experimental envelope density of 1.07 g cm -3 .

S59
Monolith Stability Testing Figure S48. Linear 77 K N2 adsorption isotherm for pristine monoHKUST-1 and a monoHKUST-1 sample stored at room temperature in a desiccator for 18 months. Closed circles represent adsorption whilst open circles represent desorption.