DOI: 10.1039/C7RE00040E
(Paper)
React. Chem. Eng., 2017, Advance Article

Kuijun Li^{a},
Priyadarshi Mahapatra^{b},
K. Sham Bhat†^{c},
David C. Miller^{d} and
David S. Mebane*^{a}
^{a}Department of Mechanical and Aerospace Engineering, West Virginia University, Morgantown, WV, USA. E-mail: david.mebane@mail.wvu.edu
^{b}National Energy Technology Laboratory, Morgantown, WV, USA
^{c}Los Alamos National Laboratory, Los Alamos, NM, USA
^{d}National Energy Technology Laboratory, Pittsburgh, PA, USA

Received
22nd March 2017
, Accepted 9th June 2017

First published on 9th June 2017

Inaccuracy in chemical kinetic models is a problem that affects the reliability of large-scale process models. Detailed and accurate kinetic models could improve this situation, but such models are often too computationally intensive to utilize at process scales. This work presents a method for bridging the kinetic and process scales while incorporating experimental data in a Bayesian framework. “Dynamic Discrepancy Reduced Modeling” (DDRM) enables the use of complex chemical kinetic models at the process scale. DDRM inserts Gaussian process (GP) stochastic functions into the first-principles dynamic model form, enabling a sharp reduction in model order. Uncertainty introduced through the order reduction is quantified through Bayesian calibration to bench-scale experimental data. The use of stochastic functions within the dynamic chemical kinetic model enables this uncertainty to be correctly propagated to the process scale. The model reduction and uncertainty quantification framework was applied to a reaction–diffusion model of mesoporous silica-supported, amine impregnated sorbents; model order reduction of over an order of magnitude was achieved through a sharp reduction in the resolution of the approximate solution for diffusion. GPs of the Bayesian smoothing splines analysis of variance (BSS-ANOVA) variety were applied to the diffusion coefficients and equilibrium constants in the model. The reduced model was calibrated using experimental dynamic thermogravimetric data, with prior probability distributions for physical model parameters derived from quantum chemical calculations. Model and parameter uncertainties represented in the posterior distribution were propagated to a bubbling fluidized-bed adsorber simulation implemented in Aspen Custom Modeler. This work is the first significant demonstration of the dynamic discrepancy technique to produce robust reduced order models for bench-to-process multi-scale modeling and serves as a proof-of-concept.

Carbon capture and storage (CCS) technology is undergoing serious consideration as a method to reduce the emission of CO_{2} from power plants.^{4–9} However, the traditional timeline to advance new technologies from laboratory to industrial-scale deployment in the power industry takes twenty to thirty years.^{10} Part of the reason for this is that the results of bench and laboratory-scale experiments are difficult to utilize for predicting performance of large-scale processes.

Advances in multi-scale modeling have great potential to significantly decrease the time needed for deployment of carbon capture technologies.^{11,12} The U.S. Department of Energy's Carbon Capture Simulation Initiative (CCSI) focuses on the development of multi-scale models and computational tools for accelerating the development and scale-up of new carbon capture technologies for large-scale power plants; the project involves five national laboratories, five universities and more than twenty industrial partners.^{13,14} The quantification of uncertainty in chemical kinetic models and its propagation to larger-scale simulations has been an important part of the project.^{10}

Mesoporous silica-supported, amine-impregnated (SSA) solid sorbents were one of the CO_{2} capture materials chosen to demonstrate the capabilities of the CCSI Toolset. These sorbents present a complex set of behaviors, especially with respect to the influence of H_{2}O on the diffusion rate and effective capacity of the sorbent for CO_{2}.^{15,16} A physical model, consistent with the known experimental behavior and quantum chemical calculations, has been developed.^{15–17} This thermodynamic and kinetic reaction–diffusion model, which considers a number of diffusive intermediates, is considerably more complex than chemical kinetic models typically employed in device models for large-scale process simulations. The challenge is to reduce the computational burden of the kinetic model to the point that it can be used in a process-scale simulation, while retaining the basic physics of the model and providing a vehicle for uncertainty quantification.

This challenge was met using dynamic discrepancy reduced modeling (DDRM), a new technique which is introduced in this article. DDRM is a methodology that reduces the order of dynamic systems by eliminating intermediate states (those that are not end points in the dynamic network) from the system. The model variability that is lost through the order reduction is replaced by incorporating Gaussian process (GP) stochastic functions into the model. The reduction is accomplished in such a manner that important features of the original dynamic system are retained. The GPs are discrepancy functions in that they contain the error of the reduced order model with respect to the original, and also potentially the error of the original model form with respect to reality. DDRM as such is a mixture of first-principles and empirical order reduction; DDRM models incorporate experimental results (or, one could imagine, the results of high-fidelity models) through the use of discrepancy functions, while retaining the important features of first-principles model forms.

DDRM is a technique that is based on the “dynamic discrepancy” method for model-based Bayesian calibration in dynamic systems.^{18} The key to the dynamic discrepancy technique is the intrusive use of stochastic discrepancy functions within the dynamic system itself. Used intrusively, discrepancy functions depend, in a simple fashion, on the states of the system. By contrast, model discrepancy functions applied outside of the dynamic model must depend on the Hilbert space of time-dependent functions of the system's inputs (such as temperature and concentrations of reactants and products) in order to properly account for extrapolations through that space; this is, generally speaking, completely impractical. Because propagation of uncertainty in a chemical kinetic model calibrated to bench-scale experiments to the process-scale almost inescapably involves such extrapolation, the dynamic discrepancy approach is essential for rigorous bench-to-process scale-bridging.

The stochastic dynamic systems resulting from the dynamic discrepancy technique are integrated during the process of model calibration to experimental datasets using Markov chain Monte Carlo (MCMC) sampling routines. Crucial to this capability is the nature of the GPs used for the discrepancy functions, which are of the Bayesian smoothing spline analysis of variance (BSS-ANOVA) variety.^{19} BSS-ANOVA GPs can be decomposed in such a manner that the stochasticity in the function resides wholly in a set of correlated parameters, the posterior distribution of which are estimated through the MCMC process.

The methodology of Bayesian calibration of computer models to data using discrepancy functions was pioneered by Kennedy and O'Hagan.^{20} This framework has been used for calibration and uncertainty quantification in many different areas, such as forest models,^{21,22} thermal problems,^{23,24} biogeochemical models,^{25} water quality,^{26,27} hydrological models,^{28,29} computational fluid dynamics models,^{30,31} and carbon capture.^{18,32–35} Work at Georgia Institute of Technology has used Bayesian calibration to quantify both epistemic and aleatory uncertainty in amine-based CO_{2} sorbents.^{36,37}

(1) |

(2) |

(3) |

(4) |

(5) |

The diffusive intermediates – amine-stabilized zwitterions, adsorbed water, and water-stabilized zwitterions – are denoted z_{1}, z_{2}, and z_{3}, respectively. Carbamates – ammonium carbamate and hydronium carbamate – are denoted x and y, respectively. Assuming physisorption reactions at the surface are in equilibrium, the site fractions of diffusive intermediates at the surface are

z_{1s} = κ_{1}P_{CO2}S_{s}^{2}
| (6) |

z_{2s} = κ_{3}P_{H2O}S_{s}
| (7) |

z_{3s} = κ_{4}P_{CO2}z_{2s}
| (8) |

The rate equations for the formation of ammonium carbamates and hydronium carbamates are

(9) |

(10) |

(11) |

(12) |

(13) |

The total number of model parameters including microstructural parameters, material properties, thermodynamic parameters, kinetic parameters and transport parameters is 23: θ = [θ_{1},…,θ_{23}]^{T}. A summary of input, output and model parameters appears in Table 1.

Functional input profile (potential system conditions) ζ(t) |
---|

ξ_{1}: temperature T |

ξ_{2}: partial pressure of CO_{2} |

ξ_{3}: partial pressure of H_{2}O |

Experimental outputs and state variables x and y |
---|

x: weight fraction of ammonium carbamate |

y: weight fraction of hydronium carbamate |

(14) |

Kennedy and O'Hagan introduced a method for model calibration that includes a quantification of error in the model itself.^{20} In the Kennedy–O'Hagan approach, the statistical model for the data is

(15) |

(16) |

Ṁ = F(M,δ;θ) | (17) |

The stochastic dynamic system can be solved and the joint posterior for (θ,δ) can be obtained using MCMC sampling if a judicious choice is made for the form of the discrepancy GP.

The BSS-ANOVA GP framework was introduced by Reich et al.^{19} The discrepancy function δ can be formulated using the BSS-ANOVA GP model^{19} and provides two critical advantages for the model (16). Most importantly, it relegates stochasticity in the GP to a set of coefficients β, which enables the deterministic evaluation of a likelihood function based on (16) when β is fixed. It also improves the computational efficiency by scaling linearly with the number of data points, as opposed to (N^{3}) for a traditional GP.^{18,38}

The discrepancy function can be subjected to an ANOVA decomposition:

(18) |

(19) |

K_{2}((u,v),(u′,v′)) = K_{1}(u,u′)K_{1}(v,v′)
| (20) |

Each functional component in (18) can be further written as an orthogonal basis expansion following a Karhunen–Loéve decomposition:^{19}

(21) |

(22) |

(23) |

Fig. 1 First nine eigenfunctions from the Karhunen–Loéve expansion for a main effect function from the BSS-ANOVA covariance. |

3.3.1 Interactions. The model as presented in sec. 2 above assumes ideal thermodynamics and kinetics; that is, it does not consider any interactions (aside from entropic interactions related to site occupancy) among the various chemical species appearing in the model.^{16} However, it is clear from the thermogravimetric analysis (TGA) experiments reported in this work that some non-idealities must be considered. One way to introduce such interactions in a continuum-level model is to explicitly consider associates among species; however, this tends to increase the complexity of the model. Discrepancy functions incorporated into the model can also treat interactions in an effective sense; for instance, the influence of interactions on the thermodynamics can be captured by applying discrepancy functions to the equilibrium constants:

such that the new, interaction-sensitive quantity κ^{E} depends now on the relevant thermodynamic quantities of partial pressure p and the temperature T. For the current study, for each of the reactions appearing in the model:‡

(24) |

(25) |

(26) |

(27) |

(28) |

(29) |

3.3.2 Diffusion. Discrepancy enables a direct reduction of the number of states resulting from the discretization of the diffusion model: discrepancy functions can be added to the transport coefficients in analogous fashion as for the equilibrium constants to replace the variability lost through the order reduction. A one-dimensional finite volume scheme consists of a set of control volumes (CVs) with fluxes into and out of each CV. In the reduced order model, discrepancy functions will be applied to each of these fluxes, which are specified at the cell boundary, as

where z_{i,j} is the site fraction of diffusive intermediate i in CV j. Discrepancies depend only on the concentration difference across the cell boundary. Because the constitutive equation for transport is the same throughout the material, the discrepancy at each CV boundary will be the same. A typical converged resolution for the model is 100 CVs; a reduction to five CVs using discrepancy thus constitutes a 20-fold decrease in the model order. This strategy was applied to each of the diffusive species:

(30) |

(31) |

(32) |

(33) |

A study showed that a balance between model variability and complexity was achieved at eight basis functions per main effect and sixteen basis functions per two-way interaction, for a total of 32 terms per discrepancy, or 256 discrepancy parameters overall. These are calibrated at the same time as the 23 physical parameters listed in Table 1, for a total of 279 parameters in the calibration overall. The reduction in model order decreased the evaluation time for the model from several minutes to less than 5 seconds, thereby enabling the calibration and uncertainty propagation reported below.

3.4.1 Monte Carlo sampling. Because of the intractability of the integral in the denominator on the right-hand side of eqn (14), a numerical simulation is required to obtain samples from the posterior distribution. MCMC sampling^{42,43} was used to obtain the posterior distributions of θ, δ and σ^{2}. In each iteration of the routine, a set of model parameters θ and model discrepancy parameters β was proposed, the likelihood was evaluated, and the proposal was either accepted or rejected using the Metropolis–Hastings (MH) criterion.^{42} Taking advantage of an inverse gamma prior that is conjugate to the white noise model for the likelihood, Gibbs updates^{44} were used to obtain samples for the observation error variance σ^{2}.

Block proposals for the parameters subject to MH sampling greatly decreased the time required for a single Monte Carlo step, which occurs when a new proposal for each parameter has been evaluated. Block proposals for the physical parameters were made using an adaptive, multivariate normal proposal distribution.^{45} A random vector was drawn from the multivariate normal distribution with a mean vector = [_{1},...,_{23}] fixed at the current point in the random walk, and covariance matrix V is the empirical covariance of a subset of the previously accepted values in the MCMC chain:

(34) |

The covariance matrix was initialized and recalculated every 1000 MCMC steps based on the last 2000 samples. It then was used as the covariance matrix for the next 1000 proposals. An adjustable multiplier on the covariance guides a reasonable acceptance rate of the proposal. After burn-in, the covariance matrix is fixed at the final proposal covariance.

For δ, all the β coefficients for each main effect and second order interaction coefficients (e.g. eqn (21)) were updated simultaneously using an adaptive multivariate normal proposal as described above. After a sufficient number of MCMC iterations, the samples of θ and β (and other parameters) will converge to its stationary distribution; this is the posterior distribution of the model parameters, model discrepancy parameters, and observation error parameters.

3.4.2 Prior distribution. To complete the model specification, prior distributions are required for θ, δ, and σ^{2}. Priors for some of the model parameters θ are derived from ab initio quantum chemistry calculations or from scientific literature and are specified using a normal distribution. The details of the quantum chemistry methods used to obtain these prior can be found in ref. 16. The model for δ in (22) requires a prior specification for τ_{j}'s from (23). A diffuse inverse gamma prior is chosen for the τ_{j}'s to promote better mixing of β and an inverse gamma prior is chosen as well for σ^{2}. The selected priors for the τ_{j}'s and σ^{2} are conjugate to the normal distribution of β's and the likelihood respectively, which ensure Gibbs sampling for for the τ_{j}'s and σ^{2} during the MCMC procedure.

Fig. 2 Schematic of the bubbling fluidized bed model.^{46} |

In this study, we assume the particles are of Geldart group A. Since the process model utilized in this study uses a more complex reaction–diffusion particle model which introduces a second spatial dimension, as described in sections 2–3 above, the original BFB model has been simplified relative to the model described in ref. 46. The gas velocity is well above the minimum fluidization velocity and within the “bubbling fluidization regime”; however, it is assumed that no bubbles are formed. The model utilizes instead a single “emulsion” region/phase where gas is direct contact with the solid particle. The particle length scale is radially discretized (spherical domain) into 30 control volumes. The device length scale is axially discretized (cylindrical domain) into 25 control volumes. Both solid sorbents and flue gas consisting of CO_{2}, H_{2}O, and N_{2} are injected at the bottom of the bed and flow upward through the absorber as the particles are transported axially through the bed, removing CO_{2} from the gas.

The model is implemented in Aspen Custom Modeler (ACM) (Aspen Technologies, Inc.). The particle length scale was solved using ACM's partial differential equation (PDE) solver utilizing a 2nd-order central finite difference method. At the process-scale, a compartment-based approach, where axial differential terms were expressed as finite differences involving a first-order forward finite difference method, was used. The physical properties of the gas phase were calculated with the cubic equation of state. The physical properties of the solids were based on experimental measurements. The process model involves 332890 equations and a significantly larger number of non-zeroes, as reported by ACM. Computations have been conducted on a Windows-based Workstation utilizing an 8-core Intel i7 Processor.

The discrepancy terms were evaluated using external procedure calls in ACM. These procedures were coded in C++ and compiled into a dynamic link library for fast evaluation of discrepancy terms and their residuals within the process models. The dynamic discrepancy evaluation framework was automated, with error handling capability, to read each posterior point from an excel spreadsheet and report the corresponding Lagrangian profiles for the solid particles. Due to the large magnitudes of the “step-changes” involved and the presence of external procedure calls, a homotopy-based approach was used.

The process model involves a steady-state evaluation with process conditions varying spatially. However, to analyze the data corresponding to the discrepancy at the particle scale, a conversion to a temporal domain is required. This involves an interpretation of process conditions surrounding individual particles as they axially traverse through the adsorber bed. Hence, a Eulerian to Lagrangian transformation,^{47} also known as a “material” or “substantive” derivative, is utilized based on the flow velocity of the solid particles in the emulsion region.

Each sample drawn from the posterior distribution involves a simulation of the absorber column within ACM, and the operating conditions and overall CO_{2} capture in the device are reported.

Coverage of the experimental data is shown by reproducing and overlaying a set of model realizations drawn from the posterior against the dynamic TGA data. Fig. 4(b)–(h) shows the coverage of calibration results to the data under dry CO_{2}, H_{2}O and humidified CO_{2} conditions. The green curve shows the temperature change during the adsorption process. The blue curve is the TGA trace of the weight fraction of adsorbed CO_{2}, and red curves are 100 model realizations drawn from the posterior. TGA data was well covered by calibration results in these seven cases. Relatively small regions of coverage gaps can mostly be attributed to systematic experimental error not considered in the model.

Bivariate scatter plots visualize the correlation between pairs of model parameters. Selected bivariate scatter plots appear in Fig. 5. In the figures, yellow represents regions of high sample density, while blue indicates relatively lower density. Strong correlations appear between certain pairs of model parameters. The correlations between ΔH_{1} and ΔS_{1}, ΔH_{3} and ΔS_{3}, and ΔH_{4} and ΔS_{4} were 0.98, 0.91 and 0.79, respectively. Generally speaking, the correlation between the formation enthalpy and entropy reflects the functional form of the equilibrium constant. The correlations between and log_{10}(ς_{b1}), and log_{10}(ς_{b2}), and log_{10}(ς_{b3}) were 0.98, 0.99 and 0.5, respectively; similarly, these correlations are expected due to the functional form of the mobility.

5.2.1 Adiabatic. The histogram and kernel smoothing density fit of carbon capture fraction from the dry gas and humidified gas adsorbers in the adiabatic case are shown in Fig. 6. In the dry case – Fig. 6(a) – the mean carbon capture fraction is 0.54 and the 95% credible region for the carbon capture fraction is between 0.14 and 0.94. In the humidified case – Fig. 6(b) – the mean carbon capture fraction is 0.56 and the 95% credible region is between 0.24 and 0.88. In both cases the significant uncertainty in the capture rate reflects the paucity of kinetic data appearing in the dynamic TGA used for calibration. The small difference in the mean capture rate reflects the fact that water is barely adsorbed when the temperature is higher than 75 °C, which is quickly exceeded in the adiabatic case (reaching in some cases temperatures as high as 85 °C) as particles move through the bed. Forty realizations from the posterior showing the operating temperature and CO_{2} partial pressure in the bed for the adiabatic case are shown in Fig. 7 and 8, respectively.§

Fig. 7 Temperature in the adsorber under dry and humidified CO_{2} conditions for the adiabatic case. |

5.2.2 External cooling. In order to observe a larger influence of water on the process, external cooling (water at 20 °C, 209 kg s^{−1}) was added to internal heat transfer tubes within the BFB model. The corresponding influence on the adsorber operating temperature is shown in Fig. 9: the significant suppression of the reactor temperature means that more water from the flue gas will adsorb on the sorbent, with attendant increase of the effective CO_{2} capacity according to the kinetic model. In this case the differences in the temperature profiles between dry and humidified cases are much more pronounced. These differences are also reflected in the CO_{2} capture fraction distributions with external cooling, as shown in Fig. 10. In the dry case – Fig. 10(a) – the mean carbon capture fraction is 0.56 (only slightly higher than the case without cooling) and the 95% credible region for the carbon capture fraction is 0.28 and 0.84. In the humidified case – Fig. 10(b) – the mean carbon capture fraction is 0.73 and the 95% credible region is between 0.47 and 0.99, reflecting a considerable increase over the case without cooling. The CO_{2} partial pressure profiles in the adsorber for both dry and humidified cases are shown in Fig. 11.

Fig. 10 Distributions of carbon capture rate under dry and humidified CO_{2} conditions with external cooling. |

The stochastic nature of the predictions is another key advantage of the DDRM method and the associated uncertainty-based paradigm for multi-scale modeling. The process-scale simulations force extrapolation in the kinetic model, as revealed by the relatively large uncertainty in the process results compared to the tightness of the fits to the dynamic TGA data. The uncertainty in the results can provide a more reliable basis for process design and scale up, in particular when combined with strategies for optimization under uncertainty.

Physically speaking, discrepancy functions represent interactions among chemical species in a model which, without discrepancy, ignores such interactions. In principle, the discrepancy functions should thus depend on all chemical states in the system plus the temperature: for the equilibrium constants, it would suffice to make the discrepancy functions depend on all partial pressures in the system; for the transport coefficients, all of the chemical species within the solid. Such an approach would be comprehensive; however, it would also increase the number of terms for calibration. As such this study limited the dependency of each discrepancy function to a single chemical constituent plus the temperature, which proved to be sufficient for this study. Future work will explore methods for generalizing the dependence of the discrepancy functions while limiting the number of variables for calibration.

In this work, the discrepancy enables fully quantitative correspondence between the model and the data. The physical model provides a semi-quantitative correspondence to the same dataset, as seen in ref. 16. The stochastic functions thus provide the increase in variability required for fully quantitative correspondence between model and experiment (as seen in this paper), along with model reduction that enables calibration and propagation of uncertainty.

In the present study, the predictions at the process-scale reflect a significant amount of uncertainty in the kinetics. This can be attributed to the relative dearth of kinetic information in the dynamic TGA experiment, thus motivating alternative bench-scale experiments (such as fixed bed) that may yield more useful data. The stochastic process-scale results present opportunities for the design of such experiments due to the dynamic results for bed temperature and gas composition: each curve depicted in Fig. 9 and 11, for instance, represents a dynamically changing set of conditions that bench-scale experiments could attempt to replicate. Adding a sampling of conditions drawn from these figures to the set of experiments used for model calibration would very likely decrease the uncertainty in the process-scale predictions. Integration of the method with bench-scale experimentation, flowsheet modeling and optimization under uncertainty presents significant opportunities for decreasing the time and costs of scale-up through more accurate modeling of industrial scale systems.

- J. T. Houghton, Climate change 1995: The science of climate change: contribution of working group I to the second assessment report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 1996 Search PubMed.
- J. H. Butler and S. Montzka, The NOAA annual greenhouse gas index, Noaa earth system research laboratory technical report, 2011 Search PubMed.
- M. K. Mondal, H. K. Balsora and P. Varshney, Energy, 2012, 46, 431–441 CrossRef CAS.
- A. B. Rao and E. S. Rubin, Environ. Sci. Technol., 2002, 36, 4467–4475 CrossRef CAS PubMed.
- J. Davison, Energy, 2007, 32, 1163–1176 CrossRef CAS.
- A. Samanta, A. Zhao, G. K. Shimizu, P. Sarkar and R. Gupta, Ind. Eng. Chem. Res., 2011, 51, 1438–1463 CrossRef.
- A. S. Bhown and B. C. Freeman, Environ. Sci. Technol., 2011, 45, 8624–8632 CrossRef CAS PubMed.
- M. Wang, A. Lawal, P. Stephenson, J. Sidders and C. Ramshaw, Chem. Eng. Res. Des., 2011, 89, 1609–1624 CrossRef CAS.
- D. C. Miller, J. T. Litynski, L. A. Brickett and B. D. Morreale, AIChE J., 2016, 62, 2–10 CrossRef CAS.
- D. C. Miller, M. Syamlal, D. S. Mebane, C. Storlie, D. Bhattacharyya, N. V. Sahinidis, D. Agarwal, C. Tong, S. E. Zitney, A. Sarkar, X. Sun, S. Sundaresan, E. Ryan, D. Engel and C. Dale, Annu. Rev. Chem. Biomol. Eng., 2014, 5, 301–323 CrossRef CAS PubMed.
- H. Simon, T. Zacharia and R. Stevens, et al., Modeling and Simulation at the Exascale for Energy and the Environment, Department of Energy Technical Report, 2007 Search PubMed.
- X. Jiang, Appl. Energy, 2011, 88, 3557–3566 CrossRef CAS.
- D. C. Miller, D. A. Agarwal, X. Sun and C. Tong, “CCSI and the role of advanced computing in accelerating the commercial deployment of carbon capture systems,” Proceedings of the SciDAC 2011 Conference, 2011 Search PubMed.
- D. C. Miller, M. Syamlal, J. Meza, D. Brown, M. Fox, M. Khaleel, R. Cottrell, J. Kress, X. Sun and S. Sundaresan, et al., “Overview of the US DOE's Carbon Capture Simulation Initiative for accelerating the commercialization of CCS technology,” 36th International Technical Conference on Clean Coal & Fuel Systems, 2011 Search PubMed.
- D. S. Mebane, J. D. Kress, C. B. Storlie, D. J. Fauth, M. L. Gray and K. Li, J. Phys. Chem. C, 2013, 117, 26617–26627 CAS.
- K. Li, J. D. Kress and D. S. Mebane, J. Phys. Chem. C, 2016, 120, 23683–23691 CAS.
- A. Lee, D. Mebane, D. Fauth and D. Miller, “A model for the adsorption kinetics of CO
_{2}on amine-impregnated mesoporous sorbents in the presence of water,” 28th International Pittsburgh Coal Conference, Pittsburgh, PA, 2011 Search PubMed. - K. S. Bhat, D. S. Mebane, C. B. Storlie and P. Mahapatra, J. Am. Stat. Assoc., 2017 DOI:10.1080/01621459.2017.1295863 , in press.
- B. Reich, C. Storlie and H. Bondell, Technometrics, 2009, 51, 110–120 CrossRef PubMed.
- M. C. Kennedy and A. O'Hagan, J. R. Stat. Soc. Ser. B Stat. Methodol., 2001, 63, 425–464 Search PubMed.
- M. Van Oijen, J. Rougier and R. Smith, Tree Physiol., 2005, 25, 915–927 CrossRef PubMed.
- M. Van Oijen, C. Reyer, F. Bohn, D. Cameron, G. Deckmyn, M. Flechsig, S. Härkönen, F. Hartig, A. Huth and A. Kiviste, et al., For. Ecol. Manage., 2013, 289, 255–268 CrossRef.
- D. Higdon, C. Nakhleh, J. Gattiker and B. Williams, Comput. Methods Appl. Mech. Eng., 2008, 197, 2431–2441 CrossRef.
- F. Liu, M. Bayarri, J. Berger, R. Paulo and J. Sacks, Comput. Methods Appl. Mech. Eng., 2008, 197, 2457–2466 CrossRef.
- G. B. Arhonditsis, D. Papantou, W. Zhang, G. Perhar, E. Massos and M. Shi, J. Mar. Syst., 2008, 73, 8–30 CrossRef.
- A. Van Griensven and T. Meixner, J. Hydroinf., 2007, 9, 277–291 CrossRef.
- W. Zhang and G. B. Arhonditsis, J. Great Lakes Res., 2008, 34, 698–720 CrossRef CAS.
- D. Kavetski, G. Kuczera and S. W. Franks, Water Resour. Res., 2006, 42, W03408 Search PubMed.
- E. Jeremiah, S. Sisson, L. Marshall, R. Mehrotra and A. Sharma, Water Resour. Res., 2011, 47, W07547 CrossRef.
- P. M. Tagade, B.-M. Jeong and H.-L. Choi, Build. Environ., 2013, 70, 232–244 CrossRef.
- J.-X. Wang, C. J. Roy and H. Xiao, 2015, arXiv preprint arXiv:1501.03189.
- D. S. Mebane, K. S. Bhat, J. D. Kress, D. J. Fauth, M. L. Gray, A. Lee and D. C. Miller, Phys. Chem. Chem. Phys., 2013, 15, 4355–4366 RSC.
- J. C. Eslick, B. Ng, Q. Gao, C. H. Tong, N. V. Sahinidis and D. C. Miller, Energy Procedia, 2014, 63, 1055–1063 CrossRef CAS.
- B. Konomi, G. Karagiannis, A. Sarkar, X. Sun and G. Lin, Technometrics, 2014, 56, 145–158 CrossRef.
- C. Lai, Z. Xu, W. Pan, X. Sun, C. Storlie, P. Marcy, J.-F. Dietiker, T. Li and J. Spenik, Powder Technol., 2016, 288, 388–406 CrossRef CAS.
- J. Kalyanaraman, Y. Kawajiri, R. P. Lively and M. J. Realff, AIChE J., 2016, 62, 3352–3368 CrossRef CAS.
- J. Kalyanaraman, Y. F. Fan, Y. Labreche, R. P. Lively, Y. Kawajiri and M. J. Realff, Comput. Chem. Eng., 2015, 81, 376–388 CrossRef CAS.
- C. B. Storlie, M. L. Fugate, D. M. Higdon, A. V. Huzurbazar, E. G. Francois and D. C. McHugh, Technometrics, 2013, 55, 436–449 CrossRef.
- A. Berlinet and C. Thomas-Agnan, Reproducing kernel Hilbert spaces in probability and statistics, Kluwer Academic Boston, 2004, vol. 3 Search PubMed.
- A. Berlinet and C. Thomas-Agnan, Reproducing kernel Hilbert spaces in probability and statistics, Springer Science & Business Media, 2011 Search PubMed.
- C. B. Storlie, W. A. Lane, E. M. Ryan, J. R. Gattiker and D. M. Higdon, J. Am. Stat. Assoc., 2015, 110, 68–82 CrossRef CAS.
- W. K. Hastings, Biometrika, 1970, 57, 97–109 CrossRef.
- W. R. Gilks, S. Richardson and D. J. Spiegelhalter, Markov Chain Monte Carlo in practice, Chapman & Hall/CRC, 1996 Search PubMed.
- S. Geman and D. Geman, IEEE Trans. Pattern Anal. Mach. Intell., 1984, 721–741 CrossRef CAS.
- H. Haario, E. Saksman and J. Tamminen, Bernoulli, 2001, 223–242 CrossRef.
- A. Lee and D. C. Miller, Ind. Eng. Chem. Res., 2012, 52, 469–484 Search PubMed.
- R. Bird, W. Stewart and E. Lightfoot, Transport Phenomena (revised second ed.), John Wiley & Sons, 2007 Search PubMed.

## Footnotes |

† The Los Alamos National Laboratory is operated by Los Alamos National Security, LLC for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396. |

‡ The functions defined in this section pertain to thermodynamic quantities and thus need not be strictly dynamic; in the particular calibration case treated here they are still dynamic because T changes with time. |

§ In all figures depicting the conditions inside the bed, time on the x-axis denotes residence time for a population-averaged particle moving through the bed. It is thus proportional to bed depth. |

This journal is © The Royal Society of Chemistry 2017 |