Representing model inadequacy: A stochastic operator approach
Rebecca E. Morrison, Todd A. Oliver, Robert D. Moser
December 8, 2016

arXiv:1604.01651v2 [cs.CE] 7 Dec 2016

Abstract
Mathematical models of physical systems are subject to many uncertainties such as measurement
errors and uncertain initial and boundary conditions. After accounting for these uncertainties, it is
often revealed that discrepancies between the model output and the observations remain; if so, the
model is said to be inadequate. In practice, the inadequate model may be the best that is available
or tractable, and so it may be necessary to use the model for prediction despite its inadequacy. In
this case, a representation of the inadequacy is necessary, so the impact of the observed discrepancy
can be determined. We investigate this problem in the context of chemical kinetics and propose a new
technique to account for model inadequacy that is both probabilistic and physically meaningful. A
stochastic inadequacy operator S is introduced which is embedded in the ODEs describing the evolution
of chemical species concentrations and which respects certain physical constraints such as conservation
laws. The parameters of S are governed by probability distributions, which in turn are characterized
by a set of hyperparameters. The model parameters and hyperparameters are calibrated using highdimensional hierarchical Bayesian inference. We apply the method to a typical problem in chemical
kinetics—the reaction mechanism of hydrogen combustion.

1

Introduction

Model inadequacy is a complex and critical issue that affects nearly all realms of computational science
and engineering. In general, models of physical systems are imperfect: they rely on abstractions and simplifications which do not perfectly represent the modeled system. Sometimes the imperfections are small
enough that any discrepancy between the model and reality is dominated by observation error, such that
the discrepancy is essentially undetectable given existing measurement technology. In contrast, a model is
demonstrably inadequate when the imperfections lead to a detectable inconsistency between the model and
observations. Such inadequacies are often detected during model validation, which is the process of assessing
whether a given mathematical model, including representations of relevant uncertainties, is consistent with
knowledge regarding the modeled system [3, 4, 38]. When an inadequacy is detected, one would generally
prefer to improve the model to remove the discrepancy prior to using it to make predictions, but such improvement is often not feasible. In this case, a representation of the inadequacy is needed to quantify the
possible impact of model inadequacy on computed quantities of interest. However, the formulation of such
model inadequacy representations is an open problem. In this work, we propose a new technique to account
for model inadequacy in the context of chemical kinetics models. Specifically, the model inadequacy in the
equations describing the evolution of the chemical system is represented by a stochastic operator which is
specifically constructed to respect certain non-negotiable constraints on the system. This approach results
in a model inadequacy representation that is both probabilistic and physically meaningful. The method is
demonstrated on a typical problem in chemical kinetics, namely the reaction mechanism model for hydrogen
combustion, where it is shown that the formulation is flexible enough to account for significant inadequacies
present in the original model.
Chemical mechanisms and kinetics models describe the process and rates of chemical reactions [45, 49].
In general, a reaction mechanism is extraordinarily complex, even when there are only two or three initial
reactants. An accurate description of the chemical processes involved in the oxidation of hydrocarbons, for
1

example, may include hundreds or thousands of reactions and fifty or more chemical species [42, 47]. At the
same time, there is significant uncertainty in the reaction rates for these reactions; recent efforts to address
this include [30, 34, 9]. Furthermore, kinetics models of these chemical mechanisms must be embedded
within a larger fluids calculation to model practical combustion systems. The chemical dynamics must then
be represented at every point in space and time. Because the computational cost of such detailed mechanisms
is so high, drastically reduced mechanisms are often used instead. Such reduced models are commonly found
in the combustion literature [11, 47, 50]. However, errors introduced by the reduced models may render the
model inadequate even if the detailed model it is based on is not. Alternatively, it may be that a highly
accurate detailed reaction mechanism is not known, in which case any available model should be viewed
as a reduced mechanism relative to the unknown reality. This work is concerned with accounting for the
inadequacy resulting from the use of a reduced chemical mechanism, though the representation proposed is
equally applicable to inadequate mechanisms.
If a chemical kinetics model is inadequate, it would be best to improve the kinetics model directly to
eliminate the inadequacy. Indeed, refinement of chemical mechanisms in combustion is an active topic of
research; for a small sample focused on H2 /O2 reactions, see [10, 11, 15, 37]. However, this type of refinement
is not always an option, because of a lack of physical insight required to develop a higher fidelity model, a
lack of detailed observational data to support such model development, or a lack of time or other resources
required for the model development process. Further, even when a higher fidelity model exists, it may be
impractical to use due to a lack of the necessary computational resources. Thus, there is a general need for
methods that account for the discrepancy between the inadequate model and the observations that do not
require traditional model improvement. To develop such a representation, one must adopt a mathematical
framework for reasoning about the model inadequacy. Here we adopt a Bayesian point of view [17, 27], and
thus, any lack of knowledge—and model inadequacy specifically—is modeled using probability. Further, the
Bayesian approach offers a natural framework for representing all uncertainties that arise in using reduced
kinetics models to make predictions—including modeling inadequacy of course but also uncertain kinetics
parameters and measurement uncertainties—as well as a natural method for updating these representations
to incorporate information from data [14, 28, 12, 41].
Treatment of model inadequacy has been a topic of research in the Bayesian statistics literature for a
number of years. A common approach is to pose and calibrate a purely statistical model of the discrepancy
between the model output and the true value of that output [29, 4, 25, 26, 46]. However, this representation
is unable to characterize the impact of inadequacy on predictions of an unobserved quantity [38], which is
often the purpose of a model, particularly in science and engineering applications. Further, we generally
have both qualitative and quantitative information about the phenomena being modeled. Incorporating
such information, which is difficult to do in a purely statistical model, is expected to lead to more reliable
predictions.
In this work, inadequacy representations are formulated not as corrections of the model output, but
as stochastic enrichments of the model itself which are specifically constructed to respect known physical
constraints. A set of chemical reactions is modeled by describing the time derivative of the species’ concentrations and temperature, so the model consists of a set of nonlinear ordinary differential equations. The
necessary constraints in this context are conservation of atoms, conservation of energy, and nonnegativity
of concentrations. The model inadequacy representation developed here is characterized by a stochastic
operator S, which is added to the formulation of the time derivatives of the species’ concentrations and
temperature. This operator is constructed such that all realizations of the solution of the enriched reduced
model—i.e. the reduced model coupled with the inadequacy representation—conserve atoms and energy
and have nonnegative concentrations at all times. The main component of the operator is additive, linear,
and probabilistic, encoded in a random matrix S. The use of the term random matrix implies that each
entry is characterized by a probability distribution. This usage is consistent with the definition of random
matrices from random matrix theory (see [19, 33]), although in that field a random matrix is usually much
less constrained than in the present case, and its properties (such as the distributions of the eigenvalues)
are found in the limit as the size of the matrix goes to infinity. A few applications of random matrices to
engineering problems are presented by Soize [43, 44]. However, this work also differs from our approach

2

in that the probability of a given matrix is characterized by properties of the entire matrix, such as the
determinant, whereas we characterize each non-zero entry. Moreover, the inadequacy operator S is more
general than a random matrix alone, including two nonlinear operators as well.
The parameters in the stochastic inadequacy operator S are uncertain because the exact form of the
inadequacy is unknown. Because of this lack of knowledge, the parameters of the inadequacy representation
are characterized by probability distributions, yielding a stochastic forward problem. Each distribution is
characterized by a corresponding set of hyperparameters, which are themselves uncertain. Thus, it is necessary to infer these hyperparameters from data from a detailed model or experiments. We use a hierarchical
Bayesian approach [6] to formulate this inverse problem.
After calibration, the process of validation checks that the model output is consistent with the observations. The particular approach used here is essentially posterior predictive assessment [22, 40], and is
described in the context of physical models in [36]. Specifically, the validity of the model is judged based on
the plausibility of observations according to the uncertain model.
The rest of the paper is organized as follows. In § 2, we give a brief overview of kinetics modeling. In § 3,
the general formulation and properties of the stochastic operator are presented. § 4 describes the Bayesian
framework for calibration and validation of the various models, including hierarchical Bayesian modeling and
validation under uncertainty. In § 5, the approach is applied to the specific case of hydrogen combustion.
Concluding remarks are given in § 6.

2

Chemical mechanism models

Chemical mechanisms and kinetics models describe the process and rates of chemical reactions. In a typical
chemical reaction, there is a set of reactant species which, after a complex series of intermediate reactions,
ultimately form the chemical products. These intermediate steps, in which chemical species react directly
with each other, are called elementary reactions. The set of elementary reactions is called the reaction
mechanism, and a typical combustion problem may include tens to thousands of elementary reactions. This
section provides a minimal introduction to the essentials of chemical kinetics necessary to understand the
development of the model inadequacy representation in § 3. For more details on chemical kinetics, see [45]
for a general text and [49] for a presentation focused on combustion.
To introduce the main concepts, consider the following example reaction set with four species and two
reversible (meaning that the reactions can proceed in either direction—e.g., from A and B to C or from C
to A and B) elementary reactions:
kf

1
−
*
A+B−
)
−
−C

k1b

k2f

−
*
A+C−
)
−
− D,
k2b

where A, B, C, and D denote the different chemical species, and k1f , k2f and k1b , k2b are the forward and
backward rate coefficients, respectively. Let x = [x1 , x2 , x3 , x4 ]T be the vector of molar concentrations
(having dimensions moles per unit volume) corresponding to species A, B, C, D. In the differential equations
describing the evolution of the species’ concentrations, each elementary reaction is a sink for the reactant
species and a sources for the product species. The rate of each reaction is often modeled as linear in the
concentration of the reactants, although this power, or order, associated with a given species may be nonunity. With the assumption that it is linear in each species, the forward rate expressions of the two reactions
are thus
r1f = k1f x1 x2
r2f = k2f x1 x3 .

3

Similarly, the backward rates are given by
r1b = k1b x3
r2b = k2b x4 .
Finally, the ODEs for the molar concentrations are
ẋ1 = −r1f + r1b − r2f + r2b
ẋ2 = −r1f + r1b
ẋ3 = +r1f − r1b − r2f + r2b
ẋ4 = +r2f − r2b .
The rate coefficients k are generally functions of temperature, and may follow a given empirical form
depending on the specific reaction. A common form is the Arrhenius Law,
k(T ) = Ae−E/(R

◦

T)

,

for some prefactor A, activation energy E, and universal gas constant R◦ . Another common form is the
modified Arrhenius,
◦
k(T ) = AT b e−E/(R T ) ,
with the additional constant b. The Arrhenius or modified Arrhenius forms are often used in the literature to
describe the forward rate coefficient, and this is true in the example problem in § 5. However, the backwards
rate coefficients are usually not specified. Instead, these are determined from the equilibrium constant which
depends only on the thermodynamics of the reaction. See [45] for details.
Given a reversible reaction and the forward rate coefficients, there are various software libraries which
will solve for the backwards rate coefficients. In this work, a chemistry software library called Antioch (A
New Templated Implementation of Chemistry for Hydrodynamics) was used to set up the chemical model,
query thermodynamic information, and solve for the reverse reaction rates [1].
To complete the specification of the system, a governing equation for temperature is required. This
equation is derived based on conservation of energy. In this work, we consider a reacting mixture of ideal
gases. Further, the reactions are assumed to occur in a constant volume that does not exchange heat or
mass with its surroundings. In this case, changes in the system temperature are due only to the difference
in chemical energy between the reactants and products. For an ideal gas, the internal energy depends only
on temperature (not v or p) and the species’ concentrations:
X
U (T, x) =
ui (T )xi .
i

Thus,
X ∂ui
X ∂T
∂xi
∂xi
∂T X
∂xi
dU
=
xi + ui
=
cvi
xi + ui
=
cvi x i + u i
.
dt
∂t
∂t
∂t
∂t
∂t i
∂t
i
i
Since the volume is constant—i.e., no work is done on the system—and no heat is added, the change in the
energy U must be zero. Setting dU/dt to zero and solving for dT /dt yields

!
X
dT
1

=− P
ui ẋi  .
dt
i cvi xi
i
Note that dT
dt is a function of both the molar concentrations and their time derivatives. The representations of the remaining functions cvi (T ) and ui (T ) may be found in the literature, generally as seven or nine

4

term polynomials. The commonly used NASA polynomials are used in this work [32].
With the time derivative of temperature, the mathematical model of the reaction mechanism is complete.
To summarize, there is an ODE for the time derivative of each species and also for temperature. These can
be written more compactly as
[ẋ1 , ẋ2 , . . . , ẋn , Ṫ ]T = F(x, ẋ, T )
where F is a nonlinear operator acting on the state (x, T ) and the time derivatives of x. Note that F
depends on ẋ only through the energy equation.

3

Formulation of the model inadequacy

In contrast to the simple example mechanism in § 2, mechanisms describing complex chemical systems
like those encountered in combustion often include hundreds of reactions. One of the standard models for
methane combustion, for example, includes fifty-three species and 325 reactions [42]. Such models are often
referred to as detailed mechanisms and will be written here as D(xD , ẋD , T ). In the context of a reacting
flow simulation, such a large mechanism may be too computationally expensive to be practical, and thus it
is common to use a reduced chemistry model consisting of a subset of the species and reactions from the
detailed model.
To be concrete, suppose that the detailed model includes nD species and mD reactions. Further suppose
that the reduced model includes nR species and mR reactions, where nR ≤ nD and mR < mD . The reduced
model always contains fewer reactions than the detailed; the number of species included in the reduced model
may or may not be smaller, although in practice it is almost always true that nR < nD . The reduced model
is denoted R(xR , ẋR , T ).
In the case that the reduced model does not adequately represent the detailed (or the real chemical
reaction), one has two options: (1) improve the reduced model directly with a more accurate mechanism,
or (2) incorporate a representation of the model error of the reduced model. As noted in § 1, it is often
impractical to improve the model and thus, the focus of this work is on developing a generally applicable
model inadequacy representation for reduced chemistry models. The model inadequacy is represented by a
stochastic operator S that is appended to the reduced model R, as indicated in Figure 1, which shows the
progression from the detailed model to the proposed stochastic model.
In the figure, the D and R operators correspond to the deterministic detailed and reduced models introduced above. The operator RG denotes the reduced mechanism used with a general stochastic inadequacy
model, Gω , where ω is a set of random variables. For example, Gω could be a purely statistical model
that is trained based on observations of the detailed model reaction rates. On the other hand, this term
could represent an augmentation of the chemical mechanism obtained by incorporating more reactions from
the detailed model. In this case, the model inadequacy representation would in fact be deterministic (with
ω = {}). This approach would necessitate more information about the true chemical reaction than we expect
to have or are willing to use and thus is not generally applicable.
Here we take an approach which is intermediate between a strictly statistical representation and a strictly
physics-based approach. Physics is incorporated into the approach in that we insist that the inadequacy
representation respect certain known features of the system, but, since our knowledge of the true dynamics is
incomplete, the model is necessarily stochastic, implying that there is a range of behaviors consistent with our
knowledge. In the chemical kinetics case, we know that errors in the model are due to unrepresented species
and unrepresented pathways by which species transform into each other. Further, both atoms and energy
must be conserved and species’ concentrations and temperature must remain non-negative. To account
for the effects of the missing species and reactions and to satisfy these constraints, a stochastic operator
representation of the inadequacy is posed, as shown on the bottom left of Figure 1. The reduced model RS
plus the stochastic inadequacy operator S is called the operator model O.
For clarity of notation, the superscripts (D, R, G, S) on either the state vector x or the reduced model
R are included to reflect the model at hand. As will become clear, to satisfy the requirements mentioned
above, the inadequacy formulation alters the state vector x and the reduced model operator R, so xS and

5

Detailed model D

Reduced model R
order reduction

D(xD , ẋD , T )

R(xR , ẋR , T )

stochastic
inadequacy
formulation

calibration
data

(RS + S)(xS , ẋS , T )

parameterization,
physical constraints

RG + G ω
General inadequacy
model G

Operator model
O = RS + S

Figure 1: Relationship between the detailed model, deterministic reduced model, reduced model with general
model inadequacy representation, and reduced model with stochastic operator model inadequacy representation.
RS differ from the corresponding xR and R.

3.1
3.1.1

Components and State Variables of the Operator Model
Components of the operator

The main action of the inadequacy operator S is to modify the time derivatives of the species’ concentrations.
The operator consists of three pieces: a random matrix S, a non-linear operator A, and a non-linear operator
W.1 The random matrix S is intended to represent the most general linear correction that can be made
to the reduced chemistry model. However, to ensure sufficient flexibility for mass to move between species,
it is necessary to introduce a non-linear operator A, as discussed in § 3.4. The final piece of the operator
accounts for conservation of energy and is denoted W.
Note that S and A act on just the concentrations, while W acts on the concentrations and their derivatives.
Moreover, it is convenient to formulate S in terms of atomic concentrations, while Ŝ denotes the corresponding
matrix in terms of molar concentrations. This section focuses on S instead of Ŝ because many of the properties
of the matrix (such as non-positivity of eigenvalues) are better expressed in terms of atoms instead of moles
(see § B). To map between S and Ŝ, the vector l is used, whose ith entry counts the number of atoms (of all
types) in one molecule of the ith species. For example, if the set of species is H2 , O2 , H, O, OH, H2 O, then
l = [2, 2, 1, 1, 2, 3]. Let L denote the square matrix with the entries of l on the diagonal. Then Ŝ = L−1 SL
applies to molar concentrations. Finally, putting the three pieces together,
S = Ŝ + A + W
= L−1 SL + A + W,

or more explicitly,

S(xS , ẋS , T ) = L−1 SLxS + A(xS ) + W(xS , ẋS , T ).
1 In general, a script letter refers to a non-linear operator, capital letters to linear operators (matrices), lowercase bold letters
to vectors, and lowercase (unbolded) letters to scalars.

6

3.1.2

Augmentation of the state vector

The reduced model tracks fewer species than the detailed model. It should be possible then, for the inadequacy formulation to represent this difference. However, we do not want to include all the extra (missing)
species in the inadequacy representation. Therefore, in order to account for the missing species in the reduced model, the state space is augmented by entries for all types of atoms. These entries are referred to
as catchall species and act as a sort of pool of each atom type, representing atoms bound in unrepresented
species. The presence of the catchall species allows the operator S to move atoms to and from these pools
instead of constraining every atom to one of the species of the reduced model, which is overly restrictive
because in the detailed model, atoms may move to species that are not part of the reduced model. Thus,
xS is of length nS = nR + nα , where nα is the number of atom types, and is of the form
xS = [x1 , . . . , xnα , xnα +1 , . . . , xnα +nR ]T .
We denote the catchall species of element X by X0 . For example, consider a reduced model that includes
H2 , O2 , OH, and H2 O. Then the catchall species are H0 and O0 , and
xS = [x1 , x2 , x3 , x4 , x5 , x6 ]T
where x1 , . . . , x6 corresponds to H0 , O0 , H2 , O2 , OH, and H2 O, in that order.
The introduction of the catchall species raises an important point about the structure of the reduced
model: it takes on a different form when used in conjunction with the stochastic operator S. Because of
this fact, the reduced model used with the stochastic operator is denoted by RS to distinguish it from the
original reduced model R. There are two differences between these operators. First, RS acts on a vector
space of dimension nS rather than nR , although it has no effect on the first nα entries of xS (i.e., the
catchall species). Second, because the effect of the catchall species on the energy equation is not additive,
the differential equation for T is removed from RS , and the entire calculation is accounted for with W.

3.2

Physical constraints and their implications

There are two non-negotiable constraints that any physically realistic model of the system at hand must respect: (I) conservation of atoms, and (II) non-negativity of concentrations. This ensures that the inadequacy
operator respects physical laws that are known to be true for the systems of interest. This section develops
the implications of these constraints for the components of the operator S with a focus on the matrix S
since, as will be shown, A and W are subsequently chosen to take forms that are known to satisfy these
constraints.
3.2.1

Conservation of atoms

To enforce (I), first let E = [eij ] be the nα × nS matrix where eij is the fraction of atoms of type i in one
molecule of species j. Then ESLxS is the rate of change of the number of each atom due to the operator
S. Because atoms are conserved ESLxS must be zero, which implies that ES = 0.
To continue the example shown in § 3.1, consider the case with atom types H and O, and species H0 , O0 ,
H2 , O2 , OH, H2 O. Then matrix E takes the form:


1 0 1 0 1/2 2/3
E=
.
0 1 0 1 1/2 1/3
To satisfy the constraint that ES = 0, the matrix S is constructed according to
S = CP,
where C is a deterministic matrix and P is probabilistic. The roll of the matrix C is to ensure conservation of
atoms, while P is constructed to ensure that the concentrations are non-negative. To guarantee conservation
7

of atoms, the columns of C must span the nullspace of E, i.e. span(C) = null(E). Thus,
ES = ECP = (EC)P = 0 · P = 0.
E is of dimension nα × nS , so the dimension of the nullspace is nS − nα = nR . Thus C is of dimension
nS × nR and P is of dimension nR × nS .
3.2.2

Non-negativity of concentrations

The second constraint (II) is that the concentrations must not be negative. To see how to enforce this,
consider the differential equation for species Xi 2 :
ẋi = (RS (x, T ))i + (S(x))i = (RS (x, T ))i + (L−1 SLx)i + (A(x))i .

(1)

We must ensure that ẋi ≥ 0 when xi = 0 for i = 1, . . . , nS . The first term of the RHS of (??) is not a
problem, as this is the non-linear part from the reaction mechanism and is thus already physically consistent
[20]. The same argument holds for A(x) because it takes the form of a standard chemical reaction model,
as will be shown in § 3.4. Note that the energy operator W is not written above because it does not modify
the derivative of xi .
Finally, the second term must satisfy the constraint. Although the constraint is naturally posed here in
terms of molar concentrations, it is helpful to rephrase this in terms of atomic concentrations. One can show
that the constraint is satisfied for moles if and only if it is satisfied for atoms:
L−1 SLx ≥ 0

⇐⇒

Sy ≥ 0.

To prove this, consider the ith entry of L−1 SLx:
(L−1 SLx)i = L−1

X

sij lj xj

j

=

1X
sij lj xj
li j

but lj xj = yj and all li > 0, i = 1, . . . , nS . Thus,
X
1X
sij lj xj ≥ 0 ⇐⇒
sij yj ≥ 0.
li j
j
But the final term is exactly the ith element of Sy.
Continuing in terms of the atomic concentrations, the ith component of Sy is given by
X
(Sy)i = sii yi +
sij yj .

(2)

j6=i

The first term from the diagonal, sii yi , automatically respects the constraint: sii may be set to be any
constant
value, since then sii yi → 0 as yi → 0. To enforce the constraint, it must be that the sum,
P
s
y
ij
j , is non-negative. But this sum does not depend on yi , so we choose to set sij ≥ 0 for all i 6= j. This
j6=i
could
be made less restrictive by incorporating information from the non-linear system, i.e. set (RS (x))i +
P
j6=i ŝij xj ≥ 0, but this would violate the linearity assumption on S. It would also necessitate using
information from the reduced model, whereas we aim to constrain the inadequacy operator independently
of R.
2 We

drop the superscript S from x here for ease of notation.

8

3.2.3

Sparsity of S

In practice, many of the entries of S are identically zero. In theory, S could be completely dense if every
species included every type of atom. However, this generally does not occur in practical combustion reactions.
The following proves which entries of S are identically zero, using an argument based on the zeros of the
matrix E.
Theorem 3.1. Consider the ith row of E. Let Ji = {j|eij 6= 0} and Jic = {j|eij = 0}. Then every element
sjk = 0 for j ∈ Ji and k ∈ Jic .
Proof. Consider the ith row of E and the kth column of S. Since ES = 0, we have
X
X
X
X
eij sjk + 0.
eij sjk =
eij sjk +
0=
eij sjk =
j

j∈Jic

j∈Ji

(3)

j∈Ji

But since j and k are in disjoint sets, the sum in line (??) does not include the diagonal term sjj . But the
diagonal term is the only negative value in the k column. Thus, all sjk = 0, where j ∈ Ji and k ∈ Jic .
For another method for determining the sparsity of S, see Appendix B. In addition to sparsity, the
constraints on S imply that it has non-positive eigenvalues. See Appendix B for more details.

3.3

Construction of the matrix S

The structure of S is now clear; the next step is to actually construct it. The challenge in this construction
is that both constraints must be simultaneously satisfied by any realization of S. This subsection presents
a method for construction of the operator. To help demonstrate the upcoming matrix decompositions and
inequality constraints, the construction will also be shown for the example set of species (H0 , O0 , H2 , O2 ,
OH, H2 O). In this case, S has the form


s1,1
0
s1,3
0
s1,5 s1,6
 0
s2,2
0
s2,4 s2,5 s2,6 


s3,1
0
s3,3
0
s3,5 s3,6 
,
S=
 0
s4,2
0
s4,4 s4,5 s4,6 


 0
0
0
0
s5,5 s5,6 
0
0
0
0
s6,5 s6,6
where the diagonal elements are non-positive and the off-diagonal elements are non-negative. Here, nR = 4,
nα = 2.
First, C is formed to span the nullspace of E. To accomplish this, let the bottom nR × nR block of C be
the identity matrix InR . The remaining top nα rows will be the negative of the last nR columns of E. Let
this matrix block be denoted E ∗ . Note that every element of E ∗ is non-positive. Thus, C has the form
#
"
E∗
C=
.
InR
Since E = [Inα | − E ∗ ], it is clear that EC = 0.
For the H2 /O2 example at hand,

1 0
E=
0 1

1
0

1/2
1/2


2/3
.
1/3

−1/2
−1/2


−2/3
−1/3

0
1

Thus,
E∗ =


−1
0

0
−1

9

and


−1
0

1
C=
0

0
0

−1/2
−1/2
0
0
1
0

0
−1
0
1
0
0


−2/3
−1/3

0 
.
0 

0 
1

Next, we construct the nR × nS random matrix P . The first step is to specify which entries are nonnegative, non-positive, or strictly zero. Then, by taking advantage of the special structure of C, it is possible
to transfer the inequalities placed on the entries of S to those of P . Let P1 contain the first nα columns of
P , and P2 the remaining nR columns. So far we have
S = CP
"
#

E∗ 
P1 P2
=
InR
"
#
E ∗ P1 E ∗ P2
=
InR P1 InR P2
#
"
E ∗ P1 |E ∗ P2
.
=
I nR P

(4)

The bottom row of (??) shows how to transfer the inequalities from matrix S to P . Since P is left-multiplied
by the identity matrix, it must be that the signs match for the corresponding elements of S. In particular,
for 1 < i ≤ nR and ∀j, then
pi,j ≤ 0 if s(i+nα ),j ≤ 0

(5)

pi,j ≥ 0 if s(i+nα ),j ≥ 0

(6)

pi,j ≡ 0 if s(i+nα ),j ≡ 0.

(7)

Thus, in the example,

p1,1
 0
P =
 0
0

0
p2,2
0
0

p1,3
0
0
0

0
p2,4
0
0

p1,5
p2,5
p3,5
p4,5


p1,6
p2,6 
,
p3,6 
p4,6

where
p1,3 , p2,4 , p3,5 , p4,6 ≤ 0
and
p1,1 , p1,5 , p1,6 , p2,2 , p2,5 , p2,6 , p3,6 , p4,5 ≥ 0.
Note that the number of non-zero elements in P is 12.
The three inequalities (??-??) are necessary but not sufficient as this only guarantees the inequalities of
the bottom row of (??) hold. The top row introduces more restrictive inequalities on a subset of the entries
of P . First consider the top left block. The only nonzero elements here are the negative entries on the
diagonal. There can be no non-zero off-diagonal elements of S in this block, because each row and column
correspond to a catchall species, and atoms can never move from one catchall to another because they are of
different types, by definition. But all the entries of E ∗ are non-positive, and all entries of P1 are non-negative
by (??) (these correspond to off-diagonal elements of S). Thus, the diagonal elements of S in this top left
block are guaranteed to be non-positive, as required.
Lastly, consider the top right block: E ∗ P2 . To guarantee that these elements are non-negative, it is

10

necessary that the negative entries of P2 (on its diagonal) are large enough in magnitude. For these elements
si,k in the top right block, 1 ≤ i ≤ nα and nα < k ≤ ns . Now
∗
P2(·,k)
0 ≤ si,k = E(i,·)
∗
P(·,k+nα )
= E(i,·)
∗
= E(i,·)
P(·,k0 )
X
=
e∗i,j pj,k0 ,
j

where k 0 = k + nα . The only positive term above in the sum is e∗i,k pk,k0 , so this implies
e∗i,k pk,k0 ≥ −

X

e∗i,j pj,k0 .

j6=k

A similar inequality is placed on the each element pk,k0 for each type of atom (each row of E ∗ that multiplies
the k 0 th column of P ). Therefore, to complete the set of inequalities on P , it is sufficient that, for i = 1, . . . , nα
and k = 1, . . . , nR :
X
1
−pk,k0 ≥
max |e∗i,j |pj,k0 ,
∗
i
mini |ei,k |
j6=k

or, in terms of the matrix C:
− pk,k0 ≥

X
1
max |ci,j |pj,k0 .
i
mini |ci,k |

(8)

j6=k

For use in the following development, denote the RHS of (??) above as qk0 .
In the example, the extra constraints from E ∗ P2 correspond to the diagonal elements of P : p1,3 , p2,4 ,
∗
∗
p3,5 , p4,6 . For example, the constraint s1,5 ≥ 0 implies E(1,·)
P(·,5) ≥ 0 and s2,5 ≥ 0 implies E(2,·)
P(·,5) ≥ 0.
These two constraints are then
1
−1p1,5 − 0p2,5 − p3,5 −
2
1
−0p1,5 − 1p2,5 − p3,5 −
2

2
p4,5 ≥ 0
3
1
p4,5 ≥ 0.
3

Using (??), the two lines above can be condensed into the following inequality which is stronger than either:
2
−p3,5 ≥ 2(p1,5 + p2,5 + p4,5 ).
3
Similarly, the constraints for the other negative elements take the form:
−p1,3 ≥ 0
−p2,4 ≥ 0
1
−p4,6 ≥ 3(p1,6 + p2,6 + p3,6 ).
2

11

ξi
ξ1
ξ2
ξ3
ξ4
ξ5
ξ6
ξ7
ξ8
ξ9
ξ10
ξ11
ξ12

=
=
=
=
=
=
=
=
=
=
=
=
=

pj,k
p1,1
−p1,3
p1,5
p1,6
p2,2
−p2,4
p2,5
p2,6
−p3,5 − 23 (p1,5 + p3,5 + 23 p4,5 )
p3,6
p4,5
−p4,6 − 3(p1,6 + p2,6 + 12 p3,6 )

Table 1: The transformed variables ξ for the example operator.
3.3.1

Transform from P to ξ

Now each element of P is of one of the following forms:
pi,k ≡ 0
pi,k ≥ 0
−pi,k ≥ qk ,

k = i + nα .
n

ξ
These variables can be transformed and reindexed to a new set {ξl }l=1
such that the inequalities take the
simple form ξl ≥ 0 for each l. This mapping also changes from a double-indexed system (pi,j ) to a single
index (ξl ). The index l is introduced because the zero elements of P are not mapped to ξ, so the mapping
is unique to every matrix. For nξ sets {l, i, k}, each ξl is of one of the following two forms:

ξl = pi,k ,

k 6= i + nα

ξl = −(pi,k + qk ),

k = i + nα .

Note that the second set is of size nR and thus the size of the first set is nξ − nR .
For the example, nξ = 12 since there are 12 non-zero elements of P . There are nR = 4 variables whose
transform depends on qk , and thus nξ − nR = 8 variables whose transform does not. The total transform is
given in table 1.
To complete the construction, it remains to specify the probability distribution that governs each variable
ξl . Since ξl ≥ 0, l = 1, . . . , nξ , let
ξl ∼ log N (µξl , ηlξ ).
(9)
The role of the hyperparameters µ and η and how to calibrate them will be explained in detail in §4.For ease
and generality of notation, let ψ be the vector of inadequacy parameters (so far, ψ = ξ but more inadequacy
parameters will be introduced in the upcoming subsections), and let ζ be the vector of all hyperparameters.
This concludes the description of S. Recall that the operator consists of three pieces:
S = Ŝ + A + W.
The next subsections continue with formulations of A and W.

12

(10)

3.4

The catchall reactions A

There is much flexibility in the matrix S with respect to how it can redistribute atoms from certain concentrations to others. In fact, it is the most flexible (or general) linear formulation. That is, at every point in
time, a certain species Xi can be redistributed among all other species Xj as long as sji > 0. Moreover, the
rates at which these processes occur are not set a priori, but are calibrated using the available data. The
random matrix S also provides the flexibility of the catchall species—allowing a place for atoms to go that
might in fact make up a species not included in R but present in D.
However, there is one serious limitation of S due entirely to the linearity: while any species can move to
the catchall species (e.g., H2 O −−→ 2 H0 + O0 ), a catchall species can only directly move to a species made up
of the same type of atom. Therefore, a reaction like the reverse of the previous, namely 2 H0 + O0 −−→ H2 O,
is not allowed. This would require a term that depends on the concentrations of both catchall species, but
in a linear operator this is not possible. In some cases, this limitation is not serious. For example, in the
case with species H2 , O2 , OH, and H2 O, the catchall species could move back to the reduced set of species
since H0 could form H2 and O0 could form O2 .
However, this movement from catchall species back to real species is not always possible. Consider a
methane combustion model that includes the species H2 , O2 , H2 O, CH4 , CO, and CO2 . Then the operator
model species set is H0 , O0 , C0 , H2 , O2 , H2 O, CH4 , CO, and CO2 . Here, S can send carbon atoms from CH4 ,
CO, and CO2 into C0 . But then they are stuck because Cn , for any n = 1, 2, . . . , is not in the reduced set
of species. To overcome the linearity limitation, we introduce a straightforward but nonlinear modification
to the operator: for any species Xi that is made up of more than one type of atom, a non-linear reaction is
included in which the product is Xi and the reactants are the corresponding catchall species. For example,
in the methane case,
κ

1
2 H0 + O0 −→
H2 O

κ

2
4 H0 + C0 −→
CH4

κ

3
O0 + C0 −→
CO

κ

4
2 O0 + C0 −→
CO2 .

This set of reactions is represented by the non-linear operator A. Note that the form is analogous to a
general reaction model. Thus, the constraints (I) and (II) are automatically satisfied.
This modification introduces nκ reaction rate coefficients κ to be calibrated. Similar to the variables ξ,
each κ is positive, by design. Thus,
κ ∼ log N (µκ , η κ ).
(11)
Then ψ is augmented to include these rate coefficients κ, and ζ is augmented to include the additional
hyperparameters µκ and η κ .

3.5

The energy operator W

The third and final component of the operator is the non-linear stochastic energy operator W. The role
of W is to account for temperature changes due to atoms moving into and out of the catchall species. In
other words, allowing for the existence of the catchall species endows them with mass; here the catchall
formulation is completed by endowing them with internal energy.
Recall the differential equation for dT /dt:

! n
S
X
1
dT

= W(x, ẋ, T ) = − PnS
ui (T )ẋi  .
dt
i cvi (T )xi
i
For nα < i ≤ nS , cvi (T ) and ui (T ) are known as functions of temperature from the literature on thermodynamic properties of chemical species [32]. The new contribution is to allow for ui (T ) and cvi (T ) for
13

i = 1, . . . , nα , that is, allow for catchall energies and specific heats and then incorporate these into the calculation of the time derivative of temperature. For actual chemical species, these properties are always given
as a function of temperature. Thus, each new coefficient will also be allowed to have a simple temperaturedependence. Consider a catchall species X0i , i = 1, . . . , nα . For the internal energy, we pose the following
form:
ui (T ) = α0i + α1i T + α2i T 2 ,
and, since cv is its derivative with respect to temperature,
cvi (T ) = α1i + 2α2i T.
Then α0 , α1 , and α2 are additional parameters to be calibrated. Furthermore, like all the other random variables introduced during the modeling of the inadequacy operator, each will in fact be represented
by a probability distribution. This is appropriate since we have incorporated some physical information
(temperature-dependence), but the true functional form is uncertain. It is known that α1 and α2 are positive, while α0 could be positive or negative. These properties are exhibited in probability densities of the
form
α
α0 ∼ N (µα
0 , η0 )

αl ∼

(12)

α
log N (µα
l , ηl ),

l = 1, 2.

(13)

Since the above applies to the nα catchall species, there are 3nα new variables. Of course, ψ and ζ are
again augmented to include the new (and final) inadequacy parameters and hyperparameters. Thus, the
inadequacy parameters are ψ = {ξ, κ, α}, and the hyperparameters are ζ = {µξ , η ξ , µκ , η κ , µα , η α }.
This concludes the description of the stochastic operator S. Many model parameters and hyperparameters
have been introduced for the formulation of the reduced and stochastic operator models; the calibration of
these parameters and validation of the models is discussed in § 4.

3.6

Mapping from the operator to typical reaction form

For given values of the inadequacy parameters ψ, the operator as constructed may be interpreted as providing
a new, enriched chemical reaction model relative to the original reduced model. This fact is important because
it allows relatively straightforward implementation of the stochastic operator model inadequacy approach in
existing software for solving chemical systems. Further, it provides an avenue for physical interpretation of
the operator and associated calibration results.
Since the nonlinear operators A and W are constructed in the usual fashion in this modeling domain,
there is nothing to show. However, the interpretation of the action of the operator Ŝ as a set of chemical
reactions may not be immediately clear. It is now demonstrated that the random matrix Ŝ = L−1 SL can
k P
be mapped to a typical chemical reaction of the form A −
→
βB.
Theorem 3.2. For every j = 1, . . . , nS , the jth column of Ŝ corresponds to the reaction
kj

Xj −→

X

βjp Xp ,

(14)

p6=j

where kj = ŝjj and βjp =

ŝjp

|ŝjj |

.

Proof. Let x be the vector of concentrations of length n (drop the subscript S for ease of notation). Let
the set of reactions above be denoted L(x) (in the same way that the reduced mechanism model is written
R(x)). We will show Ŝx = L(x), element-wise.

14

First,

ŝ1,1 x1 + ŝ1,2 x2 + · · · + ŝ1,n xn


 ŝ2,1 x1 + ŝ2,2 x2 + · · · + ŝ2,n xn 
,

Ŝx = 
..

.


ŝn,1 x1 + ŝn,2 x2 + · · · + ŝn,n xn


(15)

and for a single species Xi ,
(Ŝx)i =

X

ŝij xj .

(16)

j

Now consider L(x). The rate for a particular Xi consists of multiple terms: one in which Xi is the chemical
reactant, and n − 1 terms in which Xi is the chemical product. When Xi is a reactant, the corresponding
rate is −ki xi = ŝii xi . When Xi is a product (and Xj is the reactant, j 6= i), the rates from each reaction are
!
ŝij
+ kj βij xj = ŝjj
xj = ŝij xj , j 6= i.
(17)
ŝjj
Putting the two terms together, we have
(L(x))i = ŝii xi +

X

ŝij xj

(18)

j6=i

=

X

ŝij xj

(19)

j

= (Ŝx)i .

4

(20)

Calibration and validation

This section describes a Bayesian approach to model calibration and validation for the reduced model with
the stochastic inadequacy operator described in § 3. Since our focus is model inadequacy, we take the
reduced model, including rate parameter values, as originally specified, meaning that the rate parameters
in the reduced model are not part of the calibration procedure. Certainly, one could choose to infer the
rate parameters and model inadequacy simultaneously, but this inference is beyond the scope of the current
paper. Instead, we focus on inferring the hyperparameters of the inadequacy operator. This inference can be
accomplished in a straightforward manner using a hierarchical Bayesian approach. This approach, including
details of the data used and prior and likelihood forms, is described in § 4.1. Techniques from posterior
predictive model assessment, which are used to (in)validate both the original reduced model and the reduced
model enriched with the inadequacy operator, are discussed in § 4.2.

4.1

Calibration of the inadequacy operator model

As shown in § 3, the parameters of S, A, and W are characterized by probability distributions whose hyperparameters must be calibrated. Because the primary parameters of interest are actually hyperparameters
characterizing the probability density associated with the parameters that appear directly in the model, it
is natural to pose the calibration problem within the hierarchical Bayesian modeling framework described in
the work of Berliner [6, 48].

15

4.1.1

Hierarchical Bayesian modeling

The calibration problem is formulated as a single Bayesian update for the hyperparameters ζ and the
inadequacy model parameters ψ, given the observations d. Bayes’ theorem requires that
p(ψ, ζ|d) ∝ p(d|ψ, ζ) p(ψ, ζ).
This form can be simplified using the hierarchical structure of the model. First, the map from (ψ, ζ) to the
observables does not depend directly on ζ, the values of ζ are irrelevant after conditioning on ψ. Thus, the
likelihood becomes
p(d|ψ, ζ) = p(d|ψ).
Second, because the distribution for ψ depends on ζ, it is convenient to rewrite the joint prior as
p(ψ, ζ) = p(ψ|ζ) p(ζ).
Thus, the posterior distribution can be written as
p(ψ, ζ|d) ∝ p(d|ψ) p(ψ|ζ) p(ζ).
Clearly, the posterior represents the joint distribution for the hyperparameters ζ and the inadequacy parameters ψ conditioned on the data. However, the particular values of the inadequacy parameters ψ that
are preferred by the given data are not necessarily of interest because the goal is for the formulation to be
applicable to a broad range of problems, including scenarios outside the calibration data set. In this situation, the hyperparameters ζ are the primary target of the inference rather than ψ, and one can marginalize
over ψ to find the joint posterior distribution for hyperparameters:
Z
p(ζ|d) = p(ψ, ζ|d)dψ.
This joint posterior is equivalent to that found by formulating the following inverse problem:
p(ζ|d) ∝ p(d|ζ) p(ζ),
where the likelihood is given by
Z
p(d|ζ) =
4.1.2

p(d|ζ, ψ)p(ψ|ζ)dψ.

Inverse problem details

The data that is to be used in the Bayesian update described above are observations generated by the
detailed chemical kinetics model. Specifically, the calibration data set consists of observations of the molar
concentrations of each of the nR species tracked by the reduced model and temperature, at nt instances
in time, and for nIC initial conditions. The initial condition is given by the set {xf , xo , T }|t=0 and can
be characterized by just two quantities: the equivalence ratio φ and temperature T . The equivalence ratio
quantifies how far the initial condition deviates from the stoichiometric ratio of fuel to oxidizer and is defined
by
xf /xo
,
φ=
xfST O /xoST O
where xfST O and xoST O denote the stoichiometric concentrations of fuel and oxidizer, respectively. Thus, the
initial condition is written as the set IC = {φ, T }.
The true value of the observable (i.e., the output of the detailed model) is denoted by dt and may be

16

indexed as follows:
dtijl = xD
i (tj , ICl ),

i = 1, . . . , nR ;

dtijl = T D (tj , ICl ),

i = nR + 1.

j = 1, . . . , nt ;

l = 1, . . . , nIC ;

For the purposes of calibration, the observations are collected into the vector d, and it is assumed that data
are contaminated by additive Gaussian noise, such that
d = {dijl },
where
dijl = dtijl + ijl .
√
√
and ijl ∼ N (0, σ2ijl ), where σijl = 0.001 ≈ 0.032 if i = 1, . . . , nR and σijl = 1000 ≈ 32 if i = nR + 1.
For simplicity of exposition in the following, this can be reindexed:
di = dti + i ,

i = 1, . . . , nd ,

where nd = (nR + 1)nt nIC .
To formulate the likelihood function, consider the mapping MS from the model inadequacy parameters
ψ to the observables that is induced by the reduced chemistry model after being enriched by the stochastic
operator inadequacy representation. The model claims that
di = MSi (ψ) + i ,

i = 1, . . . , nd .

Thus, the likelihood is given by
p(d|ψ, ζ) = p(d|ψ)
=



1
(2π)nd /2 |Σ|1/2


1
S
T −1
S
exp − (d − M (ψ)) Σ (d − M (ψ)) ,
2

where Σ is the covariance matrix for i .
Following the hierarchical scheme, the joint prior for the inadequacy parameters and hyperparameters is
given by
p(ψ, ζ) = p(ψ|ζ) p(ζ).
The conditional prior distribution p(ψ|ζ) for the inadequacy parameters given the hyperparameters is implied
by the proposed structure in § 3, given in lines (??), (??), (??), (??). That is, the conditional prior
distribution of each inadequacy parameter is the following:
ξi ∼ log N (µξi , ηiξ ),

i = 1, . . . , nξ

κi ∼ log N (µκi , ηiκ ),

i = 1, . . . , nκ

α0i ∼
αli ∼

α
N (µα
i
0i , η0i ),
α
α
log N (µli , ηli ),

= 1, . . . , nα
l = 1, 2, i = 1, . . . , 2nα .

Recall that ξ are the inadequacy parameters of S, κ are those of A, and α of W.
The hyperparameters are taken to be independent in the prior with the following prior distributions:
(·)

(·)

(·)

µi ∼ N (µµi , ηiµ )
(·)

ηi ∼ J (0, ∞),

17

where (·) represents ξ, κ, or α, and J denotes the Jeffreys distribution pJ (x) ∼ J (0, ∞), which is given by
pJ (x) =

1
,
x

x ∈ (0, ∞).

Although the Jeffreys distribution is not normalizable, it is still a valid prior [18].
The posterior distribution p(ψ, ζ|d) can be sampled using Markov chain Monte Carlo sampling methods [16, 23, 24]. The algorithm used in this work is the Delayed Rejection Adaptive Metropolis (DRAM)
algorithm [24]. Specifically, the results shown in § 5 were generated using the DRAM implementation in
the QUESO (Quantification of Uncertainty for Estimation, Simulation, and Optimization) library[13, 39].
QUESO is designed to enable research in Bayesian statistics by providing sampling algorithm implementations that can be used in parallel computing environments. There are other software libraries available to
sample posterior distributions including BUGS (Bayesian Inference Using Gibbs Sampling) [31] and MUQ
(MIT Uncertainty Quantification library) [2].

4.2

Validation

Once a model has been constructed and calibrated, the next step is to validate or assess the consistency
between observations of the modeled system and the calibrated model. The validation approach used here
is that of posterior predictive assessment [22, 40].
v
Consider a set of observations of the system {vi }ni=1
. This set will in general include the data d used
for calibration and may also include additional observations of the same or different quantities not used in
calibration. However, there is an observational error  for each observation so that the observed value vi is
related to the unknown true value vit by
vi = vit + i .
The relevant validation question is whether the observations vi are consistent with the model’s claim regarding
the observation.
Denoting the model parameters as θ, the prediction of the observed value vi according to the calibrated
model is
Z

Z
t
t
p(vi |d) =
p (vi − vi )
p(vi |θ)p(θ|d) dθ dvit .
(21)
dt

θ

Here, depending on the circumstance, the parameters θ may include physical parameters as well as inadequacy parameters ψ and hyperparameters ζ.
Thus, the validation question is whether the available observations are consistent with the model distribution p(vi |d). In the case of the chemical kinetics models considered here, two different validation situations
are relevant. First, in testing the reduced model itself (without inadequacy), θ includes only the kinetic
parameters k. If these parameters are treated as deterministic and known, then p(θ|d) in (??) is the delta
distribution, and the only uncertainty about vi comes from observational error. This approach is used in
§ 5.2 to assess the reduced model alone.
Second, in developing a predictive model, one must assess whether the stochastic inadequacy representation can account for model discrepancies over a broad range of conditions, particularly for conditions not
included in the calibration. In this situation, the question is not whether there exists a ψ that can correct
the model for a given scenario but rather whether the uncertainty in ψ induced by the stochastic form
parameterized by ζ is sufficient to characterize the mismatch between the data and reduced model for the
whole range of scenarios of interest. In this case then, ζ is the target of the calibration problem and the
posterior for ψ is not used. Thus, the distribution for vit is given by
Z Z Z
p(vit |d) =
p(vit |k, ψ, ζ)p(ψ|ζ)p(k, ζ|d)dζdψdk.
k

ψ

ζ

This approach is used in § 5.4.
In each of these cases, the integral (??) yields the posterior prediction of the observation vi which can
18

be used to find the total probability of observing a value less probable than the actual observation. As
explained in [38], this probability can be used as a validation metric, which in turn makes use of highest
probability density (HPD) credibility regions [8]. The β-HPD (0 ≤ β ≤ 1) credibility region S is the set for
which the probability of belonging to S is β and the probability density for any point inside S is higher than
those outside. Define for one observation vi ,
γi = 1 − βmini ,
where βmini is the smallest value of β for which vi ∈ Si . Another way to think of γi is that it is the integral
of p(vit |d) over the domain Vi = {vit : p(vit |d) < p(vi |d)}. For samples {vij }Jj=1 of this distribution p(vit |d),
we have
Z
1X
1vij ∈Vi .
γi =
p(vit |d)dvit ≈
J j
Vi
A delicate point here is the choice of tolerance level τ for which, if γ < τ , the model is deemed inconsistent
with the observation(s). A typical value for the tolerance is 0.05, although there is an extensive discussion in
the statistics literature about how to interpret this [21, 35, 38]. When comparing multiple observations but
treating them as independent, as we will later on, the tolerance should be corrected and set lower because
with many observations of a random variable it is more likely to make a low-probability observation. The
Bonferroni correction suggests dividing the tolerance by the number of points [7]. Ideally, all data points will
be clearly consistent with the model output (the model is not invalidated), or, at least one will be clearly
inconsistent (the model is invalid and thus inadequate).

5

Hydrogen combustion

As an example, the proposed inadequacy operator approach is applied to a chemical mechanism model of
hydrogen combustion. Since there are several stages of the process, it is helpful to summarize the steps:
1. Identify the reduced kinetics model and a data source, which can be a detailed model if one exists
(§5.1).
2. Use a predictive assessment to (in)validate the reduced model (§5.2).
3. If invalid, represent the inadequacy using the stochastic operator method (§5.3).
4. Calibrate the parameters of the stochastic operator using data (§5.4).
5. Use a posterior predictive check to validate the new model (§5.4).
6. Make a prediction (if not invalidated) (§5.5).

5.1

Identification of the reduced model and data source

We investigate a reduced model of H2 /O2 combustion proposed by Williams [50]. For purposes of illustration,
data are generated according to a detailed model also proposed by Williams [50]. In the detailed model,
there are two types of atoms: hydrogen and oxygen; eight distinct species: H2 , O2 , H, O, OH, HO2 , H2 O,
H2 O2 ; and twenty-one elementary reactions. The reduced model contains only seven of these species and five
reactions. The resulting differential equations are much simpler than those given by the full model. Both
the twenty-one and five reaction mechanisms and corresponding forward reaction rate parameters are listed
in appendix A.

19

5.2

Invalidation of the reduced model

Since we choose to view the reduced model as including both the form (i.e., chosen species and reactions)
and parameter values given by the original authors, the reduced model is not calibrated. Instead, the rate
parameters are taken to be deterministic with the values assigned by [50]. Thus to assess the validity
of the reduced model, one must simply compare its predictions of the observations, which are based on
the deterministic model predictions and the stochastic description of the observation error, with data. For
purposes of this comparison, observations are generated using the detailed model for each of the seven species
tracked by the reduced model (nR = 7) plus temperature, at five instances in time (nt = 5), and for one
initial condition (nIC = 1). Thus, nd = 8 × 5 × 1 = 40. The set of time points is {20, 40, 60, 80, 100}µs, and
the initial condition
is φ = 1.0, T0 = 1300 K. The observational error is√modeled as Gaussian with standard
√
deviation σ = 0.001 mol/m3 for the species concentrations and σ = 1000 K for the temperature.
For the reduced model to be declared valid, it is required that the model output make the observations
plausible. Figure 2 shows the observations, generated by the detailed model D, compared to the reduced
model R output. It is clear that there are substantial discrepancies between the reduced model and the
observations for both species concentrations and temperature, shown in figure 3. For example, at the final
observation time, the model under-predicts the observed temperature by approximately 600K, which is far
larger than the observational uncertainty. Further, this level of error is not restricted to a single point.
The temperature is dramatically under-predicted for all observations after ignition, and many of the species
predictions have errors far larger than the observational uncertainty as well. Thus, for any reasonable
tolerance level, inspection of figures 2 and 3 shows that the reduced model alone is invalid and that a model
inadequacy representation is required if it is to be used.

5.3

Formulation of the inadequacy operator S

The elements in the reduced model are H and O. Thus, the catchall species for the inadequacy operator are
H0 and O0 . Thus, the state for the enriched model is
xS = [x1 , x2 , . . . , x9 ]T ,
with the concentrations given in the order of H0 , O0 , H2 , O2 , H, O, OH, HO2 , H2 O. Note that mR = 5,
nR = 7, nα = 2, nS = 9.
5.3.1

The random matrix S

From theorem 3.1, we know that the matrix has the

s1,1
0
s1,3
0
 0
s
0
s
2,2
2,4

s3,1
0
s3,3
0

 0
s
0
s
4,2
4,4

0
s5,3
0
S=
s5,1
 0
s
0
s
6,2
6,4

 0
0
0
0

 0
0
0
0
0
0
0
0

following structure:
s1,5
0
s3,5
0
s5,5
0
0
0
0

0
s2,6
0
s4,6
0
s6,6
0
0
0

s1,7
s2,7
s3,7
s4,7
s5,7
s6,7
s7,7
s8,7
s9,7

s1,8
s2,8
s3,8
s4,8
s5,8
s6,8
s7,8
s8,8
s9,8


s1,9
s2,9 

s3,9 

s4,9 

s5,9 
.
s6,9 

s7,9 

s8,9 
s9,9

Here, of the 81 entries of S, 42 are identically zero. Next E is the nα × nS matrix:


1 0 1 0 1 0 1/2 1/3 2/3
E=
.
0 1 0 1 0 1 1/2 2/3 1/3

20

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

1.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.3

0.8
0.2

O

H

0.6

0.4
0.1
0.2

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.3

2.0

1.5

OH

H2 O

0.2

1.0

0.1
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 2: Concentrations [mol/m3 ] versus time [s] for H2 /O2 reaction, φ = 1.0, T0 = 1300 K. Observations
(red), reduced model R (blue), plotted with 65 and 95% confidence intervals.

21

3000

Temperature

2500

2000

1500

1000

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

Figure 3: Temperature [K] versus time [s] for H2 /O2 reaction, φ = 1.0, T0 = 1300 K. Observations (red),
reduced model R (blue), plotted with 65 and 95% confidence intervals.
The first row of E corresponds to hydrogen
matrix whose columns span null(E):

−1 0
 0 −1

1
0

0
1

0
0
C=

0
0

0
0

0
0
0
0

atoms and the second to oxygen. C is the following nS × nR

P is an nR × nS matrix, where S = CP :

p1,1
0
p1,3
 0
p
0
2,2

p3,1
0
p
3,3

p4,2
0
P =
 0
 0
0
0

 0
0
0
0
0
0

−1
0
0
0
1
0
0
0
0

−1/2 −1/3
−1/2 −2/3
0
0
0
0
0
0
0
0
1
0
0
1
0
0

0
−1
0
0
0
1
0
0
0

0
p2,4
0
p4,4
0
0
0

p1,5
0
p3,5
0
0
0
0

0
p2,6
0
p4,6
0
0
0

p1,7
p2,7
p3,7
p4,7
p5,7
p6,7
p7,7

The transform is given in table 2, and the constraints are now
ξi ≥ 0,

i = 1, . . . , 33.

Finally, to complete the formulation of S,
ξi ∼ log N (µξi , ηiξ ),

22

i = 1, . . . , 33.


−2/3
−1/3

0 

0 

0 
.
0 

0 

0 
1

p1,8
p2,8
p3,8
p4,8
p5,8
p6,8
p7,8


p1,9
p2,9 

p3,9 

p4,9 
.
p5,9 

p6,9 
p7,9

ξi
ξ1
ξ2
ξ3
ξ4
ξ5
ξ6
ξ7
ξ8
ξ9
ξ10
ξ11
ξ12
ξ13
ξ14
ξ15
ξ16
ξ17
ξ18
ξ19
ξ20
ξ21
ξ22
ξ23
ξ24
ξ25
ξ26
ξ27
ξ28
ξ29
ξ30
ξ31
ξ32
ξ33

=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=

pj,k
p1,1
−(p1,3 + p3,3 )
p1,5
p1,7
p1,8
p1,9
p2,2
−(p2,4 + p4,4 )
p2,6
p2,7
p2,8
p2,9
p3,1
p3,3
−(p3,5 + p1,5 )
p3,7
p3,8
p3,9
p4,2
p4,4
−(p4,6 + p2,6 )
p4,7
p4,8
p4,9
−p5,7 − 2(p1,7 + p2,7 + p3,7 + p4,7 + (2/3)p6,7 + (2/3)p7,7 )
p5,8
p5,9
p6,7
−p6,8 − 3(p1,8 + p2,8 + p3,8 + p4,8 + (1/2)p5,8 + (2/3)p7,8 )
p6,9
p7,7
p7,8
−p7,9 − 3(p1,9 + p2,9 + p3,9 + p4,9 + (1/2)p5,9 + (2/3)p6,9 )

Table 2: The transformed variables ξ for the H2 /O2 operator.

23

5.3.2

The catchall reactions A

The catchall reactions allow the catchall species to directly form any species made up of more than one type
of atom. Otherwise, that reaction is already allowed via S (H0 −−→ H is allowed by S for example). Thus,
there are three catchall reactions:
κ

1
H0 + O0 −→
OH

κ

2
HO2
H0 + 2 O0 −→

κ

3
H2 O.
2 H0 + O0 −→

The reaction rate coefficients are denoted κ, and these are included in the set of inadequacy parameters.
Like the variables ξ, each κ will be modeled with a lognormal distribution whose hyperparameters are also
calibrated. From the reactions above, the associated rate of each is:
r10 = κ1 x1 x2
r20 = κ2 x1 x2
r30 = κ3 x1 x2 ,
and the resulting additions to the differential equations for H0 , O0 , OH, HO2 , and H2 O are
H0 :
0

O:

− r10 − r20 − 2r30
−

OH:

+

HO2 :

+

H2 O:

+

r10 −
r10
r20
r30 .

2r20

−

r30

(22)
(23)
(24)
(25)
(26)

The terms above (??) - (??) are written as A(x).
5.3.3

The energy operator W

The third and final piece of the operator S is the energy operator W. Recall

! n
S
X
1
dT

= W(x, ẋ, T ) = − PnS
ui (T )ẋi 
dt
i cvi (T )xi
i
and so a description of u(T ) and cv (T ) for each catchall species is necessary. To do so, the new parameters
α0 , α1 , and α2 are introduced. That is,
ui (T ) = α0i + α1i T + α2i T 2
cvi (T ) = α1i + 2α2i T,
where i = 1 corresponds to H0 and i = 2 to O0 .

5.4

Calibration and validation of the inadequacy operator

For the purposes of calibration and assessment of the stochastic inadequacy operator, we use nine initial
conditions given by the combinations of φ = {0.9, 1.0, 1.1} and initial temperature T0 = {1200, 1300, 1400}K.
The set of time points is again {20, 40, 60, 80, 100}µs, and the observational uncertainty is as described in
§ 5.2. The prior, likelihood, and posterior distributions exactly follow from the general form in § 4.1.

24

Specifically,

(·)

µi



N (0, 1e2) [1/s], if (·) = ξ



N (0, 1e2) [cm3 /mol/s], if (·) = κ

∼ N (0, 1e12) [J/kg/K], if (·) = α0



N (0, 1e2) [J/kg/K2 ], if (·) = α1



N (0, 1e2) [J/kg/K3 ], if (·) = α .
2

Table 3 shows the value of γ for each data point. After calibration, the agreement between the data and
the results of the stochastic operator model is generally very good. In particular, the results are substantially
better than those of the reduced model alone. For example, ten observations out of 360 have γ values equal
to or less than 0.01. One might prefer to have somewhat larger values of γ for these points, but these results
indicate a representation of the model error that should be adequate for most purposes.
Figures 4-12 show the model output of concentrations for the nine different initial conditions. Temperature results are shown in figures 13 and 14. Note that ten observations of the detailed model are shown
in each time series, instead of just the five that were used for calibration. This is simply to better demonstrate the behavior of the detailed model compared to the reduced and stochastic operator models. The
observations used in the calibration correspond to the tick marks on the x-axis.
Consistent with the γ results, the model prediction interval generally includes the observations used for
calibration. The points where the observation is outside the plotted 65 and 95% confidence intervals (and
hence γ is small) occur disproportionately at the first and last point in time in the data set. This suggests
that the operator model may have difficulty accounting for both ignition and equilibrium behavior. Further,
the mixture is initially relatively cool compared to the observations taken later in time. Thus, the difficulty
of the model capturing both the small and large time behavior may indicate that a further enrichment of
the stochastic operator to include a temperature dependence is necessary. However, such an extension is not
mandated by the data and is beyond the scope of the current paper.

5.5

Prediction

We conclude this section by predicting the concentrations and temperature for an extrapolative scenario—
i.e., an initial condition outside the range of those used to generate calibration data. The stochastic operator
approach developed in this work has been specifically designed with this use in mind, and it is formulated
such that there is the potential, with appropriate validation tests, that such an extrapolative prediction
could be supported by available information. While there are still possible improvements to the formulation,
including the temperature dependence discussed in §5.4 and other possible extensions mentioned in §6, the
fact that a validated extrapolative prediction is a possibility distinguishes this formulations from most other
inadequacy representations. Indeed, to the authors’ knowledge, this is the only existing model inadequacy
formulation for chemistry that has the potential to make a validated prediction outside of the range of the
calibration data.
The prediction scenario shown here is a higher equivalence ratio and lower initial temperature than
included in the calibration data set. Specifically, the initial condition is given by φ = 1.15 and T0 = 1150 K.
Although the corresponding output from the detailed model was not used to calibrate the model, the detailed
model output is shown analogously to the previous results. Since the detailed model data are available, this
extrapolation case can be viewed as an additional validation test of the calibrated stochastic operator model.
For completeness of the exercise, the reduced model output for the same prediction scenario is shown
first in figures 15 and 16. As expected based on the calibration results, the reduced model alone is entirely
incapable of capturing the behavior of the system, and many of the points from the detailed model fall
outside the 65 and 95% confidence intervals predicted by the reduced model.
Figures 17 and 18, show results from the stochastic operator model, and the corresponding γ-values are
shown in table 4. Unlike the reduced model alone, for the enriched model with the stochastic operator, all
observations are consistent with the prediction.

25

φ

TIC

0.9

1200

0.9

1300

0.9

1400

1.0

1200

1.0

1300

1.0

1400

1.1

1200

1.1

1300

1.1

1400

t (tj )

H2

O2

H

O

OH

HO2

H2 O

T

1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5

0.86
0.91
0.31
0.62
0.09
0.99
0.61
0.37
0.05
0.01
0.01
0.54
0.07
0.007
0.001
0.98
0.93
0.24
0.74
0.12
0.98
0.74
0.53
0.06
0.009
0.02
0.76
0.09
0.02
0.0008
0.96
0.89
0.91
0.88
0.15
0.95
0.82
0.65
0.07
0.001
0.001
0.92
0.11
0.02
0.0003

0.98
0.91
0.41
0.71
0.65
0.99
0.52
0.71
0.67
0.67
0.04
0.78
0.69
0.71
0.71
0.98
0.99
0.21
0.55
0.62
0.93
0.43
0.59
0.7
0.87
0.02
0.61
0.71
0.85
0.96
0.99
0.93
0.88
0.42
0.58
0.97
0.34
0.47
0.67
0.93
0.01
0.48
0.64
0.92
0.91

0.98
0.95
0.21
0.72
0.26
0.98
0.21
0.48
0.17
0.08
0.01
0.6
0.18
0.05
0.02
0.98
0.89
0.74
0.87
0.5
0.95
0.09
0.87
0.31
0.09
0.01
0.92
0.36
0.07
0.001
0.98
0.96
0.91
0.44
0.83
0.96
0.02
0.71
0.53
0.14
0.004
0.52
0.59
0.13
0.02

0.98
0.99
0.55
0.76
0.57
0.98
0.93
0.69
0.5
0.37
0.66
0.76
0.47
0.33
0.24
0.98
0.99
0.22
0.99
0.82
0.94
0.72
0.95
0.68
0.52
0.49
0.96
0.71
0.49
0.35
0.95
0.98
0.95
0.7
0.99
0.98
0.51
0.81
0.88
0.68
0.31
0.71
0.88
0.63
0.49

0.98
0.98
0.66
0.75
0.48
0.97
0.8
0.7
0.39
0.21
0.67
0.82
0.45
0.22
0.12
0.97
0.95
0.54
0.8
0.57
0.98
0.74
0.76
0.5
0.31
0.54
0.88
0.55
0.32
0.19
0.99
0.93
0.91
0.82
0.66
0.96
0.73
0.84
0.6
0.4
0.31
0.93
0.65
0.4
0.28

0.95
0.99
0.96
0.98
0.96
0.99
0.96
0.93
0.94
0.98
0.95
0.95
0.98
0.89
0.97
0.98
0.97
0.93
0.99
0.96
0.97
0.96
0.98
0.97
0.95
0.98
0.99
0.96
0.89
0.89
0.93
0.96
0.96
0.99
0.91
0.95
0.95
0.93
0.97
0.96
0.98
0.96
0.92
0.93
0.98

0.95
0.94
0.34
0.81
0.73
0.99
0.17
0.98
0.69
0.59
0.02
0.95
0.68
0.54
0.49
0.98
0.99
0.23
0.53
0.99
0.98
0.08
0.69
0.85
0.6
0.009
0.58
0.86
0.56
0.47
0.95
0.97
0.97
0.28
0.79
0.99
0.04
0.44
0.96
0.7
0.001
0.34
0.91
0.67
0.5

0.96
0.96
0.49
0.76
0.92
0.94
0.49
0.94
0.95
0.98
0.23
0.99
0.9
0.99
0.91
0.98
0.96
0.53
0.64
0.77
0.99
0.43
0.85
0.87
0.91
0.07
0.88
0.89
0.98
0.91
0.98
0.93
0.83
0.49
0.55
0.99
0.4
0.59
0.74
0.76
0.01
0.76
0.81
0.92
0.88

Table 3: The values of γ for each data point using the H2 /O2 operator.

26

2.0

1.5

1.5

O2

H2

1.0

1.0

0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.8

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.25

0.20
0.6

0.4

O

H

0.15

0.10

0.2
0.05

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 4: Concentrations [mol/m3 ] versus time [s], φ = 0.9, T0 = 1200 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

27

2.0

1.5

1.5

O2

H2

1.0

1.0

0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.8

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.3

0.6

0.4

O

H

0.2

0.1
0.2

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 5: Concentrations [mol/m3 ] versus time [s] for φ = 0.9, T0 = 1300 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

28

2.0

1.5

1.5

O2

H2

1.0

1.0

0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.8

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.3

0.6

0.4

O

H

0.2

0.1
0.2

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 6: Concentrations [mol/m3 ] versus time [s] for φ = 0.9, T0 = 1400 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

29

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.8

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.25

0.20
0.6

0.4

O

H

0.15

0.10

0.2
0.05

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 7: Concentrations [mol/m3 ] versus time [s] for φ = 1.0, T0 = 1200 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

30

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.8

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.3

0.6

0.4

O

H

0.2

0.1
0.2

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 8: Concentrations [mol/m3 ] versus time [s] for φ = 1.0, T0 = 1300 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

31

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.8

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.3

0.6

0.4

O

H

0.2

0.1
0.2

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 9: Concentrations [mol/m3 ] versus time [s] for φ = 1.0, T0 = 1400 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

32

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.25

0.8

0.20

0.6

0.15

0.4

0.10

0.2

0.05

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

O

H

1.0

0.00002

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 10: Concentrations [mol/m3 ] versus time [s] for φ = 1.1, T0 = 1200 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

33

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.25

0.8

0.20

0.6

0.15

0.4

0.10

0.2

0.05

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

O

H

1.0

0.00002

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 11: Concentrations [mol/m3 ] versus time [s] for φ = 1.1, T0 = 1300 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

34

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

1.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.3

0.8
0.2

O

H

0.6

0.4
0.1
0.2

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.0

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 12: Concentrations [mol/m3 ] versus time [s] for φ = 1.1, T0 = 1400 K. Observations (red), operator
model O (blue), plotted with 65 and 95% confidence intervals.

35

3000

3000
Temperature

4000

Temperature

4000

2000

1000

2000

0.00002

0.00004

0.00006
Time [s]

0.00008

1000

0.00010

5000

5000

4000

4000

3000

2000

1000

0.00008

0.00010

2000

0.00002

0.00004

0.00006
Time [s]

0.00008

1000

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

(d) φ = 1.0, T0 = 1300 K

5000

5000

4000

4000

Temperature

Temperature

0.00006
Time [s]

3000

(c) φ = 0.9, T0 = 1300 K

3000

2000

1000

0.00004

(b) φ = 1.0, T0 = 1200 K

Temperature

Temperature

(a) φ = 0.9, T0 = 1200 K

0.00002

3000

2000

0.00002

0.00004

0.00006
Time [s]

0.00008

1000

0.00010

(e) φ = 0.9, T0 = 1400 K

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

(f) φ = 1.0, T0 = 1400 K

Figure 13: Temperature [K] versus time [s] for H2 /O2 reaction, φ = {0.9, 1.0}, T0 = {1200, 1300, 1400}K.
Observations (red), operator model O (blue), plotted with 65 and 95% confidence intervals.

36

4000

Temperature

3000

2000

1000

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

(a) φ = 1.1, T0 = 1200 K
5000

Temperature

4000

3000

2000

1000

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

(b) φ = 1.1, T0 = 1300 K
5000

Temperature

4000

3000

2000

1000

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

(c) φ = 1.1, T0 = 1400 K

Figure 14: Temperature [K] versus time [s] for H2 /O2 reaction, φ = 1.1, T0 = {1200, 1300, 1400}K. Observations (red), operator model O (blue), plotted with 65 and 95% confidence intervals.

37

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.25

0.8

0.20

0.6

0.15

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

O

H

1.0

0.00002

0.4

0.10

0.2

0.05

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00

0.00010

0.25

2.0

0.20
1.5

OH

H2 O

0.15
1.0

0.10

0.5
0.05

0.00

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 15: Predicted concentrations [mol/m3 ] versus time [s] for φ = 1.15, T0 = 1150 K. Observations (red),
reduced model R (blue), plotted with 65 and 95% confidence intervals.

38

3000

Temperature

2500

2000

1500

1000

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

Figure 16: Predicted temperature [K] versus time [s] for φ = 1.15, T0 = 1150 K. Observations (red), reduced
model R (blue), plotted with 65 and 95% confidence intervals.

φ

1.15

TIC

j (tj )

H2

O2

H

O

OH

HO2

H2 O

T

1150

1
2
3
4
5

0.93
0.98
0.71
0.13
0.99

0.98
0.97
0.87
0.14
0.36

0.98
0.94
0.88
0.75
0.27

0.99
0.95
0.98
0.43
0.59

0.98
0.93
0.98
0.44
0.79

0.96
0.99
0.98
0.94
0.99

0.93
0.96
0.77
0.14
0.12

0.95
0.95
0.85
0.58
0.47

Table 4: The values of γ for the prediction case.

39

2.5

1.5

2.0
1.0

O2

H2

1.5

1.0
0.5
0.5

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

0.25

0.8

0.20

0.6

0.15

0.4

0.10

0.2

0.05

0.0

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00

0.00010

0.4

2.0

0.3

1.5

H2 O

OH

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

O

H

1.0

0.00002

0.2

0.1

0.0

1.0

0.5

0.00002

0.00004

0.00006
Time [s]

0.00008

0.0

0.00010

Figure 17: Predicted concentrations [mol/m3 ] versus time [s] for φ = 1.15, T0 = 1150 K. Observations (red),
operator model O (blue), plotted with 65 and 95% confidence intervals.

40

4000

Temperature

3000

2000

1000

0.00002

0.00004

0.00006
Time [s]

0.00008

0.00010

Figure 18: Predicted temperature [K] versus time [s] for φ = 1.15, T0 = 1150 K. Observations (red), operator
model O (blue).

6

Conclusion

This study addresses the critical problem of model inadequacy that affects nearly all mathematical models of
physical systems. In particular, we develop a novel approach to representing model inadequacy in chemical
kinetics models. The approach relies on a stochastic operator that combines the flexibility and generality
of a probabilistic model with the available deterministic physical information. In the context of predictive
models, these two properties are essential, ensuring that the model is flexible enough to adapt to fit available
observations without generating non-physical behavior.
The stochastic operator S contains three main components: 1) the random matrix S, 2) the nonlinear
catchall reactions A, and 3) the energy operator W. The random matrix S contains most of the information
in S and has some interesting properties. Typically, the matrix has many identically zero entries. It always
has a negative diagonal, is diagonally dominant, and has non-positive eigenvalues. The reactions in A allow
any species in the reduced model to be the chemical product of the corresponding catchall species (if this
is not already possible through S). Both S and A guarantee conservation of atoms and non-negativity of
concentrations. Finally, the energy operator W modifies the time derivative of temperature by endowing the
catchall species with internal energy.
The inadequacy operator is tested on an example H2 /O2 mechanism. Starting with a reduced mechanism
that was shown to be invalid for prediction, the stochastic operator model is able to account for the model
inadequacy. The observations used for calibration are plausible outcomes of the operator model—i.e., the
concentrations and temperature values from the operator model O are consistent with the data supplied by
the detailed model D. Moreover, the prediction case showed that the operator model output was consistent
with most of the observations from the detailed model for a scenario outside the range of initial conditions
used for the operator calibration.
There are many avenues for extending the work reported here. Ongoing research includes the application of the operator to a more complex chemical setting, namely, a methane-air mechanism. We are also
investigating a more complex physical setting with a more complex prediction problem: the prediction of a
hydrogen laminar flame using the calibrated operator found in this work. Consistent with the discussion in
§ 5.4, preliminary research on these topics has suggested that the temperature-independence of the operator
(save W) is too severe a limitation, and that in complex problems, the more information incorporated into
the stochastic operator, the better. Thus, immediate goals are to include more physical information in the
inadequacy representation such as realistic temperature-dependence, and to use stronger priors based on
knowledge of the chemical reactions and the physical setup. This leads to the next major opportunity for
future work: developing the connection between the stochastic operator and the actual chemistry. Mapping
41

between the random matrix and the typical chemical reactions was a first step in this direction. However, a
better understanding of what the stochastic operator means in physical terms is needed. This includes not
just the structure, but also the uncertainty in the calibrated parameters. A future goal is to infer something
about the missing chemistry from the calibrated operator. This is important when developing mechanisms
based on experimental data, when no detailed mechanism is available.
Further, it is unclear that the random matrix coupled with catchall species and reactions is the optimal
way to formulate the stochastic operator. A number of alternative formulations are being pursued. For
example, instead of using the random matrix S and the catchall reactions in A, a simplified version is the
following: given nR species in the reduced model, include nR reversible dissociation reactions where the
reactant is one of the original species and the products are the corresponding catchall species. In the reverse
combination reaction, the catchall atoms react to form any of the original species. In this fashion, the atoms
of any species could move to any other species in two steps. This representation would lose the flexibility
possible in the current formulation; that is, the current formulation allows for the most detailed and direct
linear movement from one species to another. If any such direct pathway proves significant, this would not
be captured by this simplified formulation. On the other hand, it would decrease the number of random
variables and thus would be more tractable in more complex reaction systems. With a smaller number of
additional reactions, the corresponding reaction rates could then be enriched with temperature-dependence.
Another variation could be a more complete set of nonlinear reactions. Instead of only allowing the
nonlinear catchall reactions, one could augment the reduced model with all or some subset of all possible
nonlinear terms. In contrast to the first variation, this would increase the number of random variables. A
formulation like this might only be possible with more informative priors or knowledge about the chemical
system.
Finally, it would be very informative to apply this method to new problems. For instance, the inadequacy
operator could be tested by a more realistic combustion problem. It may be that doing so requires a more
complete thermodynamic description of the catchall species. More importantly, there is no reason to restrict
the stochastic operator approach developed here to models of chemically reacting gas mixtures. Similar
model structures and inadequacies appear in many other modeling domains. The guiding principles of
this work (respecting physical constraints, maintaining flexibility, starting with a linearized version) and
developing an analogous operator (possibly random matrix) should be explored in many different physical
settings. Applications in many different domains could bring to light new challenges and common strengths
for the stochastic operator approach to representing model inadequacy.

42

Reaction
Hydrogen-oxygen chain
1. H + O2 −−→ OH + O
2. H2 + O −−→ OH + H
3. H2 + OH −−→ H2 O + H
4. H2 O + O −−→ OH + OH

A

b

E

3.52 × 1016
5.06 × 104
1.17 × 109
7.60 × 100

-0.7
2.7
1.3
3.8

71.4
26.3
15.2
53.4

Direct recombination
5. H + H + M −−→ H2 + M
6. H + OH + M −−→ H2 O + M
7. O + O + M −−→ O2 + M
8. H + O + M −−→ OH + M
9. O + OH + M −−→ HO2 + M

1.30 × 1018
4.00 × 1022
6.17 × 1015
4.71 × 1018
8.00 × 1015

-1.0
-2.0
-0.5
-1.0
0.0

0.0
0.0
0.0
0.0
0.0

Hydroperoxyl reactions
10. H + O2 + M −−→ HO2 + M
11. HO2 + H −−→ OH + OH
12. HO2 + H −−→ H2 + O2
13. HO2 + H −−→ H2 O + O
14. HO2 + O −−→ OH + O2
15. HO2 + OH −−→ H2 O + O2

5.75 × 1019
7.08 × 1013
1.66 × 1013
3.10 × 1013
2.00 × 1013
2.89 × 1013

-1.4
0.0
0.0
0.0
0.0
0.0

0.0
1.2
3.4
7.2
0.0
-2.1

Hydrogen peroxide reactions
16. OH + OH + M −−→ H2 O2 + M
17. HO2 + HO2 −−→ H2 O2 + O2
18. H2 O2 + H −−→ HO2 + H2
19. H2 O2 + H −−→ H2 O + OH
20. H2 O2 + OH −−→ H2 O + HO2
21. H2 O2 + O −−→ HO2 + OH
Units: mol, cm,

2.30 × 1018
3.02 × 1012
4.79 × 1013
1.00 × 1013
7.08 × 1012
9.63 × 106
s, kJ, K.

-0.9
0.0
0.0
0.0
0.0
2.0

-7.1
5.8
33.3
15.0
6.0
2.0

Table A.1: The detailed H2 /O2 reaction mechanism from [50].

Appendices
A

Reaction mechanisms

The 21 reactions in the detailed hydrogen-oxygen mechanism are listed in table A.1 and the five of the
◦
reduced mechanism in A.2. The associated reaction rate is k = AT b e−E/R T .

B
B.1

Properties of S
Sparsity

As noted in § 3.2.3, S is often sparse. There is a way to determine which entries of S are identically zero using
the physical restrictions about how different species’ concentrations interact with each other. To determine
the sparsity in this way, each species X is characterized by a composite number ρX . First associate a prime
number pi with each atom type i = 1, . . . nα . Each species X is made up of a collection of atom types; let
ρX be the product of prime numbers corresponding to each type of atom making up species X. For example,
if elements H and O correspond to the prime numbers 2 and 3, then ρH = 2 and ρH2 O = 6. In effect, this
yields a prime number representation of each species where multiplicity is ignored.

43

Reaction
Hydrogen-oxygen chain
1. H + O2 −−→ OH + O
2. H2 + O −−→ OH + H
3. H2 + OH −−→ H2 O + H

A

b

E

3.52 × 1016
5.06 × 104
1.17 × 109

-0.7
2.7
1.3

71.4
26.3
15.2

Hydroperoxyl reactions
10. H + O2 + M −−→ HO2 + M
5.75 × 1019
12b. H2 + O2 −−→ HO2 + H
1.4 × 1014
Units: mol, cm, s, kJ, K.

-1.4
0.0

0.0
249.5

Table A.2: The reduced H2 /O2 reaction mechanism from [50].
Next, the columns of S correspond to chemical reactants and the rows to chemical products. The entry
sij controls how many atoms move from species j- a sort of reactant- to species i- a sort of product. The
operator can only move a positive amount of species Xj to Xi if the former contains all elements that comprise
the latter. If not, then the gcd(ρXi , ρXj ) < ρXi . But in this case there can be no flow of atoms from Xj to
Xi , and thus sij ≡ 0.
This technique can also be used to count the total number of entries that are identically zero in the
matrix. Call this total number Ω and, for i = 1, . . . , nS , let λi be the number of species Xj such that
gcd(ρXi , ρXj ) < ρXi . By the argument in the previous
paragraph, λi is the number of zeros in the ith column
PnS
of S. Then the number of zeros in S is Ω = i=1
λi because the sum is taken with respect to the different
species, and each of these correspond to a different column of S.

B.2

Non-positivity of eigenvalues

Enforcing the two constraints— (I) conservation of atoms and (II) non-negativity of concentrations yields
some interesting properties of the random matrix S. The non-positivity of eigenvalues is consistent with the
constraints: no species can grow arbitrarily large over time. The proof follows:
Theorem B.1. Let S be any random matrix such that ES = 0 and the off-diagonal elements of S be nonnegative. Then (a) the columns sum to zero, (b) the diagonal is negative, (c) the matrix is weakly diagonally
dominant, and (d) the eigenvalues are non-positive.
Proof. (a) Consider ES(·,j) = 0, where S(·,j) is the jth column of S. There are nα equations:
e1,1 s1,j + e1,2 s2,j + · · · + e1,nS snS ,j = 0
e2,1 s1,j + e2,2 s2,j + · · · + e2,nS snS ,j = 0
..
.
enα ,1 s1,j + enα ,2 s2,j + · · · + enα ,nP snS ,j = 0.
Now add the lines together:
XX
i

but

P

k ek,i

ek,i si,j = 0,

(27)

si,j = 0.

(28)

k

= 1 by definition. Thus,
X
i

44

(b) In equation (??), move the diagonal term to the RHS:
X
si,j = −sj,j .

(29)

i6=j

Since all off-diagonal terms are non-negative, it must be that the diagonal element is negative.
(c) The line above also shows weak diagonal dominance, since
|sj,j | =

X

si,j

(30)

|si,j |,

(31)

i6=j

=

X
i6=j

where the second equality holds because all off-diagonal elements are non-negative.
T
(d) Since
have the same eigenvalues, we will show that the claim is true for S T . Let B = S T
P S and S P
and Bi = j6=i |bi,j | = j6=i bi,j be the sum of off-diagonals in the ith row. Now let D(bi,i , Bi ) be the closed
disc centered at bi,i with radius Bi . Then the Gershgorin theorem states that every eigenvalue of B lies
within at least one of the discs [5]. In this case, we have bi,i = si,i and Bi = |si,i |, so every eigenvalue lies
within at least one disc D(si,i , |si,i |), where si,i ≤ 0.

References
[1] A
New
Templated
Implementation
https://github.com/libantioch/antioch.

of

Chemistry

Hydrodynamics

(Antioch).

[2] MIT Uncertainty Quantification Library (MUQ). https://bitbucket.org/mituq/muq.
[3] I. Babuska, F. Nobile, and R. Tempone. A systematic approach to model validation based on Bayesian
updates and prediction related rejection criteria. Computer Methods in Applied Mechanics and Engineering, 197:2517–2539, May 2008.
[4] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C.-H. Lin, and J. Tu. A
framework for validation of computer models. Technometrics, 49(2):138–154, May 2007.
[5] H. E. Bell. Gershgorin’s theorem and the zeros of polynomials. American Mathematical Monthly,
72:292–295, 1965.
[6] M. L. Berliner. Hierarchical Bayesian time series models. In K. M. Hanson and R. N. Silver, editors,
Maximum Entropy and Bayesian Methods, volume 79 of Fundamental Theories of Physics, pages 15–22.
Springer Netherlands, 1996.
[7] C. E. Bonferroni. Il calcolo delle assicurazioni su gruppi di teste. In Studi in Onore del Professore
Salvatore Ortu Carboni, pages 13–60. Rome, 1935.
[8] G. E. P. Box and G. C. Tiao. Bayesian Inference in Statistical Analysis. Wiley-Interscience, 1 edition,
Apr. 1992.
[9] K. Braman, T. A. Oliver, and V. Raman. Bayesian analysis of syngas chemistry models. Combustion
Theory and Modelling, 17(5):858–887, 2013.

45

[10] M. P. Burke, M. Chaos, Y. Ju, F. L. Dryer, and S. J. Klippenstein. Comprehensive H2/O2 kinetic
model for high-pressure combustion. International Journal of Chemical Kinetics, 44(7):444–474, Dec.
2011.
[11] M. P. Burke, M. Chaos, Y. Ju, F. L. Dryer, and S. J. Klippenstein. Comprehensive H2/O2 kinetic
model for high-pressure combustion. International Journal of Chemical Kinetics, 44(7):444–474, July
2012.
[12] D. Calvetti and E. Somersalo. An Introduction to Bayesian Scientific Computing: Ten Lectures on
Subjective Computing, 2007.
[13] Center for Predictive Engineering and Computational Sciences. Quantification of Uncertainty for Estimation, Simulation, and Optimization (QUESO) User’s Manual, 2013.
[14] R. P. Christian. The Bayesian Choice. Springer, 2001.
[15] Conaire, H. J. Curran, J. M. Simmie, W. J. Pitz, and C. K. Westbrook. A comprehensive modeling
study of hydrogen oxidation. International Journal of Chemical Kinetics, 36(11):603–622, Nov. 2004.
[16] M. K. Cowles and B. P. Carlin. Markov chain Monte Carlo convergence diagnostics: A comparative
review. American Statistical Association, 91(434):883–904, June 1996.
[17] R. T. Cox. The Algebra of Probable Inference. Johns Hopkins University Press, 1961.
[18] M. H. DeGroot and M. J. Schervish. Probability and Statistics (4th Edition). Pearson, 4 edition, Jan.
2011.
[19] A. Edelman and N. R. Rao. Random matrix theory. Acta Numerica, 14(-1):233–297, May 2005.
[20] M. Feinberg. Lectures on Chemical Reaction Networks. 1979.
[21] A. Gelman. Exploratory Data Analysis for Complex Models. Journal of Computational and Graphical
Statistics, 13(4):755–779, Dec. 2004.
[22] A. Gelman, X.-l. Meng, and H. Stern. Posterior predictive assessment of model fitness via realized
discrepancies. In Statistica Sinica, pages 733–807, 1996.
[23] W. R. Gilks, S. Richardson, and D. Spiegelhalter. Markov Chain Monte Carlo in Practice. Chapman
and Hall/CRC, softcover reprint of the original 1st ed. 1996 edition, Dec. 1996.
[24] H. Haario, M. Laine, A. Mira, and E. Saksman. DRAM: Efficient adaptive MCMC. Statistics and
Computing, 16(4):339–354, Dec. 2006.
[25] D. Higdon, J. Gattiker, B. Williams, and M. Rightley. Computer model calibration using highdimensional output. Journal of the American Statistical Association, 103(482):570–583, June 2008.
[26] D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne. Combining field data and
computer simulations for calibration and prediction. SIAM Journal on Scientific Computing, 26(2):448–
466, Jan. 2004.
[27] E. T. Jaynes. Probability Theory: The Logic of Science. Cambridge University Press, 2003.
[28] J. Kaipio and E. Somersalo. Statistical and Computational Inverse Problems. Springer, 2005.
[29] M. C. Kennedy and A. O’Hagan. Bayesian calibration of computer models. Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 63(3):425–464, Jan. 2001.
[30] A. A. Konnov. Remaining uncertainties in the kinetic mechanism of hydrogen combustion. Combustion
and Flame, 152(4):507–528, Mar. 2008.
46

[31] D. Lunn, A. Thomas, N. Best, and D. Spiegelhalter. WinBUGS - A Bayesian modelling framework:
Concepts, structure, and extensibility. 10(4):325–337, 2000.
[32] B. J. McBride, S. Gordon, and M. A. Reno. NASA Technical Memorandum 4513: Coefficients for
Calculating Thermodynamic and Transport Properties of Individual Species. Technical report, National
Aeronautics and Space Administration, 1993.
[33] M. L. Mehta. Random Matrices, Volume 142, Third Edition (Pure and Applied Mathematics). Academic
Press, 3 edition, Nov. 2004.
[34] K. Miki, S. H. Cheung, E. E. Prudencio, and P. L. Varghese. Bayesian uncertainty quantification of
recent shock tube determinations of the rate coefficient of reaction H + O2 OH + O. International
Journal of Chemical Kinetics, 44(9):586–597, July 2012.
[35] R. G. Miller. Simultaneous Statistical Inference (Springer Series in Statistics). Springer, 2nd edition,
Mar. 1981.
[36] R. D. Moser and T. A. Oliver. Validation of physical models in the precense of uncertainty. In
R. Ghanem, D. Higdon, and H. Owhadi, editors, Handbook of Uncertainty Quantification. Springer,
2016.
[37] M. A. Mueller, T. J. Kim, R. A. Yetter, and F. L. Dryer. Flow reactor studies and kinetic modeling of
the H2/O2 reaction. Int. J. Chem. Kinet., 31(2):113–125, Jan. 1999.
[38] T. A. Oliver, G. Terejanu, C. S. Simmons, and R. D. Moser. Validating predictions of unobserved
quantities. Computational Methods in Applied Mechanics and Engineering, 283:1310–1335, 2015.
[39] E. Prudencio and S. H. Cheung. Parallel adaptive multilevel sampling algorithms for the Bayesian
analysis of mathematical models. International Journal for Uncertainty Quantification, 2(3):215–237,
2012.
[40] D. B. Rubin. More powerful randomization-based p-values in double-blind trials with non-compliance.
Statist. Med., 17(3):371–385, Feb. 1998.
[41] D. S. Sivia and J. Skilling. Data Analysis: A Bayesian Tutorial. Oxford University Press, 2006.
[42] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eiteneer, M. Goldenberg, C. T. Bowman,
R. K. Hanson, S. Song, W. C. Gardiner, V. V. Lissianski, and Q. Zhiwei. GRI-Mech. http://www.me.
berkeley.edu/gri_mech/. Accessed: 2015-01-15.
[43] C. Soize. Random matrix theory for modeling uncertainties in computational mechanics. Computer
Methods in Applied Mechanics and Engineering, 194(12-16):1333–1366, Apr. 2005.
[44] C. Soize. Generalized probabilistic approach of uncertainties in computational dynamics using random
matrices and polynomial chaos decompositions. Int. J. Numer. Meth. Engng., 81(8):939–970, Feb. 2010.
[45] J. I. Steinfeld, J. S. Francisco, and W. L. Hase. Chemical Kinetics and Dynamics. Prentice Hall, 1998.
[46] L. P. Swiler. Gaussian processes in response surface modeling. Society of Experimental Mechanics
Meeting, Jan. 2006.
[47] P. R. Westmoreland. Reduced kinetic mechanisms for applications in combustion systems. Edited by
N. Peters and B. Rogg, Spring-Verlag, New York, Lecture Notes in Physics, Monograph 15, 1993, 360
pp. AIChE J., 40(11):1926–1927, Nov. 1994.
[48] C. K. Wikle, R. F. Milliff, D. Nychka, and L. M. Berliner. Spatiotemporal hierarchical Bayesian modeling: Tropical ocean surface winds. Journal of the American Statistical Association, 96(454):382–397,
June 2001.
47

[49] F. A. Williams. Combustion Theory: Second Edition (Combustion Science and Engineering). Westview
Press, second edition edition, Mar. 1994.
[50] F. A. Williams. Detailed and reduced chemistry for hydrogen autoignition. Journal of Loss Prevention
in the Process Industries, 21(2):131–135, Mar. 2008.

48

