A simulation framework for extracting intrinsic kinetics in transient experiments

Abstract

Transient experiments provide a unique vantage point for heterogeneous catalysis, where the kinetic properties of complex industrial materials can be precisely characterized in a highly controlled manner. Dynamic variation of catalyst surface states and the changing response of chemical reactions can bring great insight, but these methods require complex analysis. Temporal Analysis of Products (TAP) is one such method, used to measure kinetic properties by separating the intrinsic reaction on a catalytic surface from the mass transport in the reactor using precisely controlled reactant pulsing under low pressure conditions. However, calculating intrinsic kinetic quantities from the exit flux measured in TAP experiments requires careful data analysis and/or modeling. In this paper, we demonstrate a virtual TAP reactor model (VTAP) that connects the observed exit flux with the reactor concentration profile and catalyst surface state evolving as a function of time. Simple adsorption processes and complex catalytic reactions are modeled and discussed. As kinetic quantities and number of active sites are changed, the presentation of distinct rate/concentration ‘fingerprints’ emerge that form the basis of benchmarking catalyst behavior. These reaction simulations are used to interpret experimental pulse response data collected on both simple, Pt/SiO₂, and complex, MoC_x/ZSM5, catalysts. The strategies for interpreting the reactor exit flux data to extract intrinsic and transient kinetics quantities using the VTAP model are discussed. Transport and reaction simulations supported by the VTAP model framework provide clear visualization of the unique reactor physics and catalyst dynamics, laying the groundwork for designing more informative experiments that advance industrial catalysis.

---

Introduction

Dynamic experiments offer a way to deconstruct complex behavior by examining the system response to stimulus. Nevertheless, the quantitative analysis of dynamic response brings its own set of challenges. Temporal Analysis of Products (TAP) is a powerful methodology that measures the intrinsic kinetics of surface reactions on complex, industrial catalysts.^1–5 The method is defined by the use of a small (10 nmol) gas pulse into an evacuated packed bed (10⁻⁸ torr) which grossly simplifies the reactor transport and enables the kinetic response to be clearly observed. This is the basic principle of the experiment but the practical implementation of a quantitative kinetic analysis, is far from simple. Moreover, the physical behavior of this diffusion reactor and the way the concentration gradients evolve in time are often misconstrued, being overshadowed by the workings of a plug flow reactor that chemical engineers are so accustomed to. Interpreting the reactor exit flux data to extract intrinsic and transient kinetics quantities requires careful application of data analysis methods.^1,2,4,6–10

In the past several decades, the TAP data analysis method mainly sought to extract intrinsic kinetic information based on established gas diffusion theory¹ and approximations leading to one-zone, multi-zone^11–14 and thin-zone TAP reactor models.^2,15 Advanced methods developed more recently include the Y-procedure (physics-based inverse diffusion)^4,9 and G-procedure (probabilistic)⁵ for extracting the time-dependent concentration and reaction rate in the catalyst zone. The Y- and G-procedures do not require an input of a kinetic model based on confirmed or hypothesized reaction mechanism and are thus also referred to as model-free data analysis methods. This model-free feature is the advantage of the Y- and G-procedures in enabling fast and easy calculation of catalytic reaction rate and concentration in catalyst zone.

Recently, microkinetic modeling has been carried out to analyze TAP experiments with many numerical methods developed for simulation and parameter estimation, each with their own advantages.^6,16–18 For example, TAPFIT is a Fortran-95 code that both simulates TAP pulses and regresses parameters.¹⁸ It forward solves using either a method-of-lines discretization via LSODE (ODEPACK) or a transfer-matrix formulation in the frequency domain for pseudo-monomolecular cases. Regression is Levenberg–Marquardt via ODRPACK. TAPsolver is an open-source Python package built on FEniCS and dolphin-adjoint.¹⁶ It treats TAP as a PDE problem, uses finite elements for the forward model, and algorithmic differentiation to address sensitivities. Those gradients drive PDE-constrained optimization for kinetic parameter estimation, with options to enforce thermodynamic-consistency constraints and sensitivity analysis. SimTAP is a purpose-built forward simulator (MATLAB) aimed at rapidly reproducing full non-steady-state TAP pulse sequences for complex, multi-step mechanisms, used in practice to validate or falsify competing kinetic models by matching entire transient experiments.⁶ The method emphasizes speed and forward simulation rather than integrated inverse modeling. Among these tools, there is no unified demonstration of calculating detailed gas-phase concentration profiles, catalyst surface state, rates of reaction, and their evolution in time within a single model.^6,16,18

In this paper, we introduce a vectorized 1-D virtual TAP model (VTAP) that represents diffusivity, site density, kinetic constants, gas concentrations, and surface coverages as position-dependent arrays on a single axial grid. This removes zone bookkeeping, lets properties vary smoothly, and outputs both exit flux and internal spatiotemporal fields. We demonstrate microkinetic modeling of TAP pulse response experiment for adsorption processes and catalytic reactions using these internal fields to extract detailed information on gas-phase concentration, catalyst surface state, and reaction. Further, VTAP handles richer cases such as multiple site types and nonlinear networks. As the magnitude of these reaction rate constants and the number of active sites change, VTAP formalizes the rate/concentration dependences that emerge serve as a kinetic fingerprint to identify and distinguish different catalysts. With these simulation results as a guide, experimentally derived rate/concentration features are interpreted for reactions on both simple, Pt/SiO₂, and complex MoC_x/ZSM-5, catalysts. Potential applications and strategies of VTAP are discussed.

Methods

The TAP reactor is designed to distinguish the time-dependent effects of mass transport from intrinsic chemical reaction on the catalyst surface. A schematic figure of the TAP reactor is demonstrated in Fig. 1A, the zones and key positions of the packed bed to be discussed in this paper are depicted in Fig. 1B. The pulse valves on top of the reactor bed precisely control the pulse amount of gas injected into the reactor in the range of several nanomoles on a time scale of microseconds (10⁻⁶ s). The reactor is connected to a vacuum chamber with a high throughput pumping system that creates very low pressure in the reactor in range of 10⁻⁸ Torr. The combination of the small pulse amount, a narrow pulse width, and very low pressure in the system guarantee the mass transport of gas in the reactor bed is in the Knudsen diffusion regime. The theory has been presented in detail and extensively discussed in literature.^2,4,5,9,16

Fig. 1 (A) Schematic figure of the TAP reactor. (B) Zones and positions in the TAP reactor bed, positions 1 through 6 are begin and middle of zones 1, 2 and 3, respectively. (C) Experimental result of exit flux of a calibration experiment of pulsing oxygen through TAP reactor filled with inert quartz. (D) Simulation result of concentration profile in the TAP reactor from the VTAP model.

The reactor model describes transport of the gas in the void volume of the TAP reactor and chemical reaction on a catalyst loaded in the middle of reactor bed. The mass transport can be described by Knudsen diffusion along the reactor axis due to the low pulse intensity (less than approximately 10 nmol) in TAP experiments. The equation of continuity for a non-reacting gas A in a packed bed reactor is given by

The initial condition for application of the above equation (eqn (1)) to a TAP pulse response experiment is given as C_A = 0 at 0 ≤ z ≤ L and t = 0, representing empty gas phase in the reactor bed. The pulse event is defined as filling the top 5% in length of the first inert zone (as shown in Fig. 1B) of the TAP reactor with concentration C_A of gas A, which corresponds to approximately 2.5% of whole reactor length. This method of setting initial condition has been used in TAPsolver simulation approach,¹⁶ in which the pulse event is defined as filling top 5% of whole reactor. TAPsolver supports constraints to ensure thermodynamic consistency and accurate descriptions of microkinetic models with exit flux as simulation output. Comparatively, the focus of VTAP is to simulate detailed transient states of gas-phase concentration and catalyst surface coverage, as well as the exit flux. The time needed for pulse injection is approximated to zero due to its shorter time scale of microseconds than the millisecond time resolution of a TAP experiments. The pulse size (N_p) in molar amount of gas injected into TAP reactor, can be calculated by

N_p = C_A × 5% L × πr². (2)

Boundary conditions are given by

C_A = 0 at z = L. (4)

The effect of the position of the catalyst on the catalytic reaction is given byThe gas flux (nmol cm⁻² s⁻¹) observed at exit of reactor is given bywhere ε_b is the fractional voidage of the packed reactor, C_A is concentration of gas A (nmol cm⁻³), z is the axial coordinate (cm), L is the length of the TAP reactor, r is the inner radius of the TAP reactor, ρ_s is the active site density, r_A reaction rate of gas A per active site of catalyst surface. The site density (ρ_s) is defined as the number of active sites referring to the number of molecules in same volume of an ideal gas at standard temperature and pressure (STP). Therefore, a standard site density (ρ_s = 1) corresponds to 4.46 × 10⁴ nmol cm⁻³.

For the adsorption of gas A on the catalyst surface,

A(g) + * ↔ A* (7)

the rate r_A is given by

r_A = k_aC_Aθ_* − k_dθ_A* (8)

where k_a and k_d are kinetic constants of adsorption and desorption, θ_A* and θ_* are surface coverage of A and empty site, respectively. We do not define the adsorption process as a reversible or irreversible reaction as an input of the model; the reversibility is determined under simulation by the kinetic constants and the surface coverages.

For reaction with a complete catalytic cycle, we employ the mechanism of a three-step reaction including adsorption of reactant, conversion and desorption of product as an example.

A(g) + * ↔ A* ↔ B* ↔ B(g) + *. (9)

The reaction rates of these steps are given by

r₁ = k_1fC_Aθ_* − k_1bθ_A* (10)

r₂ = k_2fθ_A* − k_2bθ_B* (11)

r₃ = k_3fC_B − k_3bC_Bθ_*. (12)

While keeping the theoretical foundation the same as the widely accepted TAP-2 model of Gleaves et al.,¹ we have several modifications in technical application. We vectorize the parameters i.e. diffusivity, kinetic constant, concentration, and surface coverage over the axial coordinate in the simulation code. This enables flexibility in defining zones and gradients inside a zone, and increases extendibility to more complex scenarios. Benefiting from parameter vectorization, the model does not split the reactor into zones like the multi-zone reactor model,^11–14 it is a continuous model where the catalyst can exist at any position without need to add zones and boundary conditions. Since this model does not use an approximation of uniform catalyst zone concentrations of thin-zone TAP reactor model,¹⁵ it is flexible to different experimental setups and opens the opportunity to studying effects of catalyst loading pattern, loading amount, catalyst gradients, etc. VTAP solves the same 1-D diffusion–reaction equations, but treats the whole bed as one continuous field so that properties can change smoothly along z. As such, one can resolve spatiotemporal gas-phase concentrations and evolving surface coverages. This adds flexibility for future studies, for example mimicking the variation in diffusivity and rate constants caused by temperature gradient inside the catalyst zone(s) and inert zones. In the model results described here, the TAP reactor length is set to 4.0 cm, and the catalyst is loaded in the middle of the reactor with catalyst zone length of 0.2 cm. The fractional voidage is set to 0.4 throughout the TAP reactor, and the diffusivity of gas is set to 40 cm² s⁻¹.

Results and discussion

Mass transport in the TAP reactor

In the Knudsen diffusion regime, the concentration change in the reactor as a function of time is described using eqn (1) and the exit flux is described using eqn (6). In the case of a calibration experiment, where there is no chemical reaction since the reactor is only filled with inert quartz, the exit flux is solely determined by the diffusion of the gas pulsed in the inlet area of TAP reactor. Fig. 1C shows an example of exit flux measured by the mass spectrometer located close to the outlet of the reactor bed. Since this flux curve is not affected by a chemical reaction, it can serve as measuring stick to experiments with catalysts enabling chemical reactions under same condition; it is also known as standard diffusion curve (SDC). The comparison of the flux curves and analysis of chemical reactions will be discussed later in this paper.

The exit flux curve is not linearly dependent on the concentrations at the catalyst zone (zone 2 in Fig. 1B) and other positions in the reactor because of the effect of gas diffusion. It should be noted that the exit flux is the experimental measurable, but the concentration profile in the reactor bed cannot be directly measured and this information is useful for understanding the physics of the experiment. A model or data analysis method is required for predicting concentration and its evolution in the reactor bed. An example of the concentration profile with respect to the reactor axial coordinate and time of a pulse evolution is show in Fig. 1D. Pulsing a small amount of gas into the reactor at earlier time results in relatively high concentration of gas only in the beginning of zone 1, while there is no gas concentration in the later part of the reactor bed, especially zones 2 (catalyst zone) and 3. This feature can be seen from the red line in Fig. 1D, which is corresponding to an early time (t = 0.001 s). There is a short initial period where the pulse will not be in the Knudsen regime, but this duration is insignificant and does not impact the measurement. The gas concentration in zone 1 rapidly drops and the concentration profile becomes flatter due to the expansion of the gas throughout the reactor, as depicted as the blue arrows in Fig. 1D. An animation of the diffusion behavior, such as exit flux and concentration profile during observation of pulse experiment, is provided in the SI (animation 1. standard diffusion) for better visualization and understanding in different dimensions.

The factors that may affect the SDC are tested with the VTAP model developed in this paper. Fig. 2A shows that the diffusivity affects both peak time and peak height when other parameters such as reactor length, fractional voidage and pulse size are kept constant. Large diffusivity causes early peak time and high peak height (strong signal), while the area under flux curve remains unchanged. This trend agrees well with theory and previous observations using dimensionless parameters.¹ It should be noted that there is a delay in time for the exit flux signal after the pulse injection (t = 0), as depicted in the subplot in Fig. 2A. The delay is reasonable because it represents the time needed for the gas to diffuse from the inlet area through the reactor and to be measured by the mass spectrometer.

Fig. 2 (A) Effect of diffusivity on the flux curve along time (inset for logarithmic scale), reactor length L = 4.0 cm and fractional voidage = 0.4. (B) Concentration profiles in the TAP reactor for the case of D = 40.

The concentration profile is shown in Fig. 2B. The concentration is high toward the inlet in the beginning of time (t = 0.005 s) as shown by the profile sampled before the peak (solid red). Note that at this time, the pulse has started to emerge from the reactor exit, Fig. 2A. The concentration profile is much flatter both at and after flux peak time. The diffusivity does not change this feature but affects the speed of relaxation. Larger diffusivity causes faster relaxation in the concentration profile. This observation is consistent with flux curve in Fig. 2A where higher diffusivity has earlier peak time, which means that gas diffuses through reactors faster.

To evaluate the performance of the VTAP model, the diffusivities of gases are fit to the exit flux curves in calibration experiments where gases are pulsed into the TAP reactor filled with inert quartz. As shown in Fig. 3A, the model fitted SDC aligns well with flux curve of Ar in the calibration experiment. The fitted SDC of Ar has some deviation from the experimental flux curve which is likely related to temperature gradient in the reactor. Taking temperature gradient into account of VTAP simulation is exactly what it is designed for, however, the experimental measurement of the reactor temperature gradient is currently a challenge. There are other sources of deviation, such as weak physisorption of molecules on inert quartz and tortuosity of packed bed macropores within the inert quartz. The weak physisorption on inert quartz is common for molecules containing hydrogen. The effect of physisorption is more obvious for ammonia (NH₃), which has been considered as a “sticky” molecule,¹⁹ causing deviation of calibration flux curve from SDC (Fig. 3B).

Fig. 3 (A and B) Comparison of experimental and model flux curves of Ar and NH₃ at 400 °C. (C) Fitted diffusivity (cm² s⁻¹) of molecules at 700 °C plotted against square root of molar mass of molecule. (D) Fitted diffusivity (cm² s⁻¹) of Ar at different temperatures.

Fig. 3C and D show the fitted diffusivities plotted against the square root of molar mass and temperature, respectively. The standard deviations are all within 0.5% of fitted diffusivities, indicating the high precision of the model in describing the mass transport. The error bars in Fig. 3C and D are estimated by considering the maximum and minimum values in 20–40 pulses of each experiment, showing good reproducibility in comparison of different experiments. Large diffusivity is observed with a light-weighted molecule and high temperature. Fast diffusion of light-weight molecules can be rationalized by Graham's law of diffusion stating that the rate of diffusion of a gas is inversely proportional to the square root of its molecular weight.²⁰ The diffusivities in Fig. 3C and D are fitted without considering the factors causing deviations from the ideal Knudsen diffusion as described above. The error bars in Fig. 3C and D indicate minor variations among experimental measurements. Comparatively, the diffusivity is proportional to temperature. There are several possible reasons for the outlier in Fig. 3C. Fig. 3C is a collection of experiments carried out at different times and in different reactors. Better correlation can be obtained by measurement in a single experiment to mitigate the difference in experimental settings. For hydrogen and helium, the mass spectrometer detection at low atomic mass units (AMU) is challenging, sensitivities can be nonlinear in this region.

Adsorption process

By carefully designing and executing TAP experiments, researchers can gain insight into the fundamental interactions between adsorbates and catalytic materials, which is crucial for optimizing catalysts in various chemical processes. By systematically measuring adsorption during the titration experiment, researchers gain valuable characterization of evolving kinetics and mechanisms. In this section, we discuss the transient features of the adsorption processes in TAP experiments using the developed model. The effort in TAP data analysis is to extract the information about intrinsic and transient kinetics in the reactor reflected from the flux curve. Intrinsic and transient kinetic information can be used to measure the catalyst property and reaction mechanism and their changes during the reaction. Although the TAP reactor is well designed to quickly and precisely measure large amounts of data over varying kinetic states, analyzing these TAP experimental data pulse by pulse is challenging, since all exit flux curves are very close in shape and absolute value as shown in Fig. 1C.

The adsorption and desorption process of gas molecule A(g) to a catalyst surface is added to the VTAP model. The exit flux, concentration profile and catalyst surface state (the coverage of the adsorbate) for the case of large adsorption kinetic constant (k_a = 1 × 10³ s⁻¹) and small desorption kinetic constant (k_d = 1 × 10⁻⁴ s⁻¹) are shown in Fig. 4. An animated figure is provided in the SI (animation 2. adsorption process). Due to the consumption of gas during the adsorption process, the peak height and area of the flux curve are both lower than those of SDC (Fig. 4A), which is in agreement with literature.¹ The consumption of gas causes a kink on the concentration profile in the catalyst zone (x = 1.9–2.1 cm, in Fig. 4B). Consequently, the concentration profile in the reactor is lower than the corresponding SDC in calibration experiment. As depicted in Fig. 4C, the coverage increases during the adsorption process, and the desorption process does not occur due to an insignificant desorption rate constant. The coverage of adsorbed A* on the catalyst surface has a gradient along catalyst zone. It should be pointed out that our model does not assume uniformity in catalyst zone. Therefore, the model can predict higher surface coverage at the beginning of the catalyst zone, which is reasonable because of the higher gas-phase concentration in this region.

Fig. 4 (A) The exit flux curve in a TAP simulation with adsorption (k_a = 1 × 10³ s⁻¹, k_d = 1 × 10⁻⁴ s⁻¹) in comparison with the standard diffusion curve. (B) Concentration profile of adsorbate in the gas phase. (C) Coverage of the adsorbate on the catalyst surface along the catalyst zone. (D) Average adsorption rates per total and available catalyst site plotted against the average gas-phase concentration in catalyst zone.

Fig. 4D shows the average adsorption rate per total active site (green) and per empty active site (blue) plotted against the average concentration in the catalyst zone. The adsorption rate per total catalyst active site is the apparent adsorption rate that can be observed in kinetic experiment. The apparent adsorption rate linearly increases with the increasing gas concentration in the early stage of the pulse (ascending branch), the peak of the pulse (maximum concentration) coincides with the maximum rate, then as the gas concentration decreases (descending branch) the rate is observed to return to zero at a lesser value. This curvilinear relationship has been referred to as a ‘kinetic petal’ for the flower petal shape it resembles.^3,5 The rate-concentration (RC) petal is a characteristic of the transient behavior that arises due to the decrease of available active site during adsorption process. The effect of decrease in active sites can be validated by comparison with the intrinsic rate per available empty active site, also known as turn-over frequency (TOF). As shown in Fig. 4D, the intrinsic rate follows an exact linear correlation with gas-phase concentration of A, since the number of empty sites doesn't affect the rate on each site in the tested case. The shape of the RC petal is determined by the gas phase concentration, surface coverage and kinetic constants. In this model, the kinetic constants are fixed during the adsorption process. It should be noted that the kinetic constants could also change in real experiments during the adsorption process due to the lateral interaction between adsorbates at high coverage and changes in catalyst structure. The variation of kinetic constant and binding energy can also be examined by fitting experimental results using the VTAP model.

The effects of catalyst site density on the exit flux curve and the RC petal for the case of large adsorption (1 × 10³ s⁻¹) and small desorption (1 × 10⁻⁴ s⁻¹) kinetic constants are illustrated in Fig. 5A and B. The small desorption kinetic constant guarantees the irreversibility of the adsorption process. The site density (ρ_s) is defined as the number of active sites referring to the number of molecules in same volume of an ideal gas at standard temperature and pressure (STP). Therefore, a standard site density (ρ_s = 1) corresponds to 4.46 × 10⁴ nmol cm⁻³. This quantity is not directly used in the TAP experiment, but it is a convenient definition for modeling. In a realistic experiment, the site density can be affected by the catalyst loading amount and the intrinsic catalyst properties, such as surface area and active components. Directly using the site density can avoid unnecessary complications in modeling. Fig. 5A shows that the height of the exit flux curve decreases with increasing site density, corresponding to high conversion (uptake) caused by the large catalyst loading amount. The adsorption rate per active site in the RC petal (Fig. 5B) decreases with increasing site density. In the case of high site density (ρ_s = 1), although the rate per site is lower than (approximately half) the low site density (ρ_s = 0.1), the overall adsorption rate in the reactor is higher due to the ten times larger number of available sites. In comparison of the RC petal shapes, a lower site density leads to a dilated RC petal shape. This can be predicted by a low catalyst loading relative to the pulse size (amount of gas injected into the reactor) inducing a quicker drop in number of empty sites. Therefore, for an irreversible adsorption process, a narrow RC petal is the fingerprint of low pulse size to catalyst loading amount (P/C) ratio, and broadened RC petal is the fingerprint of high P/C ratio.

Fig. 5 Effect of (A and B) site density, (C and D) adsorption kinetic constant and (E and F) desorption rate constant on the exit flux curve and the RC petal. The directions of RC petals in this figure are all clockwise.

Fig. 5C shows the effect of adsorption kinetic constant (k_a) on the exit flux curve while keeping desorption kinetic constant unchanged (k_d = 1 × 10⁻⁴ s⁻¹). It is not surprising that the exit flux peak decreases in height and area with an increasing adsorption kinetic constant, which is consistent to the high uptake caused by the fast adsorption rate. The trend in reaction rate can be seen more clearly in the RC petal in Fig. 5D, where the slope of RC petal is proportional to the adsorption kinetic constant of the irreversible adsorption process.

Fig. 5E shows the effect of the desorption kinetic constant on the exit flux curve when keeping the adsorption kinetic constant unchanged (k_a = 1 × 10³ s⁻¹). The change of the desorption kinetic constant causes the flux curves to cross each other. The general trend still holds that faster desorption kinetics resist the adsorption process and causes lower conversion, as depicted as higher flux peak height and larger area under the flux curve. The RC petals in Fig. 5F are very interesting because they show more obvious fingerprints of intrinsic and transient kinetic phenomena. The broad RC petal shape (black line) represents the case of equal adsorption and desorption kinetic constants, i.e. K_eq = 1. In this case, adsorption and desorption are balanced and the RC curve is determined by the balance between gas-phase concentration and catalyst surface coverage. Here, the rate maximum precedes the concentration maximum, with sites turning over as the gas concentration increases. When the desorption kinetic constant decreases to 1 × 10² s⁻¹ (red line), meaning the adsorption kinetic constant is now ten times greater than desorption kinetic constant, the adsorption process is favorable, yet the overall process remains reversible. A reversible feature is clearly identified by a negative rate during period of decreasing concentration (descending branch). However, the manifestation of reversibility in a negative rate depends on the balance of rate constants, i.e., the process is reversible in the first case (black line) but the negative rate is obscured. The RC petal for the case of K_eq = 10 is round. For the lowest tested desorption kinetic constant (k_d = 1 × 10⁰ s⁻¹), i.e. highest equilibrium constant (K_eq = 1 × 10³ s⁻¹), the RC petal is narrow (green line) with only a positive rate and R/C slope, which is typical for an irreversible process. Interestingly, the blue line shows the transition between reversible and irreversible in the case of an intermediate desorption kinetic constant (1 × 10¹ s⁻¹) and equilibrium constant (K_eq = 1 × 10²), characterized as existence of both positive and negative rates in the RC petal.

Overall, the above analysis of Fig. 5 demonstrates that information pertaining to the conversion rate can be obtained from the exit flux curve and more detailed information regarding the intrinsic and transient kinetic behavior is obtained from RC petal analysis. An advanced method for TAP experimental data analysis is to generate RC petals using the Y-procedure (inverse diffusion) or G-procedure (probabilistic). Prediction of RC petals from experimental data can be a challenging task. As such, the model in this paper can be used as a rational guide to interpret experimental RC petals.

TAP catalyst titration experiment measures the number of active sites available for a specific reaction by incrementally adding reactant molecules until saturation is reached. We modeled a series of multiple pulses for the cases of irreversible (k_a = 1 × 10³ s⁻¹, k_d = 1 × 10⁻² s⁻¹) and reversible (k_a = 1 × 10³ s⁻¹, k_d = 1 × 10² s⁻¹) adsorption process with a site density of 1. To assist in understanding the multiple-pulse process, an animation of the exit flux, concentration profile, surface state and RC petal with respect to time are provided in SI (animation 4 and 5). It should be noted that a molecule with irreversible adsorption on the catalyst can precisely quantify the uptake of adsorbate on the catalyst and is preferred for titrating the number of active sites.²¹ For irreversible adsorption, the exit flux curve (Fig. 6A) of the first pulse features low peak height and area due to the high uptake amount. The peak height and peak time increase with the increasing pulse number and, as the catalyst surface becomes saturated, approaches the standard diffusion curve. This is consistent with the changing RC petals (Fig. 6B) where the first pulse has a faster adsorption rate in combination with low gas phase concentration. With subsequent pulsing, the apparent adsorption rate decreases, and the gas phase concentration increases with increasing pulse number, indicating lower uptake with extended pulsing.

Fig. 6 Modeling results of (A) exit flux and (B) RC petal of irreversible (k_a = 1 × 10³, k_d = 1 × 10⁻²) and (C) exit flux and (D) RC petal of reversible adsorption (k_a = 1 × 10³, k_d = 1 × 10²) in a multiple pulse TAP experiment. The catalyst site density is 1.

The exit flux returning to zero at end of a pulse indicates the gas in TAP reactor has completely diffused out and the reactor is empty before the next pulse. In the model, the pulse repetition rate is set to 0.5 seconds, which is long enough for the exit flux of the first pulse to return to zero. However, it is noticed that the repetition rate (0.5 seconds) is not enough for later pulses to reach zero flux because of the slow adsorption rate caused by low coverage of empty surface site. This causes residual gas in reactor to accumulate and affects the flux curve of the next pulse. The effect of the residual is obvious by the fifth pulse and the flux curve before the peak (blue dots in Fig. 6A) is higher than SDC (orange dots in Fig. 6A). The direct comparison of fifth flux curve with SDC will lead to a negative rate, which is not true because the effect of residual gas from fourth pulse is ignored in such direct comparison. The RC petal from VTAP modeling in Fig. 6B indicates that the rate is positive.

The effect of residual gas is more obvious in the case of reversible adsorption. The overall adsorption rate is slow because of reverse reaction, 0.5 seconds is not long enough to allow the exit flux of the first pulse to return zero (green dots in Fig. 6C), leaving residual gas in the reactor at beginning of the second pulse. Consequently, the subsequent exit flux and concentration start from values greater than zero (purple dots in Fig. 6C and D). The effect of residual gas reaches a limit cycle after several pulses. In generating RC petals with Y- or G-procedures, the difference in flux curve compared to SDC is used to calculate the reaction rate. Without considering the residual gas concentration in the reactor, these analysis methods will very likely produce negative reaction rates for the early part of the RC petal. While such model-free data analysis methods are efficient for most cases, this example demonstrates the benefit of using VTAP to consider the complexities that may arise in experiments.

Catalytic reaction

Hereto now, we have discussed many factors affecting the TAP experiment in the case of adsorption of a gas molecule. We extend the VTAP model to a catalytic reaction using the Langmuir–Hinshelwood (L–H) mechanism. The focus of this section is to examine the general features of the TAP experiment with reaction of a complete catalytic cycle, including adsorption of reactant, surface reaction and desorption of product. We avoid adding extra complexity to the scenario by using a simple three-step reaction mechanism (eqn (9)). To reduce the complications in mass transport, the diffusivity of reactant A and product B are set to same value (40 cm² s⁻¹) to guarantee that the differences in concentration and flux curve are solely caused by catalytic reaction.

Although the same diffusivity is assigned to reactant A and product B, their flux curves (Fig. 7A) are very different. The most obvious feature is that the product flux has a later peak time than the reactant. It should be noted that the consumption of reactant A does not proceed at the same pace as the production of B, although it is a simple three-step reaction without side products and competing parallel pathways. The two gas phase adsorption and desorption steps are separated by a surface reaction step and the disconnect is more obvious in the RC petal (Fig. 7B) where A's consumption rate and B's production rate are plotted against the concentration of reactant A in the catalyst zone. While the consumption rate of A has a direct correlation with the concentration of A, the production rate of B does not. It should be noted that in our model we only apply direct numerical analysis of the reaction rate based on eqn (10)–(12) for the three steps, we do not employ a lumped kinetic model with a quasi-equilibrium assumption or steady-state approximation. A lumped kinetic model with these approximations forces correlations between reaction steps and intermediates, losing the ability to demonstrate the feature induced by transient kinetics.

Fig. 7 (A) Flux curve, (B) RC petal and (C and D) concentration profiles of reactant and product in the three-step catalytic reaction (eqn (9)).

The difference in reactant and product flux curves can be rationalized by the analysis of concentration profiles in the reactor. For reactant A, the concentration profiles at sampled times (Fig. 7C) have the same feature that we have demonstrated above for adsorption process. The kink in the middle of the concentration profile curve is caused by the consumption of A. However, the concentration profile of product B (Fig. 7D) has dramatically different features. The concentration profile at 0.02 seconds (left subplot) shows a symmetrical maximum at the reactor middle caused by the formation of product B; product B is free to diffuse toward both the inlet and outlet of the reactor. The concentration profiles at 0.1 and 0.2 second (middle two subplots) show the accumulation of production B in the inert zone on the inlet side (zone 1) which is a dead end. The overall concentration decreases after the production rate of B slows as shown as in the concentration profile at 0.5 seconds (right subplot). The source of reactant A and product B are at two different locations but both species diffuse through the same reactor volume. In the simulation, the diffusivity of A and B were made equal (this would only be expected for a simple isomerization reaction) and the difference in exit flux arises from distinct starting conditions: the reactant is initiated with a delta function filling the top 2.5% of reactor length and the product is initiated by the concentration time-dependence of the reactant through the catalyst zone at the middle of the reactor. In addition to the gas-phase diffusion, the exit flux peak time is also affected by the variation of B production rate. Especially in the case of the slow B production rate, the peak time of B is expected to be determined by the surface reaction instead of diffusion.

The RC petal for the production rate of B has the highest correlation with surface coverage of B*, followed by the concentration of B, as shown in SI (Fig. S5). This is reasonable because the reaction rate in eqn (12) indicates that the production of B is first order to surface coverage B*. The correlation of coverage, rate and concentration shows that the high coverage of B* is the driving force of fast production rate, and then is the reason for the high B concentration in the catalyst zone. In the calculated case, the rate determining steps for reactant A consumption and product B are different, because the reaction and surface state are both far from equilibrium. The RC petal features will be vastly affected by the reaction mechanism, the kinetic constants of the elementary steps and pulse number. The interplay of diffusion and kinetic behavior introduces complexity to the transient kinetic system, a data analysis method not having this considered will not be able to obtain accurate results, especially for realistic reaction with different diffusivities of molecules and complex reaction pathways. The VTAP model reliably maintains the fidelity of transient features.

Fig. 8 shows the RC petals of three sets of rate constants for the three-step reaction (eqn (9)). The set of rate constants for Fig. 8A–C has two characteristics: (i) slow backward reactions to guarantee the validation of all irreversible reaction steps and (ii) very fast forward reaction of steps 2 and 3 to guarantee the first step as rate-determining step. The catalyst surface remains empty during the VTAP simulation (Fig. 8C), since the two features mitigate the effects of coverage and reverse reaction. As a result, the rates of reactant A and product B both have linear correlations with the concentration of reactant A, as shown in Fig. 8A and B.

Fig. 8 RC petals and surface coverage of the three-step catalytic reaction (eqn (9)) with different rate constants in eqn (10)–(12). (A)–(C) are for k_1f = 1 × 10² s⁻¹, k_1b = 1 × 10⁻⁴ s⁻¹, k_2f = 1 × 10 s⁻⁴ s⁻¹, k_2b = 1 × 10⁻⁴ s⁻¹, k_3f = 1 × 10⁴ s⁻¹, k_3b = 1 × 10⁻⁴ s⁻¹. (D)–(F) are for k_1f = 1 × 10² s⁻¹, k_1b = 1 × 10¹ s⁻¹, k_2f = 1 × 10¹ s⁻¹, k_2b = 1 × 10¹ s⁻¹, k_3f = 1 × 10⁴ s⁻¹, k_3b = 1 × 10¹ s⁻¹. (G)–(I) are for k_1f = 1 × 10² s⁻¹, k_1b = 1 × 10¹ s⁻¹, k_2f = 1 × 10⁴ s⁻¹, k_2b = 1 × 10¹ s⁻¹, k_3f = 1 × 10¹ s⁻¹, k_3b = 1 × 10¹ s⁻¹.

The set of rate constants for Fig. 8D–F has relatively slow rate constant (k_2f) of second step (A* ↔ B*), making this step as rate-determining step. As shown in Fig. 8D, the rate-determining step in the middle decouples the linear correlation between consumption rate of A in the first step and production rate of B in the third step. Interestingly, the B production rate has linear correlation with surface coverage of A*, which is because (i) the second step (A* ↔ B*) being the rate-determining step and (ii) the very fast forward reaction of the third step (B* → B(g) + *, k_3f = 1 × 10⁴ s⁻¹) mitigating the effect of the reverse reaction of the second step (B* → A*). Comparatively, the RC petal in Fig. 8F has narrow shape, although the backward rate constant (k_3b = 1 × 10¹ s⁻¹) is same to the backward rate constant (k_2b = 1 × 10¹ s⁻¹) of the second step. Therefore, the shape of RC petal is not only determined by the rate constants of an elementary step; it is also affected by the competition from the other steps in the reaction mechanism.

The set of rate constants for Fig. 8G–I has relatively slow rate constant (k_3f = 1 × 10¹ s⁻¹) for forward reaction of the third step (B* → B(g) + *), making the desorption of product B from catalyst surface as the rate-determining step. This decouples the correlation of B production rate with rate, concentration and surface coverage of reactant A (Fig. 8G and H). Comparatively, the B production rate has narrow RC petal correlation with surface coverage of B*, which is the reactant of the third step.

The above discussion indicates that the shape of RC petals of catalytic reactions is largely affected by the rate constants for the same reaction mechanism. The rate-determining step, reversibility and the competition between reaction steps have been found to be nonnegligible factors in RC petal analysis. It should be noted that the three step reaction mechanism is just the simplest toy mechanism, a real catalytic reaction may have more complicated nonlinear reaction network and sometimes more than one type of active site on catalyst. VTAP has the capability of modeling such complex scenarios. We carried out VTAP simulation on an example (eqn (13)–(19)) containing two types of active sites and parallel reaction pathways, the modeling results of are shown in Fig. 9.

A(g) + * ↔ A* (13)

B(g) + * ↔ B* (14)

B(g) + # ↔ B# (15)

A* ↔ C* (16)

A* + B# ↔ D* + # (17)

C* ↔ C(g) + * (18)

D* ↔ D(g) + * (19)

Fig. 9 (A) The exit flux, (B) surface coverage, (C and D) RC petals of nonlinear reaction mechanism on catalyst containing two types of active sites. The forward rate constants for reactions 13–19 are set to 1 × 10² s⁻¹, the backward rate constant for reactions 13–17 are set to 1 × 10¹ s⁻¹, and the backward rate constants for reactions 18–19 are set to 1 × 10² s⁻¹.

In the reaction pathway of this nonlinear reaction mechanism, reactant A can adsorb on active site * and convert to product C, reactant B can adsorb on both * and # sites, but only the adsorbed B# can react with A* to form product D. In the simulation, same amount of A and B are pulsed in reactor together (pulse size = 1); the ratio between * and # sites on catalyst surface is set to 1/3. This results in the lower peak of exit flux of B than peak of A in Fig. 9A. It is interesting to point out that the yield of D is higher than C, although the rate constants of their formation steps are set to same. This can be elucidated by the coverage effect (as shown in Fig. 9B), since all species compete for * site, but # site is exclusive for B# adsorption, which is favorable for D formation.

The RC petals of reactants and products against concentration of reactant A are shown in Fig. 9C and D, respectively. Since the surface coverages of the reaction intermediates on catalyst surface and reaction rates of elementary steps in the reaction pathway transiently affect each other, and interpretation of exit flux and RC petal is not as straight forward as linear reactions or adsorption processes. Fig. 9C shows that, because of the existence of # site exclusive for B adsorption, the adsorption rate of B is faster than A and loses linearity with the concentration of A, although the rate constant of adsorption and pulse size are set to same values for A and B. The initial consumption rate of B is linearly correlated to its own concentration (Fig. S6). Fig. 9D shows that the RC petals of the rates of products C and D versus concentration of A are counterclockwise triangle shape, indicating no straight forward kinetic correlation. In comparison, the production rate of C has better correlation with surface coverage of A* (Fig. S7), and the production rate of D has better correlation with the product of surface coverages of A* and B# (Fig. S8).

Potential application

We have demonstrated detailed VTAP model analysis on several typical scenarios in adsorption processes, providing the fingerprints of reaction reversibility and dependence of catalyst surface coverage. These fingerprints can be directly used to identify the features in a chemical process in experiments and assist characterizing catalyst properties by analysis of the RC petals from Y-procedure⁹ and G-procedure.²² For example, a titration experiment using CO adsorption on a Pt/SiO₂ catalyst was conducted at 25 °C. The peak height of exit flux curves decreases with the increasing pulse number due to CO uptake, as shown in Fig. 10A. Fig. 10B demonstrates the effect of pulse number on the RC petal generated using the G-procedure.²² The adsorption rate decreases, and the gas concentration increases with the increasing pulse number. It indicates a lower apparent kinetic constant with increasing pulse number. It can be rationalized by a decrease in coverage of empty active sites on catalyst caused by CO uptake with increasing pulse number. Overall, the RC petals of this titration experiment indicate that CO adsorption on Pt/SiO₂ shows an irreversible feature; the desorption process is negligible over the duration of the experiment, and it is much slower than the adsorption process. Fig. 10C shows exit flux curves of CH₄ adsorption on a MoC_x/ZSM5 catalyst conducted at 700 °C; Fig. 10D shows the RC petal generated using the Y-procedure.⁹ The RC petal shows a fingerprint of reversible adsorption of CH₄ on MoC_x/ZSM5 as the reaction rate switches from positive to negative values. The reversible feature can be predicted by weak physisorption of CH₄ in zeolite pores. It should be noted that Y and G procedures assume uniformity in catalyst zone, which is an ideality and can lead to deviations from real experiments. VTAP modeling does not use this assumption, and one can examine where the size of the catalyst zone may exceed a certain tolerance threshold. Of course, in the real experiment, one has only very coarse control over the thickness of the catalyst zone but understanding the degree of nonuniformity to expect is very useful. The systematic comparison of RC petals from VTAP, Y and G procedure has been carried out in the recent paper.²³

Fig. 10 (A) Exit flux curves and (B) RC petals of CO titration experiments on a Pt/SiO₂ catalyst at 25 °C. (C) Exit flux curves and (D) RC petals of CH₄ pulses on a MoC_x/ZSM5 catalyst at 700 °C.

The VTAP model minimizes the use of approximations and enhances the flexibility in modeling the different settings of transient reactor experiments, such as the catalyst loading amount and location. This improves the ability to model more complex experiments. For example, Mirena et al. carried out pulse experiments under contrasting transport regimes, showing that a given amount of catalyst particles exhibits higher activity with increasing degrees of separation.²⁴ The unique feature of the VTAP model is the vectorization along the reactor coordinate. This feature removes simplifying assumptions, such as the thin zone TAP reactor model and allows for much more detailed simulations of the particle shadowing effect induced by the competition between closely neighboring active particles when catalyst particles are placed in different locations relative to each other.

Extracting intrinsic information from measured exit flux has been a challenge for not only the TAP experiments but also the general catalysis community for decades. Y-procedure and G-procedure have been developed for this purpose. It should be noted that both Y-procedure and G-procedure are under a major approximation of uniformity in catalyst zone, meaning that the concentration, coverage and reaction rate are same in all positions. Y-procedure is sensitive to data noise in practical applications and relies on significant curve smoothing treatment. The fitting of gamma distribution in G-procedure causes systematic deviation. The VTAP model can be used to extract intrinsic information by fitting experimental exit flux curves to validate rate expressions and extract the kinetic constants. Curve fitting for obtaining rate constants and rate expressions has been a widely used conventional method in chemical kinetics.

In addition to fitting the experimental exit flux curves, we have been working on a different method, the combination of VTAP model with machine-learning algorithms. This method uses artificial neural networks to learn the correlations in the high dimensional data generated by VTAP modeling. The trained artificial neural networks are then used in the back propagation process for directly predicting the rate and concentration near the catalyst from experimental exit flux curve without needing to hypothesize the reaction mechanism. This method has been applied in interpreting CO titration experiments on Pt/Al₂O₃ and Pt/SiO₂ catalysts.²⁵

The FORTRAN code TAPFIT,²⁶ the Python-based TAPSolver,¹⁶ and the MATLAB code SimTAP⁶ have been previously developed by different institutions for the purpose of modeling TAP experiments; even more variations exist privately within other TAP research groups. Each of these tools has its own advantages. While there are more intricacies to each, at a high level one may use SimTAP or VTAP for fast, high-fidelity forward simulations over rich pulse trains (using VTAP if understanding fields, gradients, or placement effects are important). Both account for any residuals in gas phase concentration and surface coverage from the previous pulse and long multipulse sequences can be studied. TAPSolver is a good choice when gradient-based parameter estimation with sensitivities and thermodynamic constraints is important. TAPFIT is a mature simulator with built-in Levenberg–Marquardt fitting and second-order statistical regression for conditioning of replicates.

VTAP adds to these more factors including the catalyst loading, different types of active sites, complex reaction pathways such as nonlinear reaction networks. VTAP also considers and obtains concentration and surface coverage gradients in catalyst zone, mitigating the error from assumption of uniformity in the catalyst zone. These features increase the practical capabilities for analysis of complexity in catalysis and surface science As an example, we have applied VTAP model to ammonia decomposition reaction, which includes the adsorption of NH₃, dehydrogenation to form H* and N* and desorption of H₂ and N₂ in the reaction pathway. Fig. 11 shows the calculated surface coverages of the intermediates with respect to time, indicating that the dominant NH_x* (x = 0–3) species is NH₃* in the beginning; it becomes N* after 0.1 second due to dehydrogenation reaction. The rate constants of the elementary steps in ammonia decomposition can be obtained by fitting the model with experimental exit flux curves and carrying out multi-scale modeling using DFT-calculated rate constants.

Fig. 11 The surface coverages of reaction intermediates in NH₃ decomposition reaction. The rate constants of adsorption processes are set to 10³ s⁻¹, and the rate constants of desorption and reaction steps are set to 10² s⁻¹.

Conclusion

The VTAP model offers new features for handling complexity in transient experiments. It allows for a more comprehensive analysis of the correlations between transport and kinetic factors influencing the reaction processes. This simulation framework enables a more detailed interpretation of experimental data, revealing how different parameters – number of active sites, reaction rate constants, and catalyst loading and distribution in the reactor – influence the observed exit flux. A webpage application (WebApp) version of VTAP is created for simulation of adsorption processes in transient experiments without the need for coding experience. A screenshot, description, and internet address of the WebApp can be found in the SI. Advanced versions of the WebApp will be continuously developed and updated for more complex reactions. Since the TAP reactor (a diffusion reactor) is so distinctive from the commonly encountered plug flow reactor, these simulations can also be used as a learning tool for understanding the complex interplay of transport and kinetics in dynamic experiments.

The VTAP model accurately captures time-dependent transport in the reactor and demonstrates how the pulsed gas quickly relaxes to flat concentration profile within several milliseconds, creating a well-mixed reactor needed for precise kinetic characterization. In adding adsorption, reaction, and desorption processes to the transport equation, the model demonstrates how rate/concentration analysis, the RC petal, can reveal clear fingerprints indicating the balance of these rate constants together with the number of active sites on a catalyst sample. The characteristics that arise from other factors, such as catalyst loading, accumulating pulses, reactor length, catalyst zone length and position, are also distinguished by the model. In the case of the three-step reaction mechanism, one of the most surprising features was the concentration profile of products predicted by the VTAP reactor model. Here the impact of different initial conditions becomes clear: diffusion of the reactant starts from the inlet pulse and diffusion of the product starts from the concentration time-dependence of the reactant.

The interplay of transport and reaction in this dynamic experiment is both complex and non-intuitive. The VTAP model provides a detailed visualization of the changing reactor concentration profiles and the time-dependent evolution of the catalyst surface during the measurement. As such, the model creates a foundation for more comprehensive analysis of multistep reaction systems and lays the groundwork for future studies aimed at resolving mechanistic details.

Author contributions

SW: conceptualization, python coding, writing original manuscript. HC: experiment, reviewing and editing manuscript. BAT: experiment, reviewing and editing manuscript. SK, JMP, DM, DS, MRK and RVR: resources, reviewing and editing manuscript. RF: funding acquisition, conceptualization, Writing and editing manuscript. All authors have given approval to the final version of the manuscript.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Supplementary information a document containing: effect of reactor length on exit flux; modeling results for case of less catalyst (site density = 0.1); effect of the thickness and position of catalyst zone; rate-concentration petals of product in A ↔ B reaction; rate-concentration petals of A + B ↔ C + D reaction; the WebApp version of VTAP; speed test. Animation figures: standard diffusion; adsorption process; adsorption process with less catalyst (site density = 0.1); multiple-pulse with irreversible adsorption and multiple-pulse with reversible adsorption. See DOI: https://doi.org/10.1039/D5RE00280J.

The data supporting this article titled A Simulation Framework for Extracting Intrinsic Kinetics in Transient Experiments have been included as part of the supplementary information (SI).

Acknowledgements

This work was supported by the U.S. Department of Energy (USDOE), Office of Energy Efficiency and Renewable Energy (EERE), Industrial Efficiency and Decarbonization Office (IEDO) Next Generation R&D Project DE-FOA-0002252-1775 under contract no. DE-AC07-05ID14517.

References

1. J. T. Gleaves, G. S. Yablonskii, P. Phungphai and Y. Schuurman, TAP-2: An interrogative kinetics approach, Appl. Catal., A, 1997, 160, 55–88.
2. S. O. Shekhtman and G. S. Yablonskii, Thin-Zone TAP Reactor versus Differential PFR: Analysis of Concentration Nonuniformity for Gas-Solid Systems, Ind. Eng. Chem. Res., 2005, 44, 6518–6522.
3. Y. Wang, M. R. Kunz, D. Constales, G. Yablonsky and R. Fushimi, Rate/Concentration Kinetic Petals: A Transient Method to Examine the Interplay of Surface Reaction Processes, J. Phys. Chem. A, 2019, 123(40), 8717–8725,  DOI:10.1021/acs.jpca.9b05911.
4. G. S. Yablonsky, D. Constales, S. O. Shekhtman and J. T. Gleaves, The Y-procedure: How to extract the chemical transformation rate from reaction–diffusion data with no assumption on the kinetic model, Chem. Eng. Sci., 2007, 62(23), 6754–6767,  DOI:10.1016/j.ces.2007.04.050.
5. M. R. Kunz, T. Borders, E. Redekop, G. S. Yablonsky, D. Constales, L. Wang and R. Fushimi, Pulse response analysis using the Y-procedure: A data science approach, Chem. Eng. Sci., 2018, 192, 46–60,  DOI:10.1016/j.ces.2018.06.078.
6. L. Brandão and C. Reece, Non-steady state validation of kinetic models for ethylene epoxidation over silver catalysts, Catal. Sci. Technol., 2024, 14(13), 3596–3608,  10.1039/d4cy00052h.
7. M. R. Kunz, A. Yonge, Z. Fang, R. Batchu, A. J. Medford, D. Constales, G. Yablonsky and R. Fushimi, Data driven reaction mechanism estimation via transient kinetics and machine learning, Chem. Eng. J., 2021, 420, 129610,  DOI:10.1016/j.cej.2021.129610.
8. J. Perezramirez and E. Kondratenko, Evolution, achievements, and perspectives of the TAP technique, Catal. Today, 2007, 121(3–4), 160–169,  DOI:10.1016/j.cattod.2007.01.001.
9. E. A. Redekop, G. S. Yablonsky, D. Constales, P. A. Ramachandran, C. Pherigo and J. T. Gleaves, The Y-Procedure methodology for the interpretation of transient kinetic data: Analysis of irreversible adsorption, Chem. Eng. Sci., 2011, 66(24), 6441–6452,  DOI:10.1016/j.ces.2011.08.055.
10. Y. Schuurman, Assessment of kinetic modeling procedures of TAP experiments, Catal. Today, 2007, 121(3–4), 187–196,  DOI:10.1016/j.cattod.2006.06.046.
11. D. Constales, G. S. Yablonsky, G. B. Marin and J. T. Gleaves, Multi-zone TAP-reactors theory and application: I. The global transfer matrix equation, Chem. Eng. Sci., 2001, 56, 133–149.
12. D. Constales, G. S. Yablonsky, G. B. Marin and J. T. Gleaves, Multi-zone TAP-reactorstheory and application: II. The three-dimensional theory, Chem. Eng. Sci., 2001, 56, 1913–1923.
13. D. Constales, S. O. Shekhtman, G. S. Yablonsky, G. B. Marin and J. T. Gleaves, Multi-zone TAP-reactors theory and application IV. Ideal and non-ideal boundary conditions, Chem. Eng. Sci., 2006, 61(6), 1878–1891,  DOI:10.1016/j.ces.2005.10.022.
14. D. Constales, G. S. Yablonsky, G. B. Marin and J. T. Gleaves, Multi-zone TAP-reactors theory and application. III Multi-response theory and criteria of instantaneousness, Chem. Eng. Sci., 2004, 59(17), 3725–3736,  DOI:10.1016/j.ces.2004.05.023.
15. S. O. Shekhtman, G. S. Yablonsky, S. Chen and J. T. Gleaves, Thin-zone TAP-reactor - theory and application, Chem. Eng. Sci., 1999, 54, 4371–4378.
16. A. Yonge, M. R. Kunz, R. Batchu, Z. Fang, T. Issac, R. Fushimi and A. J. Medford, TAPsolver: A Python package for the simulation and analysis of TAP reactor experiments, Chem. Eng. J., 2021, 420, 129377,  DOI:10.1016/j.cej.2021.129377.
17. T. S. Kim, C. R. O'Connor and C. Reece, Interrogating site dependent kinetics over SiO(2)-supported Pt nanoparticles, Nat. Commun., 2024, 15(1), 2074,  DOI:10.1038/s41467-024-46496-1.
18. R. Roelant, Mathematical determination of reaction networks from transient kinetic experiments, PhD thesis, Ghent University, Ghent, Belgium, 2011.
19. A. Goguet, C. Hardacre, B. Inceesungvorn, K. Morgan and S. O. Shekhtman, Correction for a possible reversible adsorption over an “inert” material, Catal. Sci. Technol., 2011, 1(5), 760–767,  10.1039/c0cy00075b.
20. J. Piiper and H. Worth, Value and Limits of Graham's Law for Prediction of Diffusivities of Gases in Gas Mixtures, Respir. Physiol., 1980, 41, 233–240.
21. E. A. Redekop, G. S. Yablonsky, V. V. Galvita, D. Constales, R. Fushimi, J. T. Gleaves and G. B. Marin, Momentary Equilibrium in Transient Kinetics and Its Application for Estimating the Concentration of Catalytic Sites, Ind. Eng. Chem. Res., 2013, 52(44), 15417–15427,  DOI:10.1021/ie400677b.
22. M. Ross Kunz, T. Borders, E. Redekop, G. S. Yablonsky, D. Constales, L. Wang and R. Fushimi, Pulse response analysis using the Y-procedure: A data science approach, Chem. Eng. Sci., 2018, 192, 46–60,  DOI:10.1016/j.ces.2018.06.078.
23. S. Wang, H. Chau, S. Kristy, B. A. Thompson, J. P. Malizia and R. Fushimi, Machine Learning-Assisted Recovery of Delicate Kinetic Information from Transient Reactor Experiments, ACS Eng. Au, 2025, 5(3), 298–310,  DOI:10.1021/acsengineeringau.5c00025.
24. J. I. Mirena, E. Redekop, H. Poelman, N. V. Srinath, D. Constales, G. B. Marin, G. S. Yablonsky, J. T. Gleaves and V. V. Galvita, Shadowing Effect in Catalyst Activity: Experimental Observation, ACS Catal., 2022, 12(9), 5455–5463,  DOI:10.1021/acscatal.2c00818.
25. S. Wang, H. Chau, S. Kristy, B. A. Thompson, J. P. Malizia and R. Fushimi, Machine Learning-Assisted Recovery of Delicate Kinetic Information from Transient Reactor Experiments, ACS Eng. Au, 2025, 5(3), 298–310,  DOI:10.1021/acsengineeringau.5c00025.
26. E. A. Redekop, G. S. Yablonsky, D. Constales, P. A. Ramachandran, J. T. Gleaves and G. B. Marin, Elucidating complex catalytic mechanisms based on transient pulse-response kinetic data, Chem. Eng. Sci., 2014, 110, 20–30,  DOI:10.1016/j.ces.2013.11.050.

---