title: TAPsolver: A Python package for the simulation and analysis of TAP reactor experiments
authors: Adam Yonge M. Ross Kunz Rakesh Batchu Zongtang Fang Tobin Issac Rebecca Fushimi Andrew J. Medford
text: 
  Sensitivity of the Polyakov loop and related observables to chiral symmetry restoration
      Mugdha Sarkar
      December 30, 2023
  =======================================================================================
  
  
  
  
  An open-source, Python-based Temporal Analysis of Products (TAP) reactor simulation and processing program is introduced. TAPsolver utilizes algorithmic differentiation for the calculation of highly accurate derivatives, which are used to perform sensitivity analyses and PDE-constrained optimization. The tool supports constraints to ensure thermodynamic consistency, which can lead to more accurate parameters and assist in mechanism discrimination. The mathematical and structural details of TAPsolver are outlined, as well as validation of the forward and inverse problems against well-studied prototype problems. Benchmarks of the code are presented, and a case study for extracting thermodynamically-consistent kinetic parameters from experimental TAP measurements of CO oxidation on supported platinum particles is presented. TAPsolver will act as a foundation for future development and dissemination of TAP data processing techniques.
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  Catalyst selection and optimization is a necessary step in the design of many chemical production facilities <cit.>. Determining an appropriate catalyst for a process, however, is often challenging due to the vast number of variables affecting the catalyst and its environment <cit.>. Before determining appropriate production-scale operating conditions, the choice of catalyst support <cit.>, composition <cit.>, cluster size <cit.>, and solvent <cit.> must be considered to optimize product selectivity and yield. Even mild changes to catalyst performance can result in multi-million-dollar profit shifts or substantial changes in pollutant emissions <cit.>. A robust knowledge of intrinsic kinetic parameters can lead to the rational design of both the catalyst and reactor <cit.>. For these reasons, methods to efficiently probe the intrinsic kinetic properties of catalysts have been, and continue to be, developed <cit.>.
  
  Solving for kinetic parameters through reaction model fitting is a standard approach in the field of catalysis <cit.>. The reliability of these fitted parameters, however, are limited by the availability of the experimental data points. For a steady-state experiment, the model is typically sensitive to only a few parameters because the rate-limiting step dominates the extracted kinetic information <cit.>. This means it is possible to fit multiple kinetic models with comparable quality, leading to uncertainty in the resulting parameters and reaction mechanism <cit.>. Density functional theory (DFT) <cit.> calculations, transition state theory (TST) <cit.> and micro-kinetic modeling (MKM) <cit.> are frequently used to address this challenge, leading to fundamental insights about the thermodynamics and kinetics of elementary reactions on catalysts that go beyond what can be deduced from steady-state kinetic experiments. These methods can provide users with information pertaining to the mechanism <cit.>, turnover frequency <cit.>, molecular orientation <cit.> and stability <cit.>, but running simulations with DFT requires significant computational resources, even for relatively simple systems <cit.>. Moreover, DFT simulations require an accurate representation of the relevant active sites, which is often unknown for complex catalyst formulations <cit.>, and the accuracy of the DFT calculations can also be insufficient for quantitatively accurate predictions <cit.>. Surface-science experiments can be used to measure adsorption energies and intrinsic kinetics directly, but they require even more time and resources, and are typically limited to simplified single-crystal surfaces.
  
  
  
  Transient kinetic experiments are an additional technique that can complement steady-state experiments, computational chemistry, and traditional surface science experiments <cit.>. Transient experiments can provide users with a nuanced understanding of the intrinsic kinetics of a system <cit.>, but effectively decoupling the transport and kinetics is often complex <cit.>. The temporal analysis of products (TAP) reactor system was developed to provide well-defined transport properties and enable kinetic studies of complex supported catalysts. TAP experiments consist of a series of nano-mole pulses operating at under isothermal conditions within the Knudsen diffusion regime, and data are collected at a millisecond time resolution <cit.>. This provides rich and detailed insight into the transient kinetic response of supported catalysts. However, a noteable challenge of working with TAP data is the complex data analysis required to convert the raw data to meaningful kinetic information. Steps have been taken through the years to account for these issues, with the development of the governing partial differential equations (PDEs) <cit.>, methods to solve them, and examples of solutions for different transient experiments <cit.>. This has provided a foundation for prior TAP simulators <cit.>. The number of methods have expanded over the years, as presented in Figure <ref>), with the development of approaches for the construction of mechanisms, simulation of pulses, fitting of parameters and data-driven models for analyzing results  <cit.>.
  
  Although several TAP analysis packages exist in the field, some obstacles still remain. Most codes are not easily or publicly accessible, instead requiring users to contact developers directly to gain access to code, and many make use of proprietary frameworks or programming languages that require costly licences. In addition, validation and documentation can be sparse or convoluted, making it difficult to confirm the accuracy of results without frequent discussion with the original developers. 
  
  In the case of TAP fitting programs, the optimization algorithms and loss functions are often hard-coded, making it difficult to implement emerging optimization techniques from machine learning. The TAP fitting codes also require strong initial guesses and do not provide estimates of uncertainty, making it challenging to draw strong conclusions about the fitted parameters or mechanisms. 
  To address these issues, the Python program TAPsolver, a TAP reactor simulation and analysis tool, is developed and presented. This program will be made freely available for academic research. TAPsolver is built around the open-source FEniCS <cit.> and Dolfin-Adjoint <cit.> packages, which solve PDEs via finite element methods and utilize model adjoints to enable algorithmic differentiation. TAPSolver is designed to be compatible with external optimization routines, enabling experimentation with various algorithms, and constraints can easily be added to the loss function or optimization algorithm to enforce thermodynamic consistency or other known information.
  TAPSolver flexibly accepts different reaction mechanisms and rate/reactor parameters, and interfaces with other techniques for TAP data analysis, making it a natural choice integrating future developments in TAP data analysis.
  In this work, the commands used to define TAPsolver data processing are introduced and TAPsolver is applied to both synthetic and experimental data. 
  
  
  
  § MATHEMATICAL BACKGROUND
  
  
  The fundamental PDEs and boundary conditions used for TAP simulations and constraints applied to optimize rate constants and enforce thermodynamic consistency
  are briefly outlined in this section. 
  More rigorous details can be found in publications by Van der Linde et al. <cit.>, Farrell et al. <cit.> and Mhadeshwar et al. <cit.>, or elsewhere <cit.>.
  
  
  
   §.§ Defining the Transport and Reaction Equations
  
  
  Mathematical descriptions of catalytic experiments are often complex, involving multidimensional heat, mass and momentum transfer at different temporal and spatial scales. Current limitations of computational resources often force researchers to balance simplifying assumptions with accuracy. TAP reactors, however, are setup to reduce many of these complexities. TAP primarily operates in the Knudsen diffusion regime and consists of two equally sized inert zones surrounding a thin catalytic region. Though reactions in the catalyst zone make the system inherently complex and can lead to challenging simulations, there are simplifications that can limit this burden. It has been shown that a one-dimensional model is acceptable when the length of the TAP micro-reactor is at least three and a half times the radius of the reactor <cit.>. This drastically reduces the complexity of the equations and the associated computational expense of solving them. The PDE describing the diffusion of a gas species through a TAP reactor, with a general reaction term, is
  
  
      ε_j∂ C_i(x,t)/∂ t - ∂/∂ x·(D_i,j ∇C_i(x,t)) = R(C_i(x,t),θ_i(x,t),k_j,...)
  
  
  where C_i(x,t) represents the concentration of gas species i, ε_j is the void fraction of the material in zone j, t is time, x is the length of the reactor, θ_i are surface intermediates i and D_i,j is the Knudsen diffusion coefficient of gas i in zone j. The change in surface intermediates with time is defined as 
  
  
      ∂θ_i(x,t)/∂ t = R_i(C_j(x,t),θ_j(x,t),k_l,...)
  
  
  The reaction term R is the microkinetic model, defined as <cit.>
  
  
      R_i = ∑_m s_im(k_m^+∏_j ∈ F_mθ_j(x,t) ∏_j ∈ F_m C_j(x,t) - k_m^-∏_l ∈ B_mθ_l(x,t) ∏_l ∈ B_m C_l(x,t))
  
  
  where s_im is the stoichiometry of species i in elementary step m, k_m^+/- are the forward/reverse rate constant for elementary step m, and F_m and B_m are sets of indices corresponding to the intermediate and gas-phase species in the forward and backward reactions of elementary step i. Solving this PDE requires specifying  boundary conditions. In the case of the TAP experiment, the initial condition of gas species are
  
  
      C_g(x, t_p) = δ(x, t_p)
  
  
  where δ(x, t_p) is a delta function 
  introduced at the pulse time t_p. The initial condition of surface intermediates inside the catalyst zone is
  
  
      θ(P_1 < x < P_2, t_0) = U_i
  
  
  where U_i is the initial surface concentration of adsorbate or active site i and P_1 and P_2 represent the normalized domain of the catalyst zone in the reactor. These values can be changed based on the fraction the reactor occupied by the catalyst. The initial condition of surface intermediates outside the catalyst zone is
  
  
      θ(x < P_1 and  x > P_2, t_0) = 0
  
  
   where no surface concentration is considered due to assumptions of inert interactions with non-catalytic materials. At the entrance of the reactor, a Neumann boundary condition can be imposed due to the use of a pulse valve and is written as
  
  
      ∂ C_g(0,t_p)/∂ x = 0
  
  
  while a vacuum is being applied at the exit of the reactor, leading to the Dirichlet boundary condition
  
  
      C_g(L, t) = 0
  
  
  with L being the length of the reactor. Multiple species can be introduced into the reactor system simultaneously, indicated by the index g. When a species is formed in the catalyst region but not pulsed into the reactor, then the initial intensity of the Dirac delta function associated with that species is set to zero. All gas species will have a concentration of zero at the exit of the reactor due to the applied vacuum conditions, and the outlet flux is defined as 
  
  
      F_i = D∂ C_i(x,t)/∂ x_x=L
  
  
  The corresponding variational form for the transport term of Equation <ref>, needed for solving the PDEs using finite-element methods, is written as <cit.>
  
  
      ∫( ε_j(C_i^n+1-C_i^n)v_i/t + D_i,j∇c_i^n+1· ∇v_i)  dx = 0
  
  
  where C_i^n+1 and C_i^n are the species concentrations at current and previous time steps, respectively, Δt is the size of the time step and v_i is the test function for each of the gas species. The terms test and trial functions are frequently used in variational problem notation and represent the components for performing integration by parts over the spatial domain. Reactive terms can be added to the right-hand side of Equation <ref> and consist of combinations of the concentration terms. As a specific example, consider the simple gas-phase reaction
  
  
      A → B
  
  
  with a rate constant k.  The kinetic model is defined as
  
  
      kC_Av_Adx - kC_Av_Bdx = 0
  
  
  where each term represents the consumption and generation of species A and B, respectively. When combined, the variational form of the reaction-diffusion equation is defined as
  
  
      ∫( (C_A^n+1-C_A^n)v_A/t + D_A∇c_A^n+1· ∇v_A)  dx 
          + ∫( (C_B^n+1-C_B^n)v_B/t + D_B∇c_B^n+1· ∇v_B)  dx = kC_Av_Adx - kC_Av_Bdx
  
  
  
  
   §.§ Curve Fitting and Thermodynamic Consistency
  
  
  One key advantage of TAP data is that the TAP curve will implicitly contain information about multiple elementary steps, and TAP curve fitting has been commonly used to extract intrinsic kinetic parameters from TAP data <cit.>. Fitting TAP curves to extract rate constants is a form of PDE-constrained optimization, and requires an objective function to optimize. 
  Agreement between model and experiment can be optimized by using an objective function corresponding to the summation of the difference between each of the simulated and experimental outlet flux data points, defined as
  
  
      J_data = ∑_i ( 1/2D ∫_L-x^L | u_i(x,t) - u_i,obs(x,t)|^2 dx )
  
  
  where i represents each of the points considered, L is the length of the reactor, u_i is the outlet concentration of the forward solution, u_i,obs is the experimental outlet concentration, and x is the mesh step size. Since the experimental data corresponds to the outlet flux, the inner product is evaluated over the spatial step at the reactor outlet, i.e. L - x to L. Other objective functions could be defined that involve additional experimental data or observable values (e.g. surface species, spectroscopic information, integral quantities) <cit.>. Since outlet fluxes are the most common observed and fitted values for TAP experiments, this provides a natural starting point.
  
  In addition to agreement with the data, it is also important to consider the thermochemistry of the reaction. Any set of kinetic parameters that does not obey the proper thermodynamic constraints will lead deviations in the predicted equilibrium concentrations and errors with the energy balance will arise <cit.>. In addition, Equation <ref> alone may not always be adequate for identifying the appropriate reaction mechanisms or kinetics, so the additional constraint may aid in mechanism discrimination. For example, micro-kinetic analysis through DFT and fitting results of steady-state kinetic experiments often utilizes thermodynamic consistency constraints. <cit.> These approaches have been used to help test the validity of many studies, including the partial oxidation of methane <cit.>, ethylene hydrogenation <cit.> and the water-gas shift reaction <cit.>.
  
  Variations of the thermodynamic consistency equations have been outlined by Dumesic et al. <cit.> Although the implementation of each form of thermodynamic consistency would be useful, many forms require the definition or estimation of intermediate values (i.e. the thermodynamic data per elementary processes) that may not be known. However, the thermochemistry of the overall reaction, visualized in Figure <ref>, is often well-known, so it is useful to be able to ensure that the optimum reaction parameters are consistent with this information, which will be independent of the catalyst material used.
  
  
  
  The free energy of reaction is defined as 
  
  
      ΔG_gas = ΔH - T*ΔS
  
  
  where the enthalpic and entropic terms for individual gas species are summed to define H and S, respectively. To determine the contributions of individual elementary reactions (denoted i), the equilibrium constant is defined as 
  
  
      K_i = k_f,i/k_b,i
  
  
  and the free energy of the process is calculated through
  
  
      ΔG_i = -R*Tln(K_i)
  
  
  with R representing the ideal gas constant and T the reactor temperature. The difference between the free energy of reaction and each of the elementary free energies is the additional term of the objective function defined in Equation <ref> to account for thermodynamic consistency, written as 
  
  
      J_thermo = (ΔG_gas - ∑_i ∈ MΔG_i)^2
  
  
  where M is the set of elementary steps for a given serial mechanism corresponding to the gas-phase reaction.
  A combined objective function, given by 
  
  
      J = J_data + α J_thermo
  
  
  can then be utilized to simultaneously fit experimental TAP data and satisfy thermodynamic constraints, with the parameter α controlling the relative importance of thermodynamics compared to data.
  Thermodynamics have been considered alongside TAP reactor studies previously <cit.>, but as an additional component in interpreting results rather than constraining optimized parameters. To our knowledge, no other TAP optimization code supports thermodynamic consistency constraints. 
  
  
  
  § TAPSOLVER IMPLEMENTATION
  
  
  The FEniCS package is an established code for efficient finite element simulation, and the associated Dolfin-Adjoint package utilizes model adjoints to provide algorithmic differentiation of the models. One challenge in TAP simulations is that significant modifications to the underlying PDEs are required to simulate different microkinetic mechanisms. 
  Manually changing these equations for different experimental setups or chemical systems is time consuming and error-prone. 
  However, automating the construction of these equations provides an efficient route to making FEniCS practical for TAP reactor analysis. Similarly, having to manipulate code directly to set initial conditions (e.g. pulse intensities, surface concentrations, etc.) or analysis methods (e.g. curve fitting, sensitivity analysis) can also become a bottleneck.  To overcome these issues, TAPsolver has a spreadsheet-based input file to simplify the construction of these equations, and the code is modular, facilitating facile access to various analysis methods. This section outlines each of the components of the input file and the available analysis techniques. 
  
  
  
  
  
   §.§ Defining Reactor Setup, Feed and Surface Composition
  
  
  The reactor setup, feed and surface composition, and reaction mechanism specification are the primary components of the TAPsolver input file, presented in Table <ref>, with thermodynamic consistency being an optional feature. Many of the variables defined in the input file are repeatedly used, enabling files to be copied and minimally modified. A brief discussion of the key inputs is provided here, and more details are available in the documentation.
  
  Several parameters are used to define the layout of the simulated TAP reactor, including the reactor length, radius, void fraction, catalyst location and fraction, and temperature. For operation in the Knudsen transport regime, a reference diffusion must be defined, including the reference temperature and mass in both the inert and catalyst domains. This information is used to determine the diffusion coefficient for all gaseous species based on the mass of each through scaling <cit.>.
  In addition to variables used to define the reactor setup are the variables that to control the simulation accuracy. The primary accuracy variable is the the 'Mesh Size' which defines the spatial precision.
  TAPsolver also supports mesh refinement within the catalyst domain. Simulating transport through inert zones requires a much coarser mesh than the catalyst domain where reactions happen. 
  The 'Catalyst Mesh Density' input parameter refines the precision of the mesh in the catalyst zone while keeping the mesh in both inert zones constant. For example, a mesh size to simulate transport in a catalyst free reactor would be approximately two hundred cells. If the user wanted to run a simulation with fifty cells within the catalyst zone that occupied 2% of the reactor, the simulation would require a uniform mesh size of 2500 cells. With the mesh density parameter, it is possible to have a catalyst mesh size of sixty-four (with the initial four cells in the catalyst zone being doubled four times), while the total reactor mesh size is 264. This can substantially improve accuracy with a minimal impact on simulation time.
  
  The reactant pulse intensity, pulse time, mass, and initial surface composition are defined in the input file below their associated gas names. 
  The user can also specify inert gas pulses that are not involved in the reaction process to the right of reactive gas species. Pulses can also be introduced at later times, enabling the simulation of multi-pulse or “pump-probe” experiments. The 'Output Folder' and 'Experimental Data Folder' variables specify the file structure of analysis results and directly link the experimental data to the simulations, respectively. A series of nested folders are generated following each simulation that contain all data for the processes specified by the user. 
  'thin data', 'sensitivity analysis' and 'uncertainty quantification'. 
  The input file used to run the simulation is also stored in the output folder to improve data veracity.
  
  
  
   §.§ Defining Elementary Reactions in TAPsolver
  
  
  The set of elementary reactions and list of gaseous species, adsorbed species, and active sites defined by the user are converted into a matrix defining the microkinetic model <cit.> and subsequently used to produce the variational form that is solved by FEniCS. 
  The elementary reactions defined in the input file are flexible, with reversible and irreversible reactions being defined through the arrow (either a reversible or right pointing arrow) and the values of the kinetic parameters being placed in the cells to the right of the elementary reactions, with the first cell representing the forward rate constant and the second cell representing the reverse rate constant. Three options are available for the definition of kinetic parameters in the input file, including direct rate constants, pre-exponential/activation enthalpy (noted with a ”$” symbol), or activation free energy (noted with a ”@” symbol). Users can also select which of these kinetic parameters to include or exclude in the optimization routine through the inclusion of an exclamation point following the value of the parameter, making it easier to target specific values of interest. Examples of the previous parameter options are shown in Table <ref>, with a reversible, irreversible, activation enthalpy based, free energy based and fixed parameter values being presented through each elementary reaction, respectively.   
  
  
  
  -5pt
  |c|c|c|
  Examples of how to define elementary reactions and their kinetic parameters for use by the TAPsolver input file.  
  
  Elementary Reaction      Forward     Reverse 
  
  A + * ↔ A^*     1
       2 
  
  B + * → B^*      1
       – 
  
  A* + B* → C^*      1$2     – 
  
  A* + B* → A^*     1@2     – 
  
  C^*→ C + *     1!     – 
  
  
  
  
  
  
   §.§ Running Analyses
  
  
  TAPsolver is structured to enable easy access to various features and analysis routines via Python functions.
  A snippet of TAPsolver code is presented in Figure <ref>, with examples implementations of some implemented features. Forward simulations of TAP pulses can be performed with the  command along with the name of the input file, the simulation time () and the number of pulses (). Once the simulation finishes the resulting outlet fluxes can be plotted using the  function. 
  Desired subsections of the simulated and experimental data can be specified with the parameter 'pulse', making it easier to observe trends and visualize results along with experimental data. Methods related to the inverse problem can also be programatically utilized. Parameters can be fitted with the  with similar inputs to the forward problem.  
  The gradient of the objective function with respect to fitted parameters can also be calculated through the  command. One feature of the sensitivity analysis is the option to compute the sensitivity over time (the influence of each kinetic parameter with respect to the objective function at each time step), as opposed to being defined in terms of the total objective function (used in optimization). This can be achieved using the sensType argument in , and enables the user to identify which parts of the TAP curve are controlled by specific elementary steps.
  
  
  
  
  
  § RESULTS AND DISCUSSION
  
  
  Both the forward and inverse problems are probed on simulated and experimental data. The complexity of the simulations ranges from pure diffusion to a mixed Eley-Rideal / Langmuir-Hinshelwood mechanism for the carbon monoxide oxidation reaction. TAPsolver is bench-marked and validated against analytical solutions to show how it performs for different levels of accuracy and simulation complexities. Direct comparisons between the finite difference approximation and adjoint approach for computing derivatives are provided. A case study illustrating fitting to experimental TAP data for CO oxidation over supported platinum particles is also provided.
  
  
  
   §.§ TAPSolver Validation
  
  
  
  
    §.§.§ Validation of the Forward Problem
  
  
  To confirm that TAPsolver is accurately simulating TAP reactor pulses, two examples of validation are provided: pure diffusion and irreversible adsorption. These scenarios are selected because analytical solutions are known, while exact solutions for more complex processes do not exist. 
  
  The analytical solution for the diffusion outlet flow in a TAP reactor is 
  
  
      F̅_A/N_pA = D_eAπ/ε_bL^2∑_∞^n=0(-1)^n(2n+1)exp(-(n+0.5)^2π^2t D_eA/ε_bL^2)
  
  
  where N_pA is the pulse intensity, D_eA is the diffusion coefficient of species A, ε_b is the void fraction of the material, L is the length of the reactor and t is the time <cit.>. The analytical solution for irreversible adsorption with diffusion in a TAP reactor is
  
  
      F̅_A/N_pA = D_eAπ/ε_bL^2exp(-k_a^')τ∑_∞^n=0(-1)^n(2n+1)exp(-(n+0.5)^2π^2τ)
  
  
  where -k_a^' is the rate characteristic for a unitary adsorption process and τ is the mean residence time <cit.>. Equations <ref> and <ref> are similar and the core difference lies in the inclusion of the rate constant at the beginning of Equation <ref>.
  
  The left panel of Figure <ref> shows the comparison between the analytical diffusion solution and the simulated data generated in TAPsolver, while the right panel of Figure <ref> shows the comparison between the analytical and simulated irreversible adsorption processes. The differences between the analytical and simulated data in each of these examples are negligible, confirming the accuracy of the TAPsolver simulation. 
  
  
  
  
  
    §.§.§ Validation of the Inverse Problem
  
  To test the accuracy of the adjoint-based optimization routine implemented in TAPsolver, synthetic carbon monoxide oxidation is simulated and analyzed. The simulated mechanism included reversible carbon monoxide adsorption, irreversible dissociative O_2 adsorption, and both Eley-Rideal and Langmuir-Hinshelwood reactions that lead to the formation of carbon dioxide. These elementary reactions and their values are presented in Table <ref>, as well as the initial guesses and the converged values following optimization.
  
  The initial values of the kinetic parameters used for optimization are all 1.0*10^-10. These values are intentionally far from the true values to illustrate the robustness of the optimization routine even when no prior knowledge of the kinetic parameters is available. Even though elementary reactions 2 and 4 were considered to be irreversible, reverse rate constants were considered during optimization. Exact matches between the fitted kinetic parameters and the synthetic values are not obtained, but strong agreement is observed with the reverse rate constants being notable exceptions. The convergence of the optimization problem is presented in Figure <ref>. With the initial guesses used, no reaction is initially observed and the reactants act as if they are purely diffusive. Gradually, the kinetic parameters are altered such that both the reactant and product curves match the synthetic data. 
  This result validates that TAPSolver is able to successfully fit a curve of synthetic data, and that the primary parameters that govern the kinetics are correctly identified, even with poor initial guesses. However, not all kinetic parameters are accurately identified since the objective function is not sensitive to all parameters. This highlights an inherent limitation of fitting TAP pulses, since the data from a single pulse will not be sensitive to all parameters.
  
  
   
  -5pt
  |c|c|c|c|c|c|
  List of the elementary reaction values used in the generation of the synthetic carbon monoxide oxidation data. Dashed values in the 'Actual Value' column indicate that the parameter was not included in the simulation and is therefore effectively zero.  
  
  
  Kinetic Parameter     Reaction      Actual Value       Initial Guess     Converged Value     Units
  
  1f     CO + * → CO^*     1.5*10^0       1*10^-10     1.48*10^0     cm^3/nmol s 
   
  1b     CO^*→ CO + *      1.5*10^-1       1*10^-10     3.25*10^-5     1/s 
  
  2f     O_2 + 2* → 2O^*      5.00*10^-3      1*10^-10     4.70*10^-3     cm^6/nmol^2 s 
  
  2b     2O^*→ O_2 + 2*      –      1*10^-10     1.00*10^-10     cm^3/nmol s 
   
  3f     CO^* + O* ↔ CO_2 + 2*      1.05*10^1       1*10^-10     1.04*10^1     cm^3/nmol s 
  
  3b     CO_2 + 2* ↔ CO^* + O*      1.5*10^-2      1*10^-10     6.10*10^-3     cm^6/nmol^2 s  
  
  4f     CO + O* → CO_2 + *      2.02*10^1       1*10^-10     2.50*10^1     cm^3/nmol s
  
  4b     CO_2 + * → CO + O*      –       1*10^-10     5.28*10^-3     cm^3/nmol s 
   
  
  
  
  
  
  
  
  
   §.§ TAPsolver Efficiency and Precision
  
  
  The influence of the mesh size was examined and is shown on the left in Figure <ref>. A uniform mesh can be utilized, where the inert and catalyst zones have identical mesh distributions. There is a linear increase in simulation time with an in increase in mesh size.
  
  As mentioned in Section <ref>, an alternative method of defining the mesh has been implemented in TAPsolver, where increases in the precision of the catalyst zone can be made without refinement of the inert zones. The timing for a refined mesh is also shown in Figure  <ref> alongside the equivalent uniform mesh. Refinement in this context is defined as the number of times the number of cells in the catalyst zone are doubled. If the number of cells is initially four, a refinement of two will lead to a catalyst mesh size of sixteen, while a refinement of five will lead to a catalyst mesh size of one-hundred and twenty-eight.
  The reduction in simulation time, even at smaller mesh sizes, by targeting the region of refinement is significant and can lead to the efficient and accurate simulation of reactors when a strong concentration profile within the catalyst is observed.
  In the current implementation of TAPsolver, increasing the number of time steps leads to a linear increase in simulation time.
  An explicit mixed Euler and Crank-Nicolson time stepping method is being utilized. However, future versions will explore implicit time stepping strategies  <cit.> that may be substantially more efficient, particularly in the case of stiff problems.
  
  
  
  The time requirements for evaluating the gradient and Hessian are shown on the right of Figure <ref>. Once again, simulated carbon monoxide data is used as a test. A linear increase in simulation time is observed in both figures, although the slope is substantially lower for the gradient calculation. This is a considerable improvement over numerical differentiation, where the time needed to compute derivatives scales quadratically with the number of parameters.
  In this example, the entire Hessian is calculated explicitly without considering symmetry, and therefore should be faster if used in optimization routines. Similar to inefficiencies in the forward problem, improvements to the efficiency of the simulator in TAPsolver will help make adjoint calculations less expensive, since the construction and evaluation of adjoints is directly dependent on the number of time steps.
  
  The precision of the algorithmic differentiation in TAPsolver is also explored through comparisons between derivatives computed using the finite difference (FD) and algorithmic differentiation approaches. The central difference approximation <cit.> of first-order derivatives is defined as 
  
      ∂ f/∂ k_i = f(k_i + h) - f(k_i - h)/2h
  
  
  where f is an arbitrary function of parameters k_i and h is a step size.
  
  The step size is selected to be a fraction of the parameter value. To test the convergence of these values, multiple step sizes were selected, including 1/50, 1/500 and 1/5000. Comparisons between the first order derivatives calculated by the central difference and adjoint approaches at the initial guess values of each parameter are presented in Tables <ref> and <ref>.
  
  
  
  -5pt
  |c|c|c|c|c|
  Values of the first order derivatives found through FD approximations  and the adjoint approach at the initial guess. Multiple step sizes were used for FD to show the convergence of the solution.  
  
  Kinetic  Parameter      1/50     1/500     1/5000     A.D.
  
  1f     -8.71*10^1     0.0     0.0
       -8.71*10^1 
  
  1b     0.0      0.0     0.0
       1.04*10^-9 
  
  2f      -2.28*10^4     -2.28*10^4     0.0     -2.28*10^4 
   
  2b      0.0      0.0     0.0     5.63*10^-13 
   
  3f      0.0     0.0     0.0     -3.32*10^-17 
   
  3b     0.0     0.0     0.0     -1.50*10^-15 
   
  4f      0.0      0.0     0.0     -1.03*10^-15 
   
  4b     0.0     0.0     0.0     2.00*10^-15 
   
  
  
  
  
  
  
  -5pt
  |c|c|c|c|c|
  Values of the first order derivatives found through FD approximations  and the adjoint approach at the converged value. Multiple step sizes were used for FD to show the convergence of the solution.  
  
  Kinetic  Parameter      1/50     1/500     1/5000     A.D.
  
  1f     2.82*10^0     2.62*10^0     2.60*10^0
       2.60*10^0 
    
  1b     3.10*10^-2      -4.43*10^-2     0
       -2.88*10^-2 
   
  2f      -2.20*10^3     5.16*10^1     4.63*10^1     4.57*10^1 
   
  2b      4.83*10^-1      4.83*10^-1     4.83*10^-1     4.83*10^-1 
   
  3f      7.67*10^-2     7.54*10^-2     7.52*10^-2     7.52*10^-2 
   
  3b     -1.76*10^0     -1.76*10^0     -1.76*10^0     -1.76*10^0 
   
  4f      -1.74*10^-3      -1.87*10^-3     -1.89*10^-3     -1.89*10^-3 
   
  4b     -1.76*10^-1     -1.76*10^-1     -1.76*10^-1     -1.76*10^-1 
   
  
  
  
  
  
  The results in Table <ref> indicate that when the initial guess is far from the optimum the FD approximation routinely fails, yielding derivatives of zero for many parameters. These parameters cannot be optimized through a finite difference approach unless an improved initial guess is provided. The derivatives at the final converged solution (Table <ref>) indicate that the finite difference approach is significantly more accurate near the optimum, although some small discrepancies are still observed. For example, the sign of the derivative for parameter 1b changes from positive to negative depending on the step size used, highlighting a lack of stability in the finite difference approach.
  The implication of this finding is that TAPsolver will be much more robust to weak initial guesses, enabling the use of less-biased initial guesses (e.g. setting all rate constants a constant value) or values from DFT that may capture trends in the rate constants but have absolute values that are off by many orders of magnitude <cit.>.
  In addition, the computational time required to compute algorithmic derivatives scales linearly with the number of rate constants (although the prefactor is large), suggesting that the approach will be more efficient than the FD approximation for large reaction networks. 
  
  
  
   §.§ Application: analyzing experimental carbon monoxide oxidation data
  
  
  Carbon monoxide oxidation on platinum catalysts has frequently been studied over the last decades <cit.>, including variations to reactant partial pressures <cit.>, cluster sizes <cit.>, surface cleavages <cit.> and support materials <cit.>. In this case study, a microkinetic model consisting of both Langmuir-Hinshelwood and Eley-Rideal mechanisms, with elementary steps shown in Table <ref> is used to fit the experimental TAP response. Thermodynamic constraints are applied to ensure that the rate constants obey the overall thermochemistry of the reaction. The resulting fit accurately captures the behavior of the TAP pulse response (Figure <ref>), and the resulting forward and reverse rate constants for each step are provided in Table <ref>. 
  
  
  
  This case study illustrates the power of the TAPsolver approach, yielding a full set of thermodynamically-consistent kinetic parameters from a single TAP pulse, measured in around 1 second. However, it is worth recalling that not all parameters will be well-determined, since even in the case of fitting to synthetic data it was not possible to accurately determine all kinetic parameters. Moreover, some slight discrepancies can be seen between the model and experimental data, particularly at short and long timescales for CO_2. This may be due to limitations in the microkinetic model used, such as neglecting adsorbate-adsorbate interactions or possible missing reactions such as oxygen transfer to the bulk. Addressing these challenges will require further development to the TAPsolver framework, such as integration with global optimization routines <cit.>, inclusion of additional elementary steps and lateral adsorbate interactions, or application of model reduction criteria<cit.>. Nonetheless, this brief case study illustrates TAPsolver's ability to fit even noisy data sets, and highlights the ability to directly extract intrinsic kinetic parameters from TAP data.
  
  
  
  -5pt
  |c|c|c|c|c|c|c| 
  List of kinetic parameters observed while fitting the Langmuir-Hinshelwood (L.H.) and Eley-Rideal (E.R.) reaction mechanisms with and without thermodynamic consistency constraints.  
  
      Reaction      L.H.       E.R.     L.H. Thermo.     E.R. Thermo.     Units 
   
  1f     CO + * → CO^*     3.88*10^1     7.94*10^1     3.94*10^1     4.3*10^1     cm^3/nmol s 
  
  1b     CO^*→ CO + *     1.24*10^1     2.89*10^1     6.72*10^0     1.69*10^1     1/s 
   
  2f     O_2 + * → O_2^*     1.78*10^0     1.59*10^0     5.44*10^-2     9.73*10^-1     cm^3/nmol s 
   
  2b     O_2^*→ O_2 + *     5.53*10^0     1.04*10^1     1.91*10^-1     8.79*10^1     1/s 
   
  3f     O_2^* + * → 2O*     1.03*10^0     1.26*10^0     1.08*10^1     8.68*10^-2     cm^3/nmol s 
   
  3b     2O*→ O_2^* + *     5.21*10^-2     2.46*10^-2     5.23*10^-21     8.06*10^-4     cm^3/nmol s 
   
  4f     CO^* + O^*→ CO_2 + 2*     7.16*10^-1     –     1.48*10^1     –     cm^3/nmol s
   
  4b     CO_2 + 2* → CO^* + O^*     3.19*10^-1     –     5.62*10^-20     –     cm^3/nmol s 
   
  5f     CO + O* ↔ CO_2^* + *     –     6.76*10^0     –     5.20*10^0     cm^3/nmol s
   
  5b     CO_2^* + * ↔ CO + O*     –     3.59*10^-1     –     1.68*10^0     cm^6/nmol^2 s 
  
  
  
  
  
  
  
  § CONCLUSION
  
  The analysis of TAP reactor data, an often arduous process, has been streamlined through the development of a Python program called TAPsolver that features novel analysis methods, including algorithmic differentiation and thermodynamic consistency optimization constraints. Algorithmic differentiation provides the highest accuracy derivatives and can reduce the cost of gradient estimation, particularly when the number of parameters is large. Applying constraints to the thermodynamics can reduce the number of potential reaction mechanisms and provide physically meaningful results. These improvements make the top-down identification of microkinetic models from transient kinetic data more feasible. The forward problem has been validated through comparisons with analytical solutions and the efficiency/accuracy of the current implementation of algorithmic differentiation is demonstrated. Similarly, the inverse problem was validated by fitting noiseless, simulated carbon monoxide oxidation data. Finally, the work flow was applied to experimental carbon monoxide oxidation data on a platinum catalyst. 
  
  The TAPsolver code provides a flexible new tool for TAP data analysis that can easily be integrated with emerging techniques in machine learning and optimization through a convenient Python interface. The tool has been developed in collaboration with experimentalists, facilitating easy transfer of parameters and data between experiment and simulation, and the code utilized input files that make it easy to explore different reaction mechanisms or change assumptions about reaction kinetics. The use of algorithmic differentiation enables accurate derivatives even far from the optimum, enabling convergence even with weak initial guesses, and gradient calculations scale linearly with the number of parameters. These advantages make TAPsolver particularly well-suited for problems with complex reaction mechanisms. We hope that TAPsolver provides a convenient platform for implementation of new TAP data analysis techniques, and lowers the barrier to TAP data analysis for both experimental and computational researchers.
  
  
  
  § ACKNOWLEDGEMENTS
  
  Support for this work was provided by the U.S. Department of Energy (USDOE), Office of Energy Efficiency and Renewable Energy (EERE), Advanced Manufacturing Office Next Generation R&D Projects under contract no. DE-AC07-05ID14517. The authors are also grateful to Dr. Christian Reece for discussions of simulation methods and input file construction.
  
  
  
  § CONFLICT OF INTEREST
  
  The authors declare no competing financial interest.
  
  ieeetr
  
  
  