Simultaneous Estimation of Gas Adsorption Equilibria and Kinetics of Individual Shaped Adsorbents

Shaped adsorbents (e.g., pellets, extrudates) are typically employed in several gas separation and sensing applications. The performance of these adsorbents is dictated by two key factors, their adsorption equilibrium capacity and kinetics. Often, adsorption equilibrium and textural properties are reported for materials. Adsorption kinetics are seldom presented due to the challenges associated with measuring them. The overarching goal of this work is to develop an approach to characterize the adsorption properties of individual shaped adsorbents with less than 100 mg of material. To this aim, we have developed an experimental dynamic sorption setup and complemented it with mathematical models, to describe the mass transport in the system. We embed these models into a derivative-free optimizer to predict model parameters for adsorption equilibrium and kinetics. We evaluate and independently validate the performance of our approach on three adsorbents that exhibit differences in their chemistry, synthesis, formulation, and textural properties. Further, we test the robustness of our mathematical framework using a digital twin. We show that the framework can rapidly (i.e., in a few hours) and quantitatively characterize adsorption properties at a milligram scale, making it suitable for the screening of novel porous materials.

S1 Textural Characterization of the Materials Micromeritics AutoPore IV Series Mercury Porosimeter to characterize the pore size distribution in the macropore region. The data from panels (a) and (b) is analyzed to obtain the textural characteristics, given in Table 1 of the main text, of the three materials. (c) X-ray diffraction patterns obtained using a Malvern Panalytical XRD X'pert PRO to characterize the crystal structure of the adsorbent materials.
In Figure S1a, we illustrate the raw data obtained from the N 2 sorption isotherms at 77 K for all the three materials. We can make several observations and gain insights into the nature of the pore structure within the three materials studied in this work. First, AC and 13X conform to the Type I isotherm classification. For these materials, we can observe a steep uptake of N 2 at low partial pressures followed by a plateau, characteristic of microporous materials. 1 Second, we can observe that the micropore filling occurs at a lower relative pressure (P/P 0 ) in 13X when compared to AC which corresponds to a decrease in micropore width and an increase in adsorption energy. Additionally, the range of relative pressure required to reach a plateau is smaller in 13X when compared to AC, indicating a limited range of micropore sizes as seen in Figure 2 of the main text. Finally, BN conforms to a Type IV isotherm, as observed elsewhere, 2 characteristic of a mesoporous material. This can also be observed in Figure 2 of the main text.
In Figure S1b, we illustrate the raw data obtained from Mercury Intrusion Porosimetry (MIP) for all the three materials. We can make two observations. First, for AC and BN, the pressure at which intraparticle filling is observed is the lowest among the three materials. This is indicative of a large range of macropore sizes. Second, 13X shows a sharp uptake at a higher pressure when compared to AC, indicating a smaller range of macropore sizes. These observations are also visualized in Figure 2 of the main text.
In Figure S1c, we illustrate the X-Ray diffraction (XRD) patterns for all the three materials. We can make two observations. First, both AC and BN show broad peaks in the ranges 20°to 30°and 40°to 45°, indicative of an amorphous structure. For BN these peaks are characteristic of turbostratic/amorphous BN, 3 which is the form expected to be obtained using the synthesis method we employed in the work. 4 Second, 13X exhibits sharp peaks, indicative of its crystalline structure.

S2.1 Experimental Setup
Photographs of the setup are shown in Figure S2. The details of the corresponding components used to assemble the experimental setup is described in Section 3.1 of the main text.

S2.2 Degassing Protocol
For the volumetric and MIP experiments (Section 2.2 of the main text), we degassed the AC and BN ex situ in the degassing stations of an Autosorb iQ (Quantachrome Instruments, United States) using the following protocol: (1) 1 h at 323 K; (2)  any, in the obtained MS signal due to changes to the total flow rate of the gas mixture is captured. We also performed a repeat of the MS calibration after 83 days of the first calibration. We used an in-built function fit with linear interpolation in MATLAB R2020a (The Mathworks Inc., United States) to obtain the calibration curve for the two repeats performed on two different days at different flow rates.
The calibration curve obtained from the two repeats is visualized in Figure S3. Two key observations can be made. First, the calibration curve for both the repeats are nonlinear. Second, there is a significant drift in the calibration between the two repeats. Such a drift in the calibration will have implications on the actual gas phase composition and the shape of the desorption response obtained during a dynamic sorption experiment. For instance, due to the drift, at a signal intensity ratio I He = 0.50, the differences in the esimated gas composition can be as high as 5 %. The implication of this behavior is further stark when observing low CO 2 compositions (y CO 2 < 0.10) visualized in Figure S3b. At I He = 0.90, the obtained CO 2 compositions y CO 2 can vary between 0.02 to 0.04.
For materials that exhibit a linear or moderatly nonlinear isotherm, such a difference will not have a huge impact on the estimated adsorption properties. However, for materials that exhibit a nonlinear isotherm (like 13X), it is common to observe a long tail at low gas compositions due to the shape of the isotherm (see Figure 3 of the main text and Figure S8). Therefore, an inaccurate estimation of the gas composition would lead to an inaccurate estimation of both adsorption equilibria and kinetics. As stated in the main text, this is one of the possible reasons for the disagreement observed between the equilibrium behaviour predicted for 13X by the approach presented in this work and the volumetric approach. Therefore, it is necessary to repeat MS calibrations at frequent intervals. Additionally, when exploring such broad ranges of compositions, it might be beneficial to use a more sensitive sensor at low gas phase compositions in tandem with a MS.

S2.5 Calibration of Mass Flow Meter
The mass flow meter (MFM) is originally calibrated to a pure gas (in this case He). In order to F in is the sum of the flow rate obtained from the two MFCs. For each flow rate, we obtained gas mixtures of varying CO 2 composition y in , ranging from pure He to pure CO 2 . Finally, to obtain the actual mass flow rate F act , we fit the resulting flow rate measured by the UMFM F true to the flow rate measured by the MFM F MFM , using a 3 rd order polynomial, and CO 2 gas composition y in , using a 2 nd order polynomial. We used an in-built function fit in MATLAB R2020a (The Mathworks Inc., United States) to obtain the following polynomial expression that provides the actual mass flow rate F act of the mixture

S3.2 Simultaneous Estimation of Adsorption Equilibrium and Kinetics
Prior to evaluating the objective function given in eq S2, we preprocess the experimental response from Section 4.2 of the main text and the computational response from Section 5.2 of the main text as follows: • We downsample the experimental or computational response with a significantly longer duration when compared to the ones with the shortest duration to ensure that the objective function is not dominated by the longer responses • We concatenate all the experimental or computational responses considered for the parameter estimation, resulting in a total number of N t data points • We scale the experimental responses from Section 4.2 of the main text and computational responses from Section 5.2 of the main text so as to have comparable magnitudes (each composition is subtracted by the minimum composition obtained from the concatenated data and divided by the maximum difference of the composition and the minimum composition obtained from the concatenated data) • We scale the experimental responses from Section S5.4 so as to have comparable magnitudes at high compositions and low compositions of the responses (divide the compositions into two groups at the median of the concatenated data and scale the individual groups using the methodology presented in the previous point) and subsequently minimize the sum of the objective functions obtained from eq S2 for the two groups

S4 Maximum Likelihood Estimation and Uncertainty Calculation
Based on certain underlying assumptions described elsewhere, we can show that for maximum likelihood estimation (MLE), 5 the parameter vector θ that maximizes the probability of observing a given experimental data set, is the minimizer of the function where N t is the total number of experimentally determined data points,ŷ k is the experimentally determined value of the dependent variable, and y k (θ) is the corresponding model prediction, for a model parameter vector θ.
We rate the quality of the local minimizer θ * using two approaches. First, we visually inspect the fit between the experimental responses and the model predictions. Second, we compute the confidence regions to assess how well a minimizer θ * was determined. To this aim, we estimate the positive semidefinite parameter covariance matrix V θ after approximating the measurement error covariance matrix and computing the sensitivities of the model predictions with respect to the parameters. 5,6 Subsequently, we can obtain the hyperellipsoidal confidence region in the N p -dimensional parameter space as follows is the inverse of the chi-squared cumulative distribution function with N p degrees of freedom at the probability η. Finally, we take the bounding box of the hyperellipsoidal confidence region and we report the confidence intervals for each parameter.
To visualize the confidence regions in Section 4 of the main text, we sample 1000 set of isotherm parameters from the confidence intervals for each material using a latin hypercube sampling technique (LHS function of the SMT package in Python 3.8.5 7 ). Upon evaluating the amount of gas adsorbed at a given temperature and partial pressure for the 1000 sets of isotherm parameters, using eq 2 of the main text, we find the corresponding minimum and maximum. For all the results presented in this work, the confidence region for a given material at a given temperature corresponds to the area bounded between the minimum and maximum amount of gas adsorbed over the entire partial pressure.

S5.1 Blank Volume Experiments
We performed dedicated experiments for the purpose of characterizing the blank volume of the experimental setup. In all these experiments, we removed the adsorbent, saturated the setup with pure CO 2 , and finally switched to pure inert gas (helium) using a switching valve. As explained in Section 3.1.1 of the main text, we have divided the entire blank volume of the setup into two segments. To characterize these segments, we performed two distinct sets of experiments. To characterize Segment I, we used the same setup configuration as shown in Figure 1 of the main text, and we performed the experiments without the adsorbent. To characterize Segment II, we removed all the components of Segment I and fed the gas directly to the tee connector that splits the gas into two streams, one to the MS and the other to the vent. To elucidate the impact of flow rate on the blank response, we performed these experiments at four different flow rates. We performed experiments to characterize Segment I at a gas flow rate F in = [0.2 0.3 0.5 1.1] T cm 3 s −1 and to characterize Segment II at a gas flow rate F in = [0.2 0.3 0.8 1.1] T cm 3 s −1 at room temperature and pressure.

S5.2 Blank Volume Model
Here, we formulate a generic mathematical model for the blank volume (Segment I and Segment II ) in the experimental setup shown in Figure 1 of the main text (blue and green) by modeling individual components of the setup that are lumped into sections that exhibit a given flow behavior. [8][9][10][11] To elucidate the flow behavior, we conducted a few preliminary blank experiments based on the approach described in Section S5.1. The responses obtained from these experiments exhibited either an exponential decay response (see Figure S4a) or a double exponential decay response (see Figure S4b,c). The former response is characteristic of a plug flow behavior with axial dispersion. The latter response is characteristic of two volumes in series, namely one that exhibits a plug flow behavior with axial dispersion (for the time shift) and the other that exhibits a diffusive mixing behavior (for the double exponential). We have modeled these two segments by drawing inspiration from compartment models used in chemical reaction engineering, 12 that facilitates describing plug flow, mixing, diffusion and stagnation behavior of fluids. A schematic that visualizes the volumes in the two segments as a combination of chemical reactors is shown in Figure 1 of the main text. Note that the model equations presented in this section are written for the adsorbing gas in a binary gas mixture.
Segment I: The entire segment, shown in blue in Figure 1 of the main text, has the same gas flow rate F (t) as the outlet of the CSTR that models the gas adsorption in the adsorbent. We described the gas flow in the volume that exhibits a plug flow behavior with axial dispersion using a series of N S CSTRs, with a total volume V S [m 3 ]. Note that the higher the number of CSTRs N S , the lower the axial dispersion. The mass balance around a CSTR p is given as is the concentration of the gas at the outlet of CSTR p and F (t) is the volumetric mixture gas flow rate at the outlet of the CSTR that models the gas sorption in the adsorbent, obtained by solving the system of equations given by eq 8 of the main text.
We described the gas flow in the volume that exhibits a diffusive mixing behavior as two CSTRs in parallel, characterized by a mixing volume V M [m 3 ] with a volumetric mixture gas flow rate . The mass balance around these two reactors is given by where Finally, we combine the gas from the mixing and the diffusive volume, which is the outlet from Segment I, with a gas concentration c S [mol m −3 ] defined as where c M and c D are obtained from eq S5. Note that we assumed the pressure and temperature of the gas from the two volumes to be identical. The initial conditions for eqs S4 and S5 are is the temperature and y 0 [-] is the initial gas phase mole fraction in the blank volume of the setup.
Segment II: The entire segment, shown in green in Figure 1 of the main text, has a gas flow rate and it is dictated by the vacuum pump of the MS. We described the gas flow using a series of N MS CSTRs, with a total volume V MS [m 3 ]. The mass balance around a CSTR q is given as volumetric mixture gas flow rate through the MS, which is independent of the volumetric mixture gas flow rate F (t) in Segment I. The initial condition for eq S8 is where P [Pa] is the pressure, T [K] is the temperature and y 0 [-] is the initial gas phase mole fraction in the blank volume of the MS.

S5.3 Estimating Blank Volume Model Parameters
We obtained the blank volume model parameters by following a two-step approach. First, we obtain the model parameters for Segment II, i.e., for the segment with the MS, by using the corresponding model embedded into a parameter estimator. Second, using the estimated model parameters for Segment II, we obtain the model parameters for Segment I, i.e., for the segment with the adsorption cell, the MFM, and the gas lines, by using a combined model comprised of Segment I and Segment II embedded into a parameter estimator.
Prior to performing parameter estimation, we preprocessed the raw experimental data, obtained using the protocol described in Section S5.1, through the steps outlined in Section S3. We obtained the minimizer θ * , using a maximum likelihood estimator described in Section S4, for the models corresponding to the blank volume. We solved the minimization problem, given by eq S2, to estimate the model parameters using an evolutionary optimization algorithm, i.e., a standard and elitist genetic algorithm (geneticalgorithm2 46 (v6.2.4) in Python 3.8.5). In eq S2, for all the cases,ŷ k is defined as the kth point of the experimentally measured outlet CO 2 gas composition obtained from the MS and  Table S2.

S5.4 Characterization of Blank Volume
We performed several experiments, using the protocol described in Section S5.1, to characterize the blank volume in the two segments (Segment I and Segment II ) of the experimental setup. When characterizing Segment I of the experimental setup, we considered two sub-cases. These two subcases correspond to configurations with (referred to as Segment I w/ Ball, applied for AC and 13X) or without (referred to as Segment I w/o ball, applied for BN) the stainless steel ball placed inside the adsorption cell to reduce the blank volume and to increase the heat dissipation. Subsequently, we estimated the blank volume model parameters for the corresponding experimental configuration using the methodology presented in Section S5.3. indicate the corresponding simulated response generated using the model described in Section S5.2 with model parameters given in Table S1. The stainless steel ball has been used in some experiments to reduce the blank volume of the setup and to help dissipate heat generated during gas adsorption.
The time-resolved experimental and the simulated blank responses y CO 2 for Segment I (both w/o Ball and w/ Ball) and Segment II are shown in Figure S4. The estimated blank volume model parameters are provided in Table S1. We can make three observations based on the experimental and the simulated outcomes. First, we can observe differences in the experimental response (markers) obtained from Segment I and Segment II. The response from the former, shown in Figure S4b,c, is flow dependent and is characteristic of volumes that exhibit a plug flow behavior (for the time shift) coupled with a diffusive mixing behavior (for the double exponential). The response from the latter, shown in Figure S4a, is flow independent and is characteristic of plug flow behavior with axial dispersion. 12 Second, in all the three cases, the simulated blank response (solid lines) agrees well with the experi-

S6.1 Repeatability of Experiments
We carried out dynamic sorption experiments on the three materials using the protocol described in

S6.2 Flow Dependency on Adsorption Behavior
To enable accurate estimation of both adsorption equilibria and kinetics from desorption responses, one should perform experiments at different inert gas flow rates. Low flow rate experiments enable accessing the equilibrium regime, while high flow rate experiments enable accessing the kinetic regime of adsorption. We can check the regimes being accessed by a given flow rate on a given material by visualizing the desorption response in the so-called Ft plot, 13 shown in Figure S6. In a typical Ft plot, if the curves overlap at different flow rates, we cannot extract any kinetic information as the adsorbent-gas system is under equilibrium. In Figure S6, due to the presence of blank volumes in the system, we cannot directly plot the curves at different flow rates. Therefore, once we have estimated the adsorption equilibria and kinetic parameters, we feed them into the gas adsorption model described in Gas Adsorption in the Adsorbent in Section 3.2 of the main text. Subsequently, we plot the desorption responses at the outlet of the adsorption cell. We can observe that for all the three materials, the Ft curves at the two different flow rates explored in this work do not overlap. Therefore, we can be confident that both the equilibrium and kinetic parameters are satisfactorily extracted and well-determined.  Table 2 of the main text.

S6.3 Objective Function
To understand the robustness of the optimization technique used in this work, for all the three materials we repeated the parameter estimation five times, using the methodology discussed in Section 3.3.2 of the main text. The objective function values J(θ), given by eq S2, obtained at the local minimizer θ * for the five repeats of the parameter estimation for all the three materials with the experimental and the digital twin response is provided in Figure S7. For the experimental responses, obtaining the actual value of the objective function with the true equilibria and the kinetics of a given material is not possible. However, for the digital twin responses, we can compare the objective function values J(θ) obtained through the five parameter estimation repeats with the true objective function J true . The latter can be obtained by evaluating the objective function, given by eq S2, with adsorption equilibrium and kinetic parameters, reported in Table S3, used to generate the computational responses. We can make two observations from Figure S7. First, the differences in the objective function values for the different repeats of a given material is very small. Second, for the computational case, we obtain higher values of objective function when compared to the true objective function J true . This leads to a situation in which the optimizer is unable to reach the true minima of the problem, and therefore prevents from obtaining the true adsorption equilibrium and kinetic parameters.  Table S3 is also additionally provided (J true ). Table S1: Parameters for the blank volume models, described in Section S5.2 and shown in Figure S4, of the mass spectrometer (MS) and the experimental setup. The blank volume for the experimental setup was characterized in two configurations: one with and one without a stainless steel ball in the adsorption column.

Parameter
Unit Value

S7.2 Computationally Generated and Fitted Desorption Responses
The computationally generated time-resolved CO 2 gas composition y CO 2 for the 12 in silico experiments performed at four initial gas phase composition, two inert gas flow rates, and three temperatures (using the digital twin described in Section 5.2 of the main text), along with the corresponding model fit, with the lowest objective function, for all the materials is provided in Figure S8. We can make several observations based on both the in silico experimental and the simulation outcome that are common for all three materials. First, similar to the experimental studies shown in Section 4.2 of the main text, as expected, the in silico experiment performed at lower temperature takes longer to reach a given gas phase composition at a given flow rate than an experiment performed at higher temperature, indicating a higher adsorption capacity at a lower temperature. Second, we can observe a significant time lag between the blank response (light gray) and the desorption response with the adsorbent. Finally, we can see that the overall fit of the simulated response is in good agreement with the in silico experimental response for all the three materials throughout the entire composition range. Note that we can observe minor deviations for some flow and temperature conditions in all the three materials. These minor deviations lead to the observed discrepancies in the estimated isotherms when compared to the TRUE isotherms, shown in Figure 6 of the main text. y obtained using the digital twin described in Section 5.2 of the main text using the parameters given in Table S3. The solid curves indicate the corresponding simulated response generated using the model described in Section 3.2 of the main text with the parameters corresponding to the parameter estimation repeat with the lowest objective function. The gray line in each panel represents the blank response of the setup at the given flow rate and initial gas phase composition.

S7.3 Sensitivity to Model Inputs
In this second study to evaluate the accuracy of the proposed methodology, we want to understand the sensitivity of the model inputs, e.g., porosity, blank volume model, to accurately estimate the adsorption equilibria and kinetics of the materials. To this aim, we perform a full curve fit of the computationally generated responses for AC, using the simulation framework described in Section 3.3 of the main text.
However, we provide model inputs that are different from the true parameter values used to generate the response. For example, to understand the sensitivity of the porosity T on the accuracy of the estimates of the adsorption behavior, we generate the computational responses with the digital twin using a porosity T = 0.61, but we use a porosity of either T = 0.35 or T = 0.90 (see Figure S9) for the parameter estimation routine. The former is analogous to obtaining an experimental response from a material that exhibits a given porosity, and the latter is analogous to an incorrect experimental characterization of the porosity of the material and subsequently using it for the full curve fit in the simulation framework described in Section 3.3 of the main text.
Here  Figure S10. We subsequently used the corresponding blank volume model parameters as an input for the parameter estimation routine.
The CO 2 isotherms for AC at three different temperatures as a function of its partial pressure  Figure S4). The blank volume model response for these three cases is shown in Figure S10.
timation of the porosity ( T = 0.35). We can attribute this to a significantly higher gas phase volume for T = 0.90 (V g = 9V s ), given by eq 4 of the main text, for the same volume of solid when compared to T = 0.61 (V g = 1.6V s ) and T = 0.35 (V g = 0.5V s ). Since the total amount of gas in the system is fixed for a given desorption response, an increase in the volume of the gas phase translates to a reduction of gas in the adsorbent. This leads to a reduction in the estimated adsorption capacity at that condition. Third, an overestimation of the mass of the adsorbent (m ads = 65.63 mg), leads to an underestimation of the adsorption capacity and vice versa. Since the total amount of gas in the system is fixed for a given desorption response, an overestimation of the mass of adsorbent translates to the same quantity of gas being adsorbed on a larger quantity of adsorbent when compared to the reference case. This leads to a lower estimated adsorption capacity on a unit mass basis. We can attribute this observation to a poor agreement of the blank volume model responses for TIS and TIS + D/M with the experimental response, as shown in Figure S10. In more detail, depending on the flow rate of the gas and the blank volume model chosen, we either underpredict or overpredict the blank volume of the setup. Further, due to mass balance constraints, this will lead to a corresponding overprediction or underprediction of the estimated adsorption capacity, respectively.  Table S1.
To study the impact of the blank volume model, we have identified three candidate blank volume representations, namely, the entire blank volume (Segment I and Segment II ), described in Blank Volume in Section S5.2, to exhibit either a plug flow behavior with axial dispersion (tanks in series, TIS) or a combination of a plug flow behavior with axial dispersion and a diffusive mixing behavior (TIS + D/M). As a reference to generate the computational response, we have used a combination of a plug flow behavior with axial dispersion and a diffusive mixing behavior to describe the blank volume in Segment I and a plug flow behavior with axial dispersion to describe the blank volume in Segment II (TIS + D/M + MS). We reevaluated the corresponding model parameters by refitting the experimental response, using the aforementioned models. The experimental and model responses for all the three models are visualized in Figure S10. We can observe that the overall fit of the TIS + D/M + MS model is in good agreement with the experimental response for all the flow rates. The other two models, depending on the flow rates, the time required to reach a given composition is either under-or