1 Introduction↩︎

However, in the process of model development, assumptions must be made that can degrade the trustworthiness of these models. Such assumptions include, but are not limited to, assumptions regarding invariance under transformation (e.g. rotation: isotropy, translation: homogeneity), couplings between phenomena, and dependencies (e.g., ignoring microscale dependence, assuming rate dependence).

MFU and other sources of uncertainty, such as parameter uncertainty, result in uncertainties in model predictions that must be quantified to adequately characterize the trustworthiness of resulting predictions. Being able to identify the most important sources of uncertainty is critical for determining where to invest resources to most improve prediction reliability through model improvement, uncertainty reduction, etc. To this end, this work presents a first-of-its-kind method to quantify the impact of modeling assumptions on model predictions relative to other sources of uncertainty.

The first established and most common approach was presented in [1], wherein a statistical model is appended to the output of a physical model, and its hyperparameters are calibrated to minimize misfit with data. In [2] MFU is discussed and a process to evaluate it for a given data set is provided, but no guidance is provided on how such results could be interpolated (or extrapolated) beyond that data set. In [3], MFU is characterized by augmenting uncertainty in a modeling assumption’s parameterization to account for discrepancy between data and predictions without perturbing the model assumption’s mathematical form. In [4], a data-driven modification to modeling assumptions is used to represent MFU. In [5]–[7] a statistical model for the time evolution of a correction to the assumption is inferred using time-series data. All the aforementioned approaches rely on access to observational data, which is often unavailable in scenarios where model predictions are most needed, e.g., to inform policy and engineering design decisions.

Through their parameterizations, MFU representations are constrained to respect known phenomenology and relevant physics, such as conserved quantities, with the aim of maximizing their extrapolative capability. To represent uncertainty, their parameterizations are assigned probability distributions, where relevant physical constraints and phenomenological behavior may also be encoded. The resulting MFU representation provides a range of plausible behavior for the modeled phenomenon. Because they are developed to respect known physics and phenomenology, MFU representations can be developed without access to calibration data, which is a common situation in the preliminary stages of model development and uncertainty analysis. However, if data is available, MFU representations can be informed using, e.g., Bayesian inference. Hierarchical frameworks that calibrate the hyperparameters of an MFU parameterization are common, so we demonstrate one such a framework here.

Namely, grouped Sobol’ indices measure the fraction of variance in a model prediction that can be attributed to each input factor in a model. Standard practice has focused on individual model parameters as input factors, but an input factor can also be defined as a group of parameters [8]–[10]. In the context of a group of parameters, Sobol’ indices measure the fraction of variance in a model prediction that can be attributed to all parameters in the group as well as to their interactions with each other.

Since multiple MFU representations may be valid in a given modeling scenario, we establish bounds for the differences in Sobol’ indices corresponding to varying MFU representations. Furthermore, we demonstrate how grouped VBSA can be used to analyze sensitivity to MFU when hierarchical probability densities are used, even after calibration, which often induces dependence between MFU and other model parameters.

Thus, our proposed methodology can be used to inform where to invest resources, e.g., more experimentation or computation, to drive down the most impactful sources of uncertainty on model predictions. Furthermore, following the spirit of [11], [12], it can be used to assess whether assumptions important to prediction use cases are also important in validation cases as a means of assessing the relevance of validation evidence. The key contributions of this work include:

• Presenting a novel approach to measure the impact of model-form uncertainty on model predictions using variance-based sensitivity analysis.

• Proving that grouped variance-based sensitivity analysis is robust to an MFU representation’s parameterization.

• Analyzing computation of Sobol’ indices for hierarchical distributions.

The remainder of the paper proceeds as follows. 2 presents the methodology of our proposed approach and provides the relevant background related to MFU representations (2.1) and their calibration (2.2) as well as grouped Sobol’ indices (2.3); an illustrative example is provided to build intuition. In 3.1, we establish robustness of the approach to varying MFU representations. We demonstrate our approach on a practical application problem of upscaled contaminant transport in 3. Final conclusions are presented in 4.

2 Background and methods↩︎

We first provide the mathematical framework for MFU representations along with a discussion of how this approach differs from alternative model discrepancy approaches [1] in 2.1. In 2.2 we discuss how data can be used to inform such representations through a hierarchical framework for calibration. Finally, 2.3 provides the notation and formulation of the grouped Sobol’ indices used to compute the sensitivity of the modeling assumption.

2.1 Model-form uncertainty representations↩︎

First introduced in [13] and applied in, e.g., [14]–[18], MFU representations augment the governing equations [14], [16]–[18] or replace an invalid assumption entirely [15].

Defining concepts mathematically, let a model governing the evolution of state variables \(v\) and model output mappings be denoted \[\begin{align} 0 &= \mathcal{R}_a(v), \tag{1}\\ d &= \mathcal{O}(v) + \epsilon_m, \quad \epsilon_m \sim \pi(\epsilon_m),\tag{2} \\ q &= \mathcal{Q}(v), \tag{3} \end{align}\] where 1 represents a set of governing equations in their residual form;  2 is the mapping of \(v\) to observable data \(d\) through an observation operator \(\mathcal{O}\), with measurement error \(\epsilon_m\) typically assumed to be a random variable with distribution \(\pi(\epsilon_m)\); and 3 is the mapping of \(v\) to a predicted quantity of interest (QoI), which may not be observable (e.g., a prediction of the future). In 1 , \(a\) is an aspect of the governing equations for which an assumption must be introduced, e.g., a coupling mechanism or the mathematical form for a constitutive equation. The assumption may not explicitly appear in the governing equations (e.g., neglecting a term), however, we explicitly denote its existence here to signify that a MFU representation will be introduced to characterize its impact on model predictions.

However, if the assumption is inaccurate, the governing equations will not provide a correct representation of reality; that is, the inaccurate assumption creates model-form error. The inaccuracy of the assumption will propagate through the governing equations to \(v\) and will be observed as a discrepancy between the predicted observable output \(\mathcal{O}(v)\) and \(d\) that cannot be explained merely by the presence of measurement error, a phenomenon known as model discrepancy. Mathematically, this amounts to the difference between observed data and model output not obeying the assumed probability distribution of the measurement error: \[d-\mathcal{O}(v) \not\sim \pi(\epsilon_m).\] Calibration of model parameters without accounting for model-form error can lead to biased, overconfident uncertainties, resulting in flawed representations of predictive uncertainties. Example [ex:poly] illustrates this phenomenon for a simple problem.

\[\begin{align} 0 &= \mathcal{R}_a(v),\\ d &= \mathcal{O}(v) + \epsilon_m + \delta_\phi, \quad \epsilon_m \sim \pi(\epsilon_m),\label{eq:ko} \\ q &= \mathcal{Q}(v). \end{align}\tag{6}\] While this approach is attractive as it is completely nonintrusive to the model, it is primarily useful in cases where the prediction quantity of interest is an observable output and predictions fall within the range of available calibration data, so that the statistical model is interpolating and not extrapolating [20]. Even interpolation is not guaranteed to be accurate, however [21]. Since the correction to the model output only appears in 6 , there is no way to propagate the discrepancy representation to an unobservable QoI (e.g., quantities representing the interior of modeled systems where sensors cannot be placed). Furthermore, even if prediction QoIs are observable outputs and a discrepancy representation can be calibrated, because it is purely statistical, its accuracy will be limited outside the regime where data is available.

In such cases, only embedded MFU representations, which address the model-form error at its source, can propagate to the QoIs. To this end, an MFU representation \(\xi_\gamma\) replaces or augments \(a\) to represent uncertainty in the assumption, thereby modifying the governing equations: \[\begin{align} 0 &= \mathcal{R}_{\xi_\gamma}(v), \tag{7}\\ d &= \mathcal{O}(v(\xi_\gamma)) + \epsilon_m, \quad \epsilon_m \sim \pi(\epsilon_m), \tag{8}\\ q &= \mathcal{Q}(v(\xi_\gamma)). \tag{9} \end{align}\] Here \(\gamma\) denotes the parameters of the MFU representation, and we denote the state variables \(v(\xi_\gamma)\) in 8 and 9 to emphasize the fact that \(\xi_\gamma\)’s inclusion in 7 propagates through \(v\) to the observable outputs \(\mathcal{O}(v)\) and prediction QoIs \(\mathcal{Q}(v)\). Hence, any information gained by calibration against observable data will be propagated to unobservable prediction QoIs.

Since the MFU representation \(\xi_\gamma\) is embedded in the governing equations and directly characterizes uncertainty in a modeling assumption, it can be constrained to respect physical properties of the modeled system and reflect what is known about the modeled phenomenon. Because more intuition can be brought to bear by directly characterizing uncertainty in the assumption, MFU representations can be developed without access to data, which can be especially useful early in the model development process. Furthermore, uncertainty in modeling assumptions can be encoded through the probability distributions assigned to the MFU parameters, which are defined to reflect what is known about plausible behavior for the modeled phenomenon.

MFU representations must be able to represent uncertainty arising from the fact that, often, no single mathematical model can perfectly capture a complex phenomenon that in reality depends upon quantities that aren’t resolved in the model (for example, neglected species in a chemistry model [14] or a population model [17]). No amount or quality of data can overcome such a structural inadequacy in the model; the goal must then be to use data to inform the appropriate degree of uncertainty in the model’s form.

2.2 Calibration of model-form uncertainty representations↩︎

Standard Bayesian paradigms for estimating uncertainties are often unsuitable for calibrating MFU parameterizations; in the limit of infinite data, the Bayesian posterior contracts about parameter values that best agree with the data [19], implying that uncertainty vanishes. In practice, such posterior contraction could allow one to grow increasingly confident in an incorrect mathematical model form. Thus, when informing MFU representations with data, one should employ an inference paradigm capable of avoiding posterior contraction when model-form error is present.

One approach to hyperparameterization, presented in [22], augments model parameters \(\theta\) with “discrepancy” random variables, which are stochastic: \(\theta + \delta, \; \delta \sim \pi_\delta(\cdot;\phi)\). The joint distribution of parameters and hyperparameters \(\pi(\theta,\phi)\) is then calibrated. Another common approach, as presented in [3], [13], [17], leverages hierarchical representations, where the MFU parameters are random variables described by parametric probability density functions (PDFs) whose hyperparameters are also uncertain (i.e., have corresponding PDFs). For example, if an MFU parameter \(\gamma\) is assumed to be normally distributed, its mean and standard deviation are considered uncertain with assumed PDFs as well. Denoting such hyperparameters as \(\phi\), the joint distribution of parameters and hyperparameters are defined hierarchically as \[\begin{align} \pi(\gamma,\phi) = \pi(\gamma|\phi)\pi(\phi).\label{eq:hierarchical95distribution} \end{align}\tag{10}\]

\[\begin{align} \pi(\phi|d) = \frac{\pi(d|\phi) \pi_0(\phi)}{\pi(d)}. \label{eq:hyperparameter95posterior} \end{align}\tag{11}\] An analytical expression for the likelihood with respect to hyperparameters \(\pi(d|\phi)\) is not typically available, so practical implementations perform a joint Bayesian calibration of parameters and hyperparameters: \[\begin{align} \pi(\gamma,\phi|d) = \frac{\pi(d|\gamma,\phi) \pi_0(\gamma,\phi)}{\pi(d)} =\frac{\pi(d|\gamma) \pi(\gamma|\phi)\pi(\phi)}{\pi(d)}, \end{align}\] where the model parameters are marginalized out in postprocessing.

Because the hyperparameters don’t appear directly in the model, forward propagation is carried out by sampling the hierarchical distribution 10 to generate samples of the MFU parameters. The prior pushforward of hyperparameters to MFU parameters uses the hyperprior, \(\pi_0(\phi)\). The posterior pushforward of hyperparameters to MFU parameters uses the posterior distribution in hyperparameters, \(\pi(\phi|d)\). An example of the hierarchical Bayesian approach to inform MFU is presented in Example [ex:hier].

2.3 Grouped Sobol’ indices↩︎

Specifically, Sobol’ indices are a global sensitivity analysis method that quantify the fraction of variance in a model output that can be attributed to an uncertain model input (or group of inputs). Sobol’ indices were originally defined in the context of the functional analysis of variance (ANOVA) decomposition [8], which we do not review here for the sake of brevity. Extensive discussion regarding Sobol’ indices and their computation can be found in, e.g., [8]–[10], [23]. Here we define grouped Sobol’ indices and methods for their computation.

First we define the mathematical notation for the remainder of the paper. Let \(\inpRV_i\) be a random variable with support \(\Omega_i \subset \mathbb{R}\), for \(i = 1, 2, \dots, \inpdim\). Consider \(f \in \mathcal{L}^2\) where \(f: \Omega \to \mathbb{R}\) for \(\Omega = \Omega_1 \times \Omega_2 \times \dots \times \Omega_{\inpdim}\). We denote the indices for parameter group \(\ug = \{i_1, i_2, \dots, i_k\} \subset \{1, 2, \dots, \inpdim\}\) and the group of inputs \(\binpRV_{\ug} = \left( \inpRV_{i_1}, \inpRV_{i_2}, \dots, \inpRV_{i_k}\right)\). We denote the complement of the group \(\binpRV_{\sim \ug} = \binpRV_{ \{1,2,\dots,\inpdim \} \setminus \{i_1, i_2, \dots, i_k\}}\). The joint probability density for \(\X_{\boldsymbol{u}}\) is denoted \(\pi(\x_{\boldsymbol{u}})\) (for notational simplicity, we omit the distribution subscripts in \(\pi_{\boldsymbol{u}}(\x_{\boldsymbol{u}})\)).

The parameters within a given group may be statistically dependent on each other. However, groups are independent of one another and other inputs. Hence, the joint probability density for \(\X\) can be written as the product of each group, e.g., \(\pi(\x) = \pi(\x_{\boldsymbol{u}})\pi(\x_{\sim \ug})\). As is common, we assume that the parameters are statistically independent to facilitate the derivation of the Sobol’ indices. However, we will relax this assumption in 2.3.3 to accommodate dependencies that arise as a result of calibration.

The grouped first-order Sobol’ index (or grouped main effect index) for parameter group \(\ub\) is defined as \[\begin{align} S_{\ub}^{g} = \frac{\Var_{\X_\ub}\left(\E_{\X_{\sim \ub}}\left[f(\X)|\X_\ub\right] \right)}{\Var(f(\X))}. \label{eq:grouped95main95effect} \end{align}\tag{14}\] This grouped index is equivalent to the closed Sobol’ index (often denoted \(S_\ub^{\text{clo}}\)) and measures the fraction of variance in \(f\) that can be attributed to each parameter in the group \(\ub\) and to the interactions between the parameters in the group. Note that if the cardinality of the group is \(1\), 14 reduces to the traditional first-order Sobol’ index for a single input.

The grouped total Sobol’ index (or grouped total effect index) for the group of parameters \(\ub\) is defined as \[\begin{align} T_\ub = \frac{\E_{\X_{\sim \ub}}\left[\Var_{\X_\ub}(f(\X)|\X_{\sim \ub})\right]}{\Var(f(\X))} = 1 - \frac{\Var_{\X_{\sim \ub}} \left( \E_{\X_{\ub}} \left[ f(\X) | \X_{\sim \ub}\right] \right) }{\Var(f(\X))}. \label{eq:grouped95total95effect} \end{align}\tag{15}\] The total Sobol’ index represents the contribution to the total variance due to \(\X_{\ub}\), including interactions with other model parameters. As with the grouped first-order Sobol’ index, if the cardinality of group \(\ub\) is \(1\), 15 reduces to the traditional total Sobol’ index for a single input.

2.3.1 Computation of grouped Sobol’ indices↩︎

We employ the estimators presented in [10], [24] for the main and total effect indices defined in 2.3. Derivations for the estimator of the numerator in the grouped first-order index can be found in [10], while Appendix A presents a similar derivation corresponding to the grouped total index. Let \(\pi({\boldsymbol{x}})\) be the probability density function (PDF) associated with \(\binpRV\), and let \({\boldsymbol{x}},{\boldsymbol{x}}' \sim \pi\) denote replicate random vectors whose PDF is \(\pi\), i.e., \(\pi({\boldsymbol{x}}) \equiv \pi({\boldsymbol{x}}')\). The numerator for the grouped first-order Sobol’ index is computed as \[\begin{align} \Var_{\X_\ub}\left(\E_{\X_{\sim \ub}}\left[f(\X)|\X_\ub\right]\right) &\approx \frac{1}{N}\sum_{i=1}^N f(\x^{(i)})\left( f(\x'^{(i)}_{\sim \ub}, \x^{(i)}_\ub) - f(\x'^{(i)}) \right), \quad \x^{(i)}, \x'^{(i)} \sim \pi(\x). \end{align}\] The numerator of the grouped total Sobol’ index is computed as \[\begin{align} \E_{\binpRV_{\sim \ug}} \left[ \V_{\binpRV_{\ug}} \left( f\left(\binpRV \right) | \binpRV_{\sim \ug} \right) \right] \approx \frac{1}{2N} \sum_{i=1}^N \left( f(\x^{(i)}) - f(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub) \right)^2, \quad \x^{(i)}, \x'^{(i)} \sim \pi(\x). \end{align}\] Finally, the total variance appearing in the denominators of the indices is computed as \[\begin{align} V(f(\X)) &\approx \frac{1}{2N}\sum_{i=1}^N \left[\left(f(\x^{(i)}) - \hat{\mu}_f \right)^2 + \left(f(\x'^{(i)}) - \hat{\mu}'_f \right)^2\right], \\ \hat{\mu}_f &\approx \frac{1}{N}\sum_{i=1}^N f(\x^{(i)}), \quad \hat{\mu}'_f \approx \frac{1}{N}\sum_{i=1}^N f(\x'^{(i)}). \nonumber \end{align}\]

2.3.2 Treatment of hierarchical distributions↩︎

Previous works have considered differing interpretations of Sobol’ indices in the context of hierarchical uncertainty characterizations. The work [25] formulates randomized Sobol indices, which compute an average (with respect to the hyperparameter distribution) Sobol’ index. This approach is based on the interpretation that model parameter uncertainty is reducible (epistemic), while randomness in the hyperparameters reflects variability across experiments (irreducible). Thus, the authors intend for the standard Sobol’ index to reflect sensitivity for a particular experiment, with the randomized indices representing aggregate sensitivity over all experiments. Similarly, the work of [26] proposes that for computer models with both reducible and irreducible sources of uncertainty, one should compute an expected Sobol’ index; however, alternative to [25], the expectation is with respect to the epistemic variables.

Here we wish to characterize the aggregate influence of MFU parameters and hyperparameters on model outputs. We thus define the MFU group in terms of the MFU parameters and their hyperparameters: \(\X_{\ug} = (\theta, \phi)\). All estimators presented in 2.3.1 and derived in 5 hold for samples drawn from a hierarchical posterior since the estimators only require sampling from the joint distribution over the group and make no assumptions about the structure of the distribution. This choice to represent the joint parameters and hyperparameters as a group ensures the resulting Sobol’ index reflects uncertainty associated with the MFU representation as a whole.

2.3.3 Treatment of dependent groups and inputs↩︎

Thus, the group of MFU parameters may be correlated to the other model parameters in the posterior distribution. This violates the assumption of statistical independence between the groups. Generalizations of Sobol’ indices in the presence of dependent inputs are considered in [27]–[30]. However, the traditional interpretation of Sobol’ indices attributing a relative variance contributing does not apply when the inputs are dependent, e.g., when they are correlated. Various alternative global sensitivity analysis approaches have been considered [31]–[33].

In this article, we focus on the total Sobol’ indices thanks to their alternative approximation theoretic interpretation that generalizes to the case with dependent inputs [34]. Specifically, the total Sobol’ index for the group of inputs \(\X_{\ug}\) corresponds to the relative approximation error if \(f\) is optimally approximated (in the \(L^2\) sense) by a function that only depends on \(\X_{\sim \ug}\). From this perspective, a small total index \(T_{\ug}\) implies that \(f\) can be approximated well without information from \(\X_{\ug}\), and hence \(\X_{\ug}\) is not influential.

The presence of input correlations typically results in smaller total indices because some information about parameter variability may be captured through correlations with other parameters. In the context of comparing two groups, model parameters and MFU parameters, a larger total index implies that information about variability from one group cannot be captured through correlations with the other group, and hence contributes more significantly to uncertainty in the function’s output. Therefore, if there are only two groups of parameters, it is possible to reason about the relative importance of the groups based on their total effects indices. However, if there are more than two groups of input parameters, this reasoning breaks down, since it is possible that a group has a large Sobol’ index because it is strongly correlated with a group that is highly influential, rather than because it is influential in its own right.

It would not be meaningful to directly compare the magnitudes of total effects indices between prior and posterior distributions—the dependence structure between the input groups has changed, and the posterior predictive variance is typically smaller than the prior predictive variance. However, total effect indices’ numerators can be meaningfully compared since the indices for prior and posterior both correspond to absolute function approximation errors with the same units. The reduction in magnitude from prior to posterior indicates the decrease in variability due to calibration. By comparing the total effect index numerators for model parameters and MFU parameters, we may understand how the relative importance of the two groups changed from prior to posterior. This is only possible, however, because we have split the parameters into two groups, so that there is no confounding dependence on a third group of parameters.  Example 2.3 demonstrates computation and interpretation of these indices for the MFU representation introduced in Example [ex:hier].

3 Results↩︎

In 3.1 we provide theoretical proof of robustness of grouped sensitivity analysis to MFU parameterization. In 3.2 we present numerical results demonstrating our method and corroborating our theoretical result. Our numerical studies are carried out in the context of a model for contaminant transport through a heterogeneous porous medium. Code to reproduce all numerical experiments in this section is publicly available at https://github.com/sandialabs/MFUQ-MFU-GSA.

3.1 Theoretical result: robustness of grouped sensitivity analysis to MFU parameterization↩︎

Here, we provide the assumptions under which the differences in Sobol’ indices for varying MFU representations will be bounded.

Let us consider two mathematical models, \(f:\Omega_{\X} \to \mathbb{R}\) and \(q:\Omega_{\Xtilde} \to \mathbb{R}\), whose difference lies in their MFU representation. The inputs for \(f\) are given as \(\X = \{\X_\ub, \X_\vb\}\), where \(\X_\ub\) represents the MFU parameters. Similarly, the inputs for \(q\) are given as \(\Xtilde = \{\Xtilde_\ubtilde, \X_\vb\}\), where \(\Xtilde_\ubtilde\) represents the alternative MFU parameters. Note that parameters other than MFU parameters are the same between both models and given as \(\X_\vb\). The Sobol’ indices for groups \(\ub\) and \(\ubtilde\) are \[\begin{align} S_\ub^g = \frac{\Var_{\Xu}\left(\E_{\X_\vb}[f(\X)|\Xu]\right)}{\Var(f(\X))} \quad \quad S_\ubtilde^g = \frac{\Var_{\Xtilde_\ubtilde}\left( \E_{\X_\vb}[q(\Xtilde) | \Xtilde_\ubtilde]\right)}{\Var\left(q(\Xtilde)\right)}, \end{align}\] and the grouped total Sobol’ indices are \[\begin{align} T_\ub = \frac{\E_{\Xu}\left[\Var_{\X_\vb}\left(f(\X)|\Xu\right)\right]}{\Var(f(\X))} \quad \quad T_\ubtilde = \frac{\E_{\Xtilde_\ubtilde}\left[ \Var_{\X_\vb}\left(q(\Xtilde) | \Xtilde_\ubtilde\right)\right]}{\Var\left(q(\Xtilde)\right)}. \end{align}\] We assume statistical independence between parameter groups, i.e., \(\pi(\X_\ub, \X_\vb)=\pi(\X_\ub)\pi(\X_\vb)\).

Assumptions (1) and (2) are motivated by the belief that two valid MFU representations should have a similar effect in the output statistics of the QoI prediction. Assumption (3) is for mathematical convenience and physically immaterial as it corresponds to a change in units (for which Sobol’ indices are invariant). Note that both \(f\) and \(q\) are rescaled in the same manner to maintain consistency in units. For a proof of Propositions [prop:bounds] and [prop:bounds952], see Appendix 6.

3.2 Numerical results: contaminant transport through a heterogeneous porous medium↩︎

Transport through porous media is governed by the advection-diffusion equation, where the flow field (fluid velocity) depends on the porous medium’s permeability and porosity fields through Darcy’s law [35]. Limitations in sensing technology make it impossible to perfectly characterize the heterogeneous properties of the subsurface. While it is possible to generate random instantiations of the subsurface and collect summary statistics to address missing information about the subsurface properties, this can be computationally prohibitive. Instead, it is common to assume statistical homogeneity in the subsurface (i.e., the statistical properties do not change throughout the domain) and to average the equations with the aim of predicting mean behavior.

\[\begin{align} \diffp[]{\meanc}{t}(x,t) + \mean{u}\diffp[]{\meanc}{x}(x,t) + \diffp[]{\mean{u'c'}}{x}(x,t) = \nu_p \diffp[2]{\mean{c}}{x}(x,t), \\ \mean{c}(0,t) = \mean{c}(L_x,t), \\ \mean{c}(x,0) = \exp\left(-\frac{(x - s)^2}{2\ell^2}\right), \end{align} \quad \begin{align} x \in [0,L_x], \end{align} \label{eq:averaged}\tag{16}\] where \(\mean{u}\) is the mean velocity, \(\partial_x \mean{u'c'}\) is the dispersion resulting from fluctuations in the velocity field, and \(\nu_p \partial_{xx} \meanc\) is the pore-scale diffusion of the contaminant. The initial condition is a Gaussian pulse with width \(\ell=0.1\) and mode \(s\). We assume periodic boundary conditions, which is valid provided the velocity fluctuations are homogeneous with correlation lengths small compared to \(L_x=4\).

An assumption about the mathematical form of the dispersion term must be introduced—this would be represented by the \(a\) in our definition of the governing equations in 1 . For heterogeneous porous media, significant fluctuations in the velocity lead to nonlocal transport of the contaminant, often called anomalous diffusion [36]–[38]. This effect is shown in 1, which compares the velocity fluctuations and evolution of the Gaussian initial condition in 16 for a homogeneous and a heterogeneous porous medium. The higher-magnitude and larger-lengthscale fluctuations of the velocity in the heterogeneous case transports the contaminant downstream at drastically different rates across the domain, resulting in larger concentrations downstream. In contrast to the homogeneous case, this effect is considered “anomalous” since its effect on the concentration does not mimic typical diffusion.

However, in general, the exact nature of the nonlocality in the dispersion is uncertain due to lack of detailed information about the subsurface. This leads to uncertainty in the assumed form of the dispersion.

3.2.1 Potential model-form uncertainty representations for dispersion↩︎

The most general form of such an operator is denoted as \[\begin{align} \diffp{\mean{u'c'}}{x} &\approx -\Lcal \mean{c}. \label{eq:general95linear} \end{align}\tag{17}\] Under the mild assumption that \(\Lcal\) is shift invariant (does not depend on absolute location, just relative distances), the Fourier modes are the eigenfunctions of the operator \(\Lcal\): \[\begin{align} \Lcal e^{i a_k} = \lambda_k e^{i a_k}, \quad a_k = \frac{2 \pi k }{L_x}, \quad k\in \mathbb{Z}. \end{align}\] 16 thus admits a Fourier series solution, and the action of \(\Lcal\) on \(\mean{c}\) can be entirely described in terms of its eigenvalues \(\lambda_k\). The eigenvalues may be complex-valued, with real parts influencing the magnitudes and imaginary parts influencing the phases of associated Fourier coefficients [39]. In this work, we consider three MFU representations for dispersion. First, a general linear operator parametrized by its eigenvalues, with as many eigenvalues as there are terms in the Fourier series solution: \[\begin{align} \Lcal_g [e^{ia_k}] &= \lambda_k e^{i a_k}, \quad k=1, \ldots, N_k. \label{oisaguyq} \end{align}\tag{18}\] Second, a Reisz fractional derivative operator, parameterized by its scaling coefficient \(\nu_m\) and fractional power \(\alpha\): \[\begin{align} \Lcal_{f} [e^{i a_k }] &= -\nu_m\diffp[\alpha]{}{x} [e^{i a_k} ] = - \nu_m (i a_k)^\alpha e^{i a_k}. \label{eq:fractional} \end{align}\tag{19}\] Finally, a complex operator inspired by the Reisz fractional derivative operator, parameterized by the complex and imaginary counterparts of its eigenvalues: \[\begin{align} \Lcal_{c} [e^{i a_k}] &= \left( - \nu_m^r (a_k)^{\alpha_r} + i \nu_m^i (a_k)^{\alpha_i} \right) e^{i a_k}. \label{eq:complex95fractional} \end{align}\tag{20}\]

The general linear operator given by 17 is the most flexible, since each eigenvalue can vary independent of the others. However, the number of parameters that must be constrained and informed by data could be large for a fine spatial discretization. A parametrization with dimensionality in the hundreds could easily result in computational intractability. Fractional derivative operators like 19 yield nonlocal behavior for fractional powers \(\alpha \in (1,2)\), making them an attractive option to represent MFU in dispersion for this application problem. Their simple parametrization in terms of the fractional power and scaling coefficient also makes them more tractable than the general linear operator; however, they may not be flexible enough to reproduce dispersion effects in all problems—for example, in practice the real and imaginary parts may need to vary independently to represent differing influences on diffusion and advection by dispersion. This inflexibility motivated the form of the complex operator inspired by the fractional derivative given in 20 . While having relatively few parameters, it is more flexible in the behavior of the real and imaginary parts of the eigenvalues.

For example, for the general linear operator \(\Re[\lambda_k] < 0\) ensures diffusive behavior in \(\mean{c}\), and \(\lambda_0=0\) conserves mass. For additional constraints that can be placed on the parametrization of the general linear operator, see [39]. For the fractional derivative and complex fractional operators, all scaling coefficients must be positive to ensure diffusive behavior (for real parts) and advection downstream (for imaginary parts). Additionally, the fractional powers must remain in the range \((1,2)\) to express nonlocality of the dispersion. A benefit specific to the fractional derivative operator is that it preserves positivity of \(\mean{c}\) by construction, which is not as easily imposed on the other operators.

3.2.2 Model-form uncertainty propagation without data↩︎

In all following discussions, the quantity of interest is the mean concentration \(\mean{c}\) at the outflow boundary at time \(t=1.5\). We first characterize uncertainty in the model parameters bulk velocity \(\mean{u}\), pore-scale diffusion coefficient \(\nu_p\), and contaminant source location \(s\): \[\begin{align} \mean{u} &\sim 1 \cdot \left( 1 + \mathcal{U}[ -0.1, 0.1]\right),\\ \nu_p &\sim 0.01 \cdot \left(1 + \mathcal{U}[-0.2, 0.2 ] \right), \\ s &\sim \mathcal{U}[0.2, 1.5]. \end{align} \label{eq:other95param95distributions}\tag{21}\] The distributions for \(\mean{u}\) and \(\nu_p\) represent a \(\pm10\) and \(20\) percent uniform uncertainty, respectively, about a nominal value. The distribution for \(s\) reflects the assumption that the source location is in the upstream section of the computational domain without being too close to the upstream boundary.

Encoding physical properties (see 3.2.1) into the probabilistic representation of the MFU parameters of 19 provides the following assumed distributions: \[\begin{align} \tag{22} \nu_m &\sim \scriptU[0.05, 0.15], \\ \tag{23} \alpha &\sim \text{Triangular}[1,2,1.5], \end{align}\] where 22 is based on the assumptions that dispersion dominates pore-scale diffusion (and thus the scaling coefficient \(\nu_m\) should be larger than \(\nu_p\)) and that dispersion is not so great that it completely diffuses the solution by the time it has been transported the length of the domain. The expression given by 23 is based on the knowledge that the fractional power should fall in the range \([1,2]\), with the mode set at \(1.5\) to weakly concentrate fractional powers away from the boundary values, which result in pure advection or pure diffusion.

Forward propagation of these uncertainties is shown in 2, where 2 (a) displays 10 sample evolutions of \(\mean{c}\) over the computational domain. Note the varying degrees of nonlocality in the concentration field, as evidenced by some samples being right-skewed with a heavier tail downstream. A histogram of our QoI, the outflow concentration at \(t=1.5\), is shown in 2 (b). Note that all samples are positive due to the fractional derivative MFU representation enforcing positivity of concentrations by construction.

3.2.3 Verification of robustness of Sobol’ indices for different MFU parametrizations↩︎

Since the fractional derivative is a specialization of the general linear operator, we can derive a parametrization of \(\Lcal_g\) that produces very similar output statistics to those induced by \(\Lcal_f\), which enables us to evaluate the theoretical results numerically. Following the notation of 3.1, we compare the grouped Sobol’ index for the fractional derivative MFU parameters \(\X_\ub=[\nu_m,\alpha]\) and the general linear operator MFU parameters \(\Xtilde_\ubtilde=\lb\) relative to the other uncertain model parameters in 16 , \(\X_{\vg} = [\langle u\rangle, \nu_p, s]\). We additionally denote the QoI obtained with the fractional derivative MFU representation as \(f(\X)\) and with the general linear operator as \(g(\Xtilde)\). We adopt the probability distributions for \(\X_\ub\) and \(\X_{\vg}\) discussed in 3.2.2.

DCI is a measure-theoretic framework for solving stochastic inverse problems that identifies parameter distributions whose predicted outputs reproduce distribution on the observed data. Therefore, DCI is leveraged to determine a distribution on the eigenvalues \(\lambda_k\) that results in the the distribution—and hence statistics–of \(g(\X)\) being consistent with\(f(\X)\) (the QoI corresponding to the fractional derivative representation). More details on this procedure can be found in 6.

50 replicate Sobol’ indices are computed using 5e4 independent samples. For fairness in comparison, the same set of random samples of the auxiliary parameters \(\langle u\rangle, \nu_p\) and \(s\) are used to estimate indices across the two MFU representations. For each replicate sample, the QoI is rescaled so that its variance for the fractional derivative MFU representation is equal to \(1\). The mean variance over replicates of the QoI for the general linear operator representation is \(1.02\) with a standard deviation of \(0.01\), indicating an approximately 2% difference in the variance from the fractional derivative parameterization.

The Sobol’ index values are qualitatively similar; the conclusions would be identical for a ranking problem (ranking the most important sources of uncertainty).

As discussed above, the norm of differences in total variances \(\errtwo\) is 0.02 on average. We compute the norm of differences in conditional variances \(\errone\) with quadrature when possible and Monte Carlo sampling otherwise, again computing 50 replicate samples. On average, \(\errone=0.12\) with a standard deviation of \(0.02\). The resulting bounds on the differences between the main and total effects indices for the MFU representation parameters are therefore approximately \(0.16\) and \(0.14\), respectively. The results in 3 indicate the actual difference between the main and total effects indices to be \(0.008\) and \(0.025\) on average, which could indicate the bounds derived in 3.1 are not tight.

Future work could aim to derive tighter bounds; however, for the purpose of this work, this numerical example is sufficient to confirm the intuitive understanding that MFU parametrizations resulting in similar QoI statistics should have similar Sobol’ indices.

3.2.4 Calibrating model-form uncertainty with data↩︎

This work serves as an application of the methods demonstrated on a simple polynomial example in 2, applied to the more complex application problem of contaminant transport. We found in preliminary studies that the fractional derivative MFU representation defined in 19 was not adequately flexible to capture MFU for the calibration case, so we use the complex fractional MFU representation defined in 20 in these results.

Prior densities↩︎

Since each of the model parameters is known to be positive with a specified nominal value, a log-normal prior density is assumed where the mode falls at the nominal value, and 95% probability mass falls below 120% of the nominal value. Nominal values are 1, 0.01, and 1 for \(\mean{u}, \nu_p,\) and \(s\) respectively.

We assign the same hyperpriors for both scaling parameters as there is no reason to distinguish between the real and imaginary parts of the eigenvalues a priori. Recall that the PDF for a log-normal random variable \(x\) is defined as \[\begin{align} \pi(x) &= \frac{1}{x\sigma\sqrt{2\pi}}\exp\left( -\frac{(\ln(x)-\mu)^2}{2\sigma^2} \right). \end{align}\] Therefore, hyperpriors must be defined for \(\mu\) and \(\sigma\).

Namely, we anticipate the scaling parameters fall in the range \([0.1,0.5]\), with weaker confidence in the lower vs. the upper bound. We express this belief by stating that \(\mathbb{P}(\nu_m\leq 0.1)\approx 0.1\) and \(\mathbb{P}(\nu_m\leq 0.5)\approx 0.99\). We map this belief about plausible values for the scaling parameters into nominal values for the hyperparameters \(\mu\) and \(\sigma\) using the quantile function for log-normal random variables, which returns a value \(x\) associated with a probability level \(p\): \[\begin{align} Q(p) = \exp\left( \mu + \sigma\Phi^{-1}(p)\right), \end{align}\] where \(\Phi\) is the CDF of the standard normal distribution. Substituting \(Q(0.1)=0.1\) and \(Q(0.99)=0.5\) into the equation above, we arrive at two equations to solve for the two unknowns, \(\mu\) and \(\sigma\), which we will use as our nominal values \(\mu_n\), \(\sigma_n\). Since \(\mu\) needn’t be positive, it is assumed to be a normal random variable with mean \(\mu_n\) and standard deviation \(0.5|\mu_n|\). Since \(\sigma\) must be positive, it is assumed to be log-normal with mode at \(\sigma_n\) and 99% probability below \(1.5\sigma_n\).

As in 3.2.2, we assume the fractional powers follow a triangular density on the range \([1,2]\), but the mode is inferred as a hyperparameter. With no intuition about the mode beyond that it should fall within \([1,2]\), we assume a uniform prior on this range.

Likelihood↩︎

\[\begin{align} \diffp[]{\meanc}{t}(x,t) + \diffp[]{\mean{u'c'}}{x}(x,t) = \nu_p \diffp[2]{\mean{c}}{x}(x,t) + \Lcal\mean{c}, \\ \mean{c}(0,t) = \mean{c}(L_x,t), \\ \mean{c}(x,0) = \exp\left(-\frac{(x - s)^2}{2\ell^2}\right), \end{align} \quad \begin{align} x \in [0,L_x], \end{align} \label{eq:likelihood95model}\tag{24}\] with \(\mean{u}=1.05\), \(\nu_p=0.0095\), \(s=0.9\), and \(\ell, L_x\) defined identically to 16 . The linear operator is parameterized by its eigenvalues, where high-fidelity eigenvalues are computed via Monte Carlo sampling and ensemble averaging of the detailed governing equations, as described in [15]. The problem is discretized using a uniform grid with \(N_x=512\). Since the averaged equation admits a Fourier series solution, it can be evaluated at any point in time.

The calibration data is thus collected at \(x=1.4\) with observations at \(\mathbf{t}=[0.01, 0.02, \ldots, 0.2]\) (20 uniformly-spaced time observations up to \(t=0.2\)). Because concentration data must always be positive, measurement error is assumed to be multiplicative with log-normal measurement error: \[\begin{align} d_i = \epsilon_m^{(i)} \cdot \mean{c}(1.4, t_i). \end{align}\] We thus pose the likelihood in terms of the log-transformed concentrations, i.e., \[\begin{align} \log(d_i) = \log(\mean{c}(1.4, t_i)) + \log(\epsilon_m^{(i)}), \quad \log(\epsilon_m^{(i}))\sim\Ncal(0, (10^{-2})^2). \end{align}\] Note that this formulation of the measurement uncertainty is posed such that the median of the untransformed error distribution is 1.W

This is done in a rejection-sampling fashion, where if negative concentrations are observed for a given ordered pair of parameter values, a likelihood of \(0\) is assigned. Since it would be computationally intractable to check for positivity over all times, we evaluate the model and check for negative concentrations at all spatial locations for times \(0.1, 0.5, 1, 1.5,\) and \(2\). Note that since this check is only performed for a small number of time points during calibration, this approach can’t penalize parameters that produce negative concentrations at other times. Therefore, forward uncertainty propagation is not guaranteed to result in positive concentrations for all model outputs.

Note that the times and location of calibration data differ significantly from the prediction QoI’s, so proper calibration of uncertainties is critical to meaningfully extrapolate.

Inference results↩︎

The source location for the contaminant, \(s\), exhibits significant bias in the posterior relative to its true value, indicating that for this example, hierarchical calibration leveraging an MFU representation did not completely mitigate posterior bias and overconfidence. We expand upon challenges associated with hierarchical inference that impact its application to models leveraging an MFU representation in [40].

Additionally, although the support of the prior and posterior pushforward densities does not visually differ significantly, the prior pushfoward has heavier tails that result in significantly higher variance (\(10^{-3}\) and \(3\cdot 10^{-4}\) for prior and posterior, respectively).

This is because the complex fractional MFU representation doesn’t preserve positivity of concentration by construction. However, due to the Bayesian calibration, the posterior pushforward has a lower incidence of negative values—2% for posterior vs. 7% for prior.

3.2.5 Forward propagation of MFU post-calibration↩︎

However, as discussed in 2.3.3, it is still meaningful to compute the total effects indices for the two groups of parameters, as is depicted in 6 (b). Clearly, uncertainty in the dispersion assumption dominates uncertainty in the QoI post-calibration.

The magnitude of the numerators for the posterior pushforward to the QoI has dropped relative to the prior pushforward. For both the prior and posterior, the MFU parameters appear to be more influential, although the ratio of the numerators between the groups becomes more pronounced post-calibration, indicating increased significance of MFU relative to model parameter uncertainty. In the case of the posterior sensitivities, the ratio between MFU and model parameters is so great that the importance of model parameters could be considered negligible, indicating that further data collection or model improvement efforts should be focused on the dispersion assumption.

4 Conclusions↩︎

We leverage model-form uncertainty (MFU) representations to represent uncertain modeling assumptions while grouped variance-based sensitivity analysis (VBSA) is used to quantify the impact of MFU on model predictions, which may be extrapolative. This approach enables model assumptions and other sources of uncertainty to be holistically ranked in order of importance to model predictions, thereby facilitating more quantitative prioritization of modeling and data acquisition efforts towards the aspects of the model that most impact predictions.

Thus, MFU representations that produce “similar” effects on model outputs will have “similar” Sobol’ index values. We numerically demonstrate the robustness of the approach on an upscaled contaminant transport problem and find that for two MFU representations of dispersion, the difference in Sobol’ indices is negligible. Future work on the robustness of Sobol’ indices to MFU representation will focus on the derivation of tighter bounds.

In the case where data is used to inform the representation via Bayesian inference, we show how VBSA can still be applied and interpreted despite violation of the common assumption that input distributions are statistically independent. We present numerical demonstrations in the form of a simple polynomial example as well as for the upscaled contaminant transport problem. While the approach was applied to the hierarchical MFU formulation presented in 10 , it is amenable to other formulations, such as the one posed by [22], so long as the formulation yields an MFU parameterization.

This enables assessment of assumption importance to extrapolative model predictions and across a range of modeling use cases. Future work could exploit this capability to measure the importance of assumptions to model outputs used in validation contexts to ensure that validation tests exhibit similar sensitivity to assumptions as those observed for target prediction use cases. Application of our approach to assess validation test relevance would follow in the spirit of previous works which employ sensitivity analysis to assess the utility of validation tests [11], [12]—however, this work allows for the novel incorporation of MFU into such assessments.

A general approach to develop an MFU representation does not currently exist, and given the problem-specific nature of model assumptions, a general approach may be infeasible. However, methods that enable rapid prototyping of MFU representations with minimal intrusion on simulation codes would significantly increase ease of adoption of our approach. Additionally, pick-freeze methods to compute Sobol’ indices are computationally costly, often requiring \(\geq 10^5\) evaluations, with costs that scale linearly with the number of inputs/groups. Methods that mitigate this cost are needed to increase tractability for computationally expensive models. Given-data [41]–[43] and multifidelity methods [44] for Sobol’ index computation are promising avenues of investigation in this regard. Finally, methods are needed for joint calibration of model and MFU parameters that mitigate confounding effects between the parameter groups.

Acknowledgements↩︎

We thank Luis Damiano for helpful conversations on hierarchical Bayesian problems and Tian Yu Yen for helpful conversations on data-consistent inversion. This work was supported by the Laboratory Directed Research and Development program (Project 233072) at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc. for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC. The employee owns all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan.

5 Grouped Sobol’ index estimator derivations↩︎

Let us first derive the estimator for the numerator for the grouped total Sobol’ index, \[\E_{\binpRV_{\sim \ug}} \left[ \V_{\binpRV_{\ug}} \left( f\left(\binpRV \right) | \binpRV_{\sim \ug} \right) \right].\] Let \(\pi({\boldsymbol{x}})\) be the probability density function (PDF) associated with \(\binpRV\), which is assumed to exist; similarly, \(\pi({\boldsymbol{x}}_{\ug})\) and \(\pi({\boldsymbol{x}}_{\sim \ug})\) are the joint PDFs associated with \(\binpRV_{\ug}\) and \(\binpRV_{\sim \ug}\), respectively. Furthermore, we let \(\pi({\boldsymbol{x}}) \equiv \pi({\boldsymbol{x}}')\), and use prime notation to denote a replicate set of samples, e.g. \({\boldsymbol{x}} \sim \pi({\boldsymbol{x}}) \implies {\boldsymbol{x}}' \sim \pi({\boldsymbol{x}})\).

Next, consider that \[\begin{align} \V_{\binpRV_{\ug}} \left( f\left(\binpRV \right) | \binpRV_{\sim \ug} \right) = \E_{\binpRV_{\ug}} \left[f^2\left( \binpRV \right) | \binpRV_{\sim \ug} \right] - \E^2_{\binpRV_{\ug}} \left[f\left( \binpRV \right) | \binpRV_{\sim \ug} \right]. \end{align}\] Therefore, \[\begin{align} \label{eq:g95tot95derv} \E_{\binpRV_{\sim \ug}} \left[ \V_{\binpRV_{\ug}} \left( f\left(\binpRV \right) | \binpRV_{\sim \ug} \right) \right] = \E_{\binpRV_{\sim \ug}} \left[ \E_{\binpRV_{\ug}} \left[f^2\left( \binpRV \right) | \binpRV_{\sim \ug} \right] \right] - \E_{\binpRV_{\sim \ug}} \left[ \E^2_{\binpRV_{\ug}} \left[f\left( \binpRV \right) | \binpRV_{\sim \ug} \right] \right]. \end{align}\tag{25}\] Now consider that the second term of 25 can be expanded as \[\begin{align} \nonumber &&\E_{\binpRV_{\sim \ug}} \left[ \E_{\binpRV_{\ug}} \left[f\left( \binpRV \right) | \binpRV_{\sim \ug} \right] \E_{\binpRV_{\ug}} \left[f\left( \binpRV \right) | \binpRV_{\sim \ug} \right] \right] \\\nonumber = && \E_{\binpRV_{\sim \ug}} \left[ \int \int f\left( {\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}_{\ug}\right) f\left( {\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}'_{\ug}\right) \pi({\boldsymbol{x}}_{\ug}) d{\boldsymbol{x}}_{\ug} \pi({\boldsymbol{x}}'_{\ug}) d{\boldsymbol{x}}'_{\ug} \right] \\ = && \int \left( \int \int f\left( {\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}_{\ug}\right) f\left( {\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}'_{\ug}\right) \pi({\boldsymbol{x}}_{\ug}) d{\boldsymbol{x}}_{\ug} \pi({\boldsymbol{x}}'_{\ug}) d{\boldsymbol{x}}'_{\ug} \right) \pi({\boldsymbol{x}}_{\sim \ug})d{\boldsymbol{x}}_{\sim \ug}\\\label{eq:term2} =&& \int \int f\left( {\boldsymbol{x}} \right) f\left( {\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}'_{\ug}\right) \pi({\boldsymbol{x}}) d{\boldsymbol{x}} \pi({\boldsymbol{x}}'_{\ug})d{\boldsymbol{x}}'_{\ug}. \end{align}\tag{26}\] Now consider the first term in 25 , can be expanded as \[\begin{align} \nonumber \E_{\binpRV_{\sim \ug}} \left[ \E_{\binpRV_{\ug}} \left[f^2\left( \binpRV \right) | \binpRV_{\sim \ug} \right] \right] &=& \E_{\binpRV_{\sim \ug}} \left[ \int f^2({\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}_{\ug})\pi({\boldsymbol{x}}_{\ug})d{\boldsymbol{x}}_{\ug} \right]\\\nonumber &=& \int \int f^2({\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}_{\ug})\pi({\boldsymbol{x}}_{\ug})d{\boldsymbol{x}}_{\ug} \pi({\boldsymbol{x}}_{\sim \ug})d{\boldsymbol{x}}_{\sim \ug} \\\label{eq:term1} &=& \int f^2({\boldsymbol{x}})\pi({\boldsymbol{x}})d{\boldsymbol{x}}. \end{align}\tag{27}\] We can then add a constant into 27 as \[\begin{align} \label{eq:term1952} \int \int f^2({\boldsymbol{x}})\pi({\boldsymbol{x}})d{\boldsymbol{x}} \pi({\boldsymbol{x}}'_{\ug})d{\boldsymbol{x}}'_{\ug}. \end{align}\tag{28}\] Finally, we add 26 and 28 to represent the numerator of the grouped total Sobol’ index as \[\begin{align} \nonumber \E_{\binpRV_{\sim \ug}} \left[ \V_{\binpRV_{\ug}} \left( f\left(\binpRV \right) | \binpRV_{\sim \ug} \right) \right] &=& \int \int \left( f^2\left({\boldsymbol{x}} \right) - f\left( {\boldsymbol{x}}\right) f\left( {\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}'_{\ug}\right) \right) \pi({\boldsymbol{x}})d{\boldsymbol{x}} \pi({\boldsymbol{x}}'_{\ug})d{\boldsymbol{x}}'_{\ug} \\\label{eq:total95num95est} &\approx& \frac{1}{2N} \sum_{i=1}^N \left( f(\x^{(i)}) - f(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub) \right)^2, \quad \x^{(i)}, \x'^{(i)} \sim \pi(\x). \end{align}\tag{29}\] Note that 29 results from the fact that \[\begin{align} \frac{1}{2}\left( f(\x^{(i)}) - f(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub) \right)^2 = \frac{1}{2} \left( f^2(\x^{(i)}) - 2f(\x^{(i)})f(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub) + f^2(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub) \right), \end{align}\] where \(\frac{1}{2}\left(f^2(\x^{(i)}) + f^2(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub) \right)\) is an estimator for \(f^2({\boldsymbol{x}})\), and \(f(\x^{(i)})f(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub)\) an estimator for \(f({\boldsymbol{x}})f({\boldsymbol{x}}_{\sim \ug}, {\boldsymbol{x}}'_{\boldsymbol{u}})\). Furthermore, from 29 one can understand the term pick-and-freeze as \(f(\x^{(i)}_{\sim \ub},\x'^{(i)}_\ub)\) differs from \(f(\x^{(i)}) = f(\x_{\sim \ug}^{(i)}, \x_{\ug}^{(i)})\) in that the group of inputs \(\ug\) are from the replicate samples \({\x}'\) rather than \({\x}\); we’ve “picked” the group \(\sim \ug\) to be “frozen”, i.e., from the same samples \({\x}\).

A similar derivation produces the estimator for the numerator of the grouped first-order index given in 14 : \[\begin{align} \Var_{\X_\ub}\left(\E_{\X_{\sim \ub}}\left[f(\X)|\X_\ub\right]\right) &\approx \frac{1}{N}\sum_{i=1}^N f(\x^{(i)})\left( f(\x'^{(i)}_{\sim \ub}, \x^{(i)}_\ub) - f(\x'^{(i)}) \right), \quad \x^{(i)}, \x'^{(i)} \sim \pi(\x). \end{align}\] See [10] for details. Finally, we use the two replicate sample sets \(\left\{\x^{(i)}, f(\x^{(i)})\right\}_{i=1}^N\) and \(\left\{\x'^{(i)}, f(\x'^{(i)})\right\}_{i=1}^N\) to estimate the total variance used in the denominators of the Sobol’ index estimators as \[\begin{align} \hat{\mu}_f &\approx \frac{1}{N}\sum_{i=1}^N f(\x^{(i)}), \nonumber \\ \hat{\mu}'_f &\approx \frac{1}{N}\sum_{i=1}^N f(\x'^{(i)}), \nonumber \\ V(f(\X)) &\approx \frac{1}{2N}\sum_{i=1}^N \left(f(\x^{(i)}) - \hat{\mu}_f \right)^2 + \left(f(\x'^{(i)}) - \hat{\mu}'_f \right)^2. \end{align}\]

6 Proofs of Propositions [prop:bounds] and [prop:bounds952]↩︎

Here, we prove Propositions [prop:bounds] and [prop:bounds952].

We now prove Proposition [prop:bounds] as follows:

Next, we prove Proposition [prop:bounds952] as follows:

7 Data-consistent inversion↩︎

Here, details are provided on the use of data-consistent inversion (DCI) to formulate the sampling distribution on the linear operator MFU parameters – needed to evaluate robustness (see, 2.3.1) for the contaminant transport problem. Given an initial distribution on parameters \(\pi_0(\Xtu)\), the DCI solution is defined as \[\begin{align} \pi_{update}(\Xtu) = \pi_0(\Xtu) \frac{\pi_{target}(q(\Xtu))}{\pi_{predict}(q(\Xtu))}, \label{eq:dataconsistent95update} \end{align}\tag{34}\] where \(\pi_{target}\) is the target probability density of the model output using the fractional derivative MFU representation, and \(\pi_{predict}\) is the probability density of the pushforward of the initial density using the general linear operator representation. The solution to 34 , \(\pi_{update}(\Xtu)\), has the property that, for any set \(A\) in the Borel \(\sigma\)-algebra defined on the output space, \(\mathcal{B}_D\), \[\begin{align} \mathbb{P}_{update}(q^{-1}(A)) \equiv \int_{q^{-1}(A)} \pi_{update}(\Xtu) \d \Xtu = \int_A \pi_{target}( q ) \d q \equiv \mathbb{P}_{target}(A), \quad \forall A\in \mathcal{B}_D. \end{align}\] That is, the pushforward of the updated density is consistent with the target density. However, in practice, it is necessary to approximate the update using, e.g., rejection sampling, as described in Algorithm 2 of [45]. As a result, the output statistics produced using \(\pi_{update}(\Xtu)\) will be similar (but not identical) to the output statistics produced using the fractional derivative MFU representation.

Implementation of DCI requires two key elements: approximating \(\pi_{target}\) and determining \(\pi_0(\Xtu)\). To approximate \(\pi_{target}\), we construct a kernel-density estimate (KDE) using \(10^3\) samples of \(f(\X)\) with \(\X_{\vg}\) set to their mean values and \(\X_\ub\) sampled according to [eq:nu_dist,eq:alpha_dist]. We defined \(\pi_0(\Xtu)\) using a multivariate normal approximation of the fractional derivative eigenvalues as follows. Given samples of \(\nu_m\) and \(\alpha\), eigenvalues of the fractional derivative operator are computed according to 19 . The distribution of fractional derivative eigenvalues can then be approximated by a multivariate normal distribution using the sample mean and covariance of the samples, where the multivariate normal distribution is constructed with respect to the following unraveled and transformed values of the eigenvalues: \(\big[\log(-\Re[\lambda_1]), \log(-\Re[\lambda_2]), \ldots, \log(\Im[\lambda_1]), \log(\Im[\lambda_2]), \ldots]\big]\). Thus, we generate samples from \(\pi_0(\Xtu)\), as follows:

• Sample \(R_1, R_2, \ldots, I_1, I_2, \ldots \sim \mathcal{N}(\hat{\mu}, \hat{\Sigma})\).

• \(\lambda_k = -\exp(R_k) + i \exp(I_k)\).

Histograms of the QoI from the target distribution and before and after DCI are shown in 8. Note that the histograms of \(\pi_{target}\) and the pushforward of the data-consistent update, \(\pi_{update,Q}\) are very similar. Additionally, the updated density is computed with the other sources of uncertainty in the model (the model parameters) set at their mean values; as they vary during VBSA, we expect there to be differences between the model outputs statistics for the two MFU representations due to sampling the model parameters. The goal of this DCI procedure is to produce similar output distributions for the purpose of numerically verifying our robustness bounds, so this is desirable.

References↩︎