Causality and Bayesian network PDEs for multiscale representations of porous media

TL;DR

This paper develops a Bayesian network-based multiscale modeling approach for porous media, incorporating causal relationships among microscopic properties to improve predictions of macroscopic behavior and uncertainty quantification.

Contribution

It introduces a systematic method to embed causal relationships into stochastic multiscale models using Bayesian networks, enhancing the accuracy of macroscopic predictions.

Findings

Causal relationships significantly influence macroscopic QoIs.

Including correlations alters uncertainty propagation results.

The approach improves physical interpretability of multiscale models.

Abstract

Microscopic (pore-scale) properties of porous media affect and often determine their macroscopic (continuum- or Darcy-scale) counterparts. Understanding the relationship between processes on these two scales is essential to both the derivation of macroscopic models of, e.g., transport phenomena in natural porous media, and the design of novel materials, e.g., for energy storage. Most microscopic properties exhibit complex statistical correlations and geometric constraints, which presents challenges for the estimation of macroscopic quantities of interest (QoIs), e.g., in the context of global sensitivity analysis (GSA) of macroscopic QoIs with respect to microscopic material properties. We present a systematic way of building correlations into stochastic multiscale models through Bayesian networks. This allows us to construct the joint probability density function (PDF) of model parameters through causal relationships that emulate engineering processes, e.g., the design of hierarchical nanoporous materials. Such PDFs also serve as input for the forward propagation of parametric uncertainty; our findings indicate that the inclusion of causal relationships impacts predictions of macroscopic QoIs. To assess the impact of correlations and causal relationships between microscopic parameters on macroscopic material properties, we use a moment-independent GSA based on the differential mutual information. Our GSA accounts for the correlated inputs and complex non-Gaussian QoIs. The global sensitivity indices are used to rank the effect of uncertainty in microscopic parameters on macroscopic QoIs, to quantify the impact of causality on the multiscale model's predictions, and to provide physical interpretations of these results for hierarchical nanoporous materials.

Tables

Table 1. Table 1 . Narrow hyperparameter ranges that allow for comparison with results in [ 47 ] .

	$p=R$ ($\mathrm{nm}\text{/}$)	$p=\theta$ ($\mathrm{rad}\text{/}$)	$p=d$ ($\mathrm{nm}\text{/}$)	$p=l$ ($\mathrm{nm}\text{/}$)
(maximum) $p_{+}$	$60$	$0.7$	$8$	$18$
(minimum) $p_{-}$	$10$	$0.07$	$4$	$8$

Table 2. Table 2 . Physical hyperparameter ranges.

	$p=R$ ($\mathrm{nm}\text{/}$)	$p=\theta$ ($\mathrm{rad}\text{/}$)	$p=d$ ($\mathrm{nm}\text{/}$)	$p=l$ ($\mathrm{nm}\text{/}$)
(maximum) $p_{+}$	$60$	$0.4\pi$	$60$	$60$
(minimum) $p_{-}$	$10$	$0.05\pi$	$5$	$1$

Table 3. Table 3 . Results for a two-way nonparametric Cramér test ( [ 4 ] ) on the hypothesis of equality of the empirical distirubtions for comparable variables displayed in Figures   7(a) , 7(b) , 8(a) and  8(b) each based on 2000 2000 2000 values sampled using a gPCE.

Variable	Cramér-statistic	Critical value	Conf. interval	$p$-value	Result	Figures
$D_L$	$10.16$	$0.3116$	$0.95$		reject	\rdelim}30mm[7(a) vs. 7(b)]
$D_T$	$1.154$	$0.1049$	$0.95$		reject
${\gamma }_{\mathrm{eff}}$	$0.2232$	$0.3835$	$0.95$	$0.154$	accept
$(D_L,D_T)$	$10.2$	$0.3154$	$0.95$		reject	\rdelim}30mm[8(a) vs. 8(b)]
$(D_T,{\gamma }_{\mathrm{eff}})$	$1.3$	$0.3854$	$0.95$		reject
$({\gamma }_{\mathrm{eff}},D_L)$	$7.525$	$0.4533$	$0.95$		reject

Table 4. Table 4 . The global sensitivity index S 𝑆 S in 31 based on the mutual information quantifies the effect of pore-scale uncertainties in terms of how additional knowledge of the input Θ ∈ 𝚯 Θ 𝚯 \Theta\in\boldsymbol{\Theta} reduces uncertainty in our prediction of the macroscopic Darcy-scale variables D L subscript 𝐷 𝐿 {D_{L}} , D T subscript 𝐷 𝑇 {D_{T}} , and γ eff subscript 𝛾 eff \gamma_{\mathrm{eff}} . Below, the estimator S ^ ^ 𝑆 \widehat{S} in 33 and the associated ranking r ^ ^ 𝑟 \widehat{r} in 35 is given for the probabilistic models P 1 subscript 𝑃 1 P_{1} in 20 over both the hyperparameter ranges in Table   1 and Table   2 (cf. graphical representation of rankings in Figure   10 ).

		Mutual information 33				Ranking effects 35			Hyperparameters
$\Theta$		${\hat{S}}_{\Theta }\left(D_L\right)$	${\hat{S}}_{\Theta }\left(D_T\right)$	${\hat{S}}_{\Theta }\left({\gamma }_{\mathrm{eff}}\right)$		${\hat{r}}_{\Theta }\left(D_L\right)$	${\hat{r}}_{\Theta }\left(D_T\right)$	${\hat{r}}_{\Theta }\left({\gamma }_{\mathrm{eff}}\right)$
${\Theta }_R$		$0.0419$	$0.0424$	$0.8955$		$0.0724$	$0.1177$	$0.8509$	\rdelim}40mm[Table 1]
${\Theta }_{\theta }$		$0.5074$	$0.2049$	$0.0366$		$0.8770$	$0.5692$	$0.0348$
${\Theta }_d$		$0.0150$	$0.0878$	$0.0271$		$0.0259$	$0.2438$	$0.0258$
${\Theta }_l$		$0.0143$	$0.0249$	$0.0932$		$0.0247$	$0.0692$	$0.0885$
${\Theta }_R$		$0.0655$	$0.0714$	$0.2878$		$0.1354$	$0.2470$	$0.5242$	\rdelim}40mm[Table 2]
${\Theta }_{\theta }$		$0.3539$	$0.0312$	$0.0251$		$0.7312$	$0.1079$	$0.0457$
${\Theta }_d$		$0.0225$	$0.1646$	$0.0823$		$0.0465$	$0.5697$	$0.1499$
${\Theta }_l$		$0.0420$	$0.0218$	$0.1539$		$0.0868$	$0.0754$	$0.2802$

Table 5. Table 5 . Two-way nonparametric Cramér test ( [ 4 ] ) rejects the hypothesis of equality of the empirical marginal and joint distributions for comparable variables in Figures   11(c) , 11(d) , 7(c) and  8(c) with high statistical significance (each test is based on 2000 2000 2000 values sampled using a gPCE).

Variable	Cramér-statistic	Critical value	Conf. interval	$p$-value	Result	Figures
$D_L$	$4.688$	$0.3389$	$0.95$		reject	\rdelim}30mm[11(c) vs. 7(c)]
$D_T$	$79.08$	$0.7396$	$0.95$		reject
${\gamma }_{\mathrm{eff}}$	$0.4619$	$0.3958$	$0.95$	$0.029$	reject
$(D_L,D_T)$	$72.31$	$0.7827$	$0.95$		reject	\rdelim}30mm[11(d) vs. 8(c)]
$(D_T,{\gamma }_{\mathrm{eff}})$	$72.63$	$0.8284$	$0.95$		reject
$({\gamma }_{\mathrm{eff}},D_L)$	$4.734$	$0.4882$	$0.95$		reject

Equations

$$\frac{\partial c}{\partial t}=\nabla \cdot (D\nabla c),\mathbf{x}\in \mathcal{P},t>0,$$

$$c(\mathbf{x},0)=c_{\mathrm{in}},\mathbf{x}\in \mathcal{P},$$

$$-D\mathbf{n}\cdot \nabla c=q_{\mathrm{m}}\frac{\partial s}{\partial t},\mathbf{x}\in \Gamma ,t>0,$$

$$\frac{\mathrm{d}s}{\mathrm{d}t}=\gamma (s_{\mathrm{eq}}-s),$$

$$s_{\mathrm{eq}}=\frac{Kc}{1+Kc},$$

$$u(\mathbf{x},t):=\frac{1}{\Vert \breve{\mathcal{V}}\Vert }{\int }_{\breve{\mathcal{V}}(\mathbf{x})}c(\mathbf{\xi },t)\mathrm{d}\mathbf{\xi }=\frac{1}{\Vert \breve{\mathcal{V}}\Vert }{\int }_{\breve{\mathcal{P}}(\mathbf{x})}c(\mathbf{\xi },t)\mathrm{d}\mathbf{\xi }=\frac{\varphi }{\Vert \breve{\mathcal{P}}\Vert }{\int }_{\breve{\mathcal{P}}(\mathbf{x})}c(\mathbf{\xi },t)\mathrm{d}\mathbf{\xi },\mathbf{x}\in \mathcal{V},$$

$$\varphi \frac{\partial u}{\partial t}=\nabla \cdot ({\mathbf{D}}_{\mathrm{eff}}\nabla u)-\varphi q_{\mathrm{m}}{\gamma }_{\mathrm{eff}}\frac{Ku}{1+Ku},\mathbf{x}\in \mathcal{V}.$$

$${\gamma }_{\mathrm{eff}}=\frac{\Vert \breve{\Gamma }\Vert }{\Vert \breve{\mathcal{P}}\Vert },$$

$${\mathbf{D}}_{\mathrm{eff}}=\frac{1}{\Vert \breve{\mathcal{V}}\Vert }{\int }_{\breve{\mathcal{P}}}(\mathbf{I}+{\nabla }_{\mathbf{\xi }}{\mathbf{\chi }}^{\top })\mathrm{d}\mathbf{\xi },$$

$${\nabla }_{\mathbf{\xi }}\cdot (D{\nabla }_{\mathbf{\xi }}\mathbf{\chi })=0,\mathbf{\xi }\in \breve{\mathcal{P}},$$

$$⟨\mathbf{\chi }⟩:=\frac{1}{\Vert \breve{\mathcal{V}}\Vert }{\int }_{\breve{\mathcal{P}}}\mathbf{\chi }(\mathbf{\xi })\mathrm{d}\mathbf{\xi }=0,$$

$$\mathbf{n}\cdot {\nabla }_{\mathbf{\xi }}\mathbf{\chi }=-\mathbf{n}\cdot \mathbf{I},\mathbf{\xi }\in \breve{\Gamma },$$

$${\chi }_1(-a,{\xi }_2)={\chi }_1(a,{\xi }_2)=0,\frac{\partial {\chi }_1}{\partial {\xi }_2}({\xi }_1,0)=\frac{\partial {\chi }_1}{\partial {\xi }_2}({\xi }_1,b)=0,$$

$${\chi }_2({\xi }_1,0)={\chi }_2({\xi }_1,b)=0,\frac{\partial {\chi }_2}{\partial {\xi }_1}(-a,{\xi }_2)=\frac{\partial {\chi }_2}{\partial {\xi }_1}(a,{\xi }_2)=0,$$

$$a=R\cos \theta \text{and}b=2R\cos ({\sin }^{-1}(\frac{d}{2R}))+l=2R\sqrt{1-\frac{d^2}{4R^2}}+l$$

$$P(\mathbf{X}\mid \mathbf{\Theta })={\delta }_{\mathbf{X}}(\mathbf{X}-\mathbf{\chi }(\mathbf{\xi };\mathbf{\Theta })).$$

$$P(U\mid \mathbf{X})={\delta }_U(U-u(\mathbf{x},t;\mathbf{X})).$$

$$\begin{array}{r}P:=P(\mathbf{\Theta },\mathbf{X},U)=P(U\mid \mathbf{X})P(\mathbf{X}\mid \mathbf{\Theta })P(\mathbf{\Theta })={\delta }_U(U-u(\mathbf{x},t;\mathbf{X})){\delta }_{\mathbf{X}}(\mathbf{X}-\mathbf{\chi }(\mathbf{\xi },t;\mathbf{\Theta }))P(\mathbf{\Theta }).\end{array}$$

$$P({\Theta }_i,{\Theta }_j)=P({\Theta }_i)P({\Theta }_j),$$

$$P(\mathbf{\Theta })=\prod _{i=1}^{n}P({\Theta }_i).$$

$$P_0:={\delta }_U(U-u(\mathbf{x},t;\mathbf{X})){\delta }_{\mathbf{X}}(\mathbf{X}-\mathbf{\chi }(\mathbf{\xi },t;{\Theta }_1,\dots ,{\Theta }_n))\prod _{i=1}^{n}P({\Theta }_i).$$

$${\Theta }_R\sim \mathrm{unif}(R_{-},R_{+})\text{and}{\Theta }_{\theta }\sim \mathrm{unif}({\theta }_{-},{\theta }_{+}),$$

$${\Theta }_d\mid {\Theta }_R,{\Theta }_{\theta }\sim \mathrm{unif}(d_{-},\min \{2{\Theta }_R\cos {\Theta }_{\theta },d_{+}\}).$$

$${\Theta }_l\mid {\Theta }_R,{\Theta }_d\sim \mathrm{unif}(\max \{l_{-},2{\Theta }_R-\sqrt{4{\Theta }_R^2-{\Theta }_d^2}\},l_{+}),$$

$$P(\mathbf{\Theta })=P({\Theta }_R,{\Theta }_{\theta },{\Theta }_l,{\Theta }_d)=P({\Theta }_l\mid {\Theta }_R,{\Theta }_d)P({\Theta }_d\mid {\Theta }_R,{\Theta }_{\theta })P({\Theta }_R)P({\Theta }_{\theta })$$

$$P_1:={\delta }_U(U-u(\mathbf{x},t;\mathbf{X})){\delta }_{\mathbf{X}}(\mathbf{X}-\mathbf{\chi }(\mathbf{\xi },t;\mathbf{\Theta }))P({\Theta }_l\mid {\Theta }_R,{\Theta }_{\theta },{\Theta }_d)P({\Theta }_d\mid {\Theta }_R,{\Theta }_{\theta })P({\Theta }_R)P({\Theta }_{\theta }),$$

$$D_L:={\mathbf{D}}_{\mathrm{eff}}^{11}\text{and}{D}_T:={\mathbf{D}}_{\mathrm{eff}}^{22},$$

$$P(g\mid \mathbf{\Theta })=P(g\mid \mathbf{X})P(\mathbf{X}\mid \mathbf{\Theta })P(\mathbf{\Theta }).$$

$$P(D_L\mid \mathbf{\Theta })=P(D_L\mid \mathbf{X}){\delta }_{\mathbf{X}}(\mathbf{X}-\mathbf{\chi }(\mathbf{\xi },t;\mathbf{\Theta }))P({\Theta }_l\mid {\Theta }_R,{\Theta }_{\theta },{\Theta }_d)P({\Theta }_d\mid {\Theta }_R,{\Theta }_{\theta })P({\Theta }_R)P({\Theta }_{\theta }).$$

$$g(\mathbf{\Theta })=\sum _{i=0}^{\infty }{G}_i{\Psi }_i(\mathbf{\Theta })\approx \sum _{i=0}^{N_{PC}}{G}_i{\Psi }_i(\mathbf{\Theta })=:\hat{g},$$

Full text

∗Both authors contributed equally to this work.

Key words. Bayesian networks, causality, multiscale modeling, porous media, energy storage, uncertainty quantification, global sensitivity analysis, mutual information

1. Introduction

Understanding statistical and causal relations between properties/model parameters at various scales is essential for science-based predictions in general and for forecasts of transport phenomena in porous media in particular. For example, the design of materials for energy storage devices aims to optimize macroscopic material properties (quantities of interest or QoIs), such as effective diffusion coefficient and capacitance, through engineered pore structures [52, 51]. Quantification of both uncertainty in predictions of these macroscopic QoIs and their sensitivity to variability and uncertainty in microscopic features are crucial for informing such decision tasks as optimal experimental design and reliability engineering.

The simulation-assisted approach to the optimal design of porous meta-materials takes advantage of the availability of microscopic (pore-scale) and macroscopic (continuum or Darcy-scale) models, as well as of a bridge between the two provided by various upscaling techniques [52, 51, 24]. Such a bridge also facilitates analysis of both uncertainty propagation from the microscopic scale to the macroscopic scale and sensitivity of macroscopic material properties to the microscopic ones [47]. Predictions resulting from this multiscale approach can be made more robust by incorporating information about correlations and causal relationships between scales and within a single scale. Causality stems for example from physical, chemical, and/or engineering design constraints. Our primary objective is to bring a Bayesian network perspective to the incorporation of causality into modeling process.

Bayesian networks are a special class of hierarchical probabilistic graphical models (PGMs) with a directed acyclic graph structure that represent causal relationships among random variables [37, 38, 19]. Bayesian networks and, more generally, PGMs provide a rich framework for encoding distributions over large, complex domains of interacting random variables that can include causal relationships and expert knowledge.

In our application, they provide a coherent framework for representing causal relationships both between problem scales and among the space of parameters representing pore-scale features. Incorporating Bayesian networks into random PDEs is a novel approach to the modeling of multiscale porous media; it breaks down the stochastic modeling and statistical inference tasks into smaller, controllable parts enabling us to

1. (i)

build systematically informed parameter priors that include physical constraints and/or correlations; 2. (ii)

mirror engineering processes related to the design of hierarchical nanoporous media networks; 3. (iii)

construct Bayesian network (random) Darcy-scale PDE models informed by (possibly uncertain) pore-scale data, parameters, and constraints; and 4. (iv)

carry out global sensitivity analysis (GSA) and uncertainty quantification (UQ).

This framework for incorporating causal relationships through structured priors demands GSA tools that differ from the standard variance-based GSA methods. These typically assume unstructured (i.e., mutually independent) priors and are neither easy to interpret nor cheap to compute for correlated inputs [30, 18]. In our application, causal relationships exist not just between parameters but also between scales. The latter is important since our predictive PDFs are not necessarily Gaussian and/or do not have a known analytical form. This suggests moving away from a fixed number of moments to a moment-independent quantity such as mutual information. We employ a moment-independent GSA relying on mutual information [12, 45] and empirical distributions acquired through simulations. We demonstrate that differential mutual information provides a measure of input effects that is suited to tackling the twin challenges of structured or correlated inputs and non-Gaussian QoIs.

Design of (nano)porous metamaterials for energy storage provides an ideal setting to illustrate the Bayesian network PDE approach. Macroscopic material properties are dependent on a set of microscopic parameters characterizing the pore geometry, e.g., pore radius or pore connectivity. These microscopic parameters are typically correlated due to the presence of geometric and topological constraints and uncertain due to natural and/or manufacturing variability. This setting gives rise to a number of theoretical and practical questions: How does uncertainty in microscopic properties (quantified, e.g., in terms of a pore-size distribution) propagate to the macroscopic scale (expressed in terms of the PDF of, e.g., the effective diffusion coefficient)? How sensitive are a material’s macroscopic properties to its microscopic counterparts? Etc. Our computational strategy, which makes exhaustive sampling for prediction and uncertainty quantification feasible, has three ingredients: i) Rosenblatt transforms to decorrelate inputs for non-intrusive scientific computing, ii) generalized polynomial chaos expansions obtained using the DAKOTA software [1], and iii) kernel density estimation techniques.

In Section 2, we formulate a macroscopic (Darcy-scale) model of reactive transport in hierarchical nanoporous media; the model parameters are expressed in terms of microscopic (pore-scale) material properties by means of homogenization [47], which facilitates multiscale UQ and GSA. Section 3 contains a description of our Bayesian network-based approach to linking the components of this model across the two scales. Section 4 contains details of its implementation highlighting the use of the inverse Rosenblatt transform to non-intrusive utilization of existing software, such as DAKOTA. This section also collates results of our numerical experiments, which demonstrate the importance of causality. In Section 5, we adopt sensitivity indices based on differential mutual information and provide ranking for the impact of uncertainty in the microscopic parameters on uncertainty in their macroscopic counterparts. Our examples illustrate how the inclusion of causal relationships encoding structural constraints provides rankings more consistent with the physics anticipated for a simple hierarchical nanoporous material. In Section 6, we present an alternative Bayesian network to highlight the method’s flexibility and ability to mirror distinct engineering and design processes. Major conclusions drawn from our study are summarized in Section 7.

2. Models of flow and transport in nanoporous materials

A volume $\mathcal{V}=\mathcal{P}\cup\mathcal{S}$ of a hierarchical porous material is comprised of a fluid-filled pore space $\mathcal{P}$ and solid matrix $\mathcal{S}$, with a (multi-connected) fluid-solid interface denoted by $\Gamma=\mathcal{P}\cap\mathcal{S}$. To mimic a manufacturing process and to make subsequent use of the homogenization theory, we assume that the volume $\mathcal{V}$ consists of a periodic arrangement of unit cells ${\breve{\mathcal{V}}}={\breve{\mathcal{P}}}\cap{\breve{\mathcal{S}}}$ with pore space ${\breve{\mathcal{P}}}\subset\mathcal{P}$, solid matrix ${\breve{\mathcal{S}}}\subset\mathcal{S}$, and fluid-solid interface ${\breve{\Gamma}}={\breve{\mathcal{P}}}\cap{\breve{\mathcal{S}}}$. For example, the hierarchical nanoporous material in Figure 1 consists of mesopores that are connected longitudinally (horizontally) through nanotunnels and transversely (vertically) by a series of nanotubes. These features are described by a set of parameters $\{R,\theta,d,l\}$, where $R$ is the mesopore radius; $\theta$ is the angle of overlap between adjacent mesopores in a nanotunnel; and $d$ and $l$ are the diameter and length of the nanotubes, which serve as nano-bridges between adjacent mesopores/nanotunnels.

The design of novel materials calls for a systematic analysis of the sensitivity of desired macroscopic (Darcy-scale) properties to imperfections (natural variability) in microscopic (pore-scale) parameters and/or their distributions. Following [47], we use the homogenization theory to map uncertainty in the microscopic parameters and processes to their macroscopic counterparts. Sources of uncertainty, as well as representations of randomness in the microscopic and macroscopic models are described below.

2.1. Pore-scale model

At the pore-scale, the evolution of a solute concentration $c(\boldsymbol{x},t)$ ($\mathrm{mol}\text{/}{\mathrm{\ell}}^{3}$), at point $\boldsymbol{x}\in\mathcal{P}$ and time $t>0$, is governed by the evolution equation

[TABLE]

where $D(\boldsymbol{x})$ (${\mathrm{L}}^{2}\text{/}\mathrm{T}$), $\boldsymbol{x}\in\mathcal{P}$, is the pore-scale diffusion coefficient. The spatial variability of $D$ allows for Fickian diffusion through mesopores and Knudson diffusion through nanotubes. This equation is subject to the uniform initial condition

[TABLE]

and the boundary condition

[TABLE]

where $q_{\mathrm{m}}$ and $s(\boldsymbol{x},t)$ are related to the sorption properties of the material surface $\Gamma$. Specifically, $q(\boldsymbol{x},t)=q_{\mathrm{m}}\cdot s(\boldsymbol{x},t)$ is the adsorption amount per unit area of $\Gamma$ ($\mathrm{mol}\text{/}{\mathrm{\ell}}^{2}$), $q_{\mathrm{m}}$ ($\mathrm{mol}\text{/}{\mathrm{\ell}}^{2}$) is the maximal adsorption amount, and $s(\boldsymbol{x},t)$ is the fractional coverage of $\Gamma$. The fractional coverage is assumed to follow Lagergren’s pseudo-first-order rate equation,

[TABLE]

where $\gamma$ ($\text{/}\mathrm{T}$) is the adsorption rate constant and the equilibrium adsorption coverage fraction $s_{\mathrm{eq}}$ satisfied Langmuir’s adsorption isotherm,

[TABLE]

with the adsorption equilibrium constant $K$ (${\mathrm{\ell}}^{3}\text{/}\mathrm{mol}$).

2.2. Darcy-scale model

At the macroscopic scale, the volume-averaged solute concentration $u$,

[TABLE]

treats a porous material as a continuum. Here, $\|\cdot\|$ denotes the volume of a domain and $\phi\mathrel{\mathop{:}}=\|\mathcal{P}\|/\|\mathcal{V}\|=\|{\breve{\mathcal{P}}}\|/\|{\breve{\mathcal{V}}}\|$ is the material porosity. Using homogenization via multiple-scale expansions [52], one can show that $u$ satisfies a reaction-diffusion equation

[TABLE]

The effective diffusion coefficient $\mathbf{D}_{\mathrm{eff}}$ and the effective rate constant $\gamma_{\mathrm{eff}}$ are random, stemming from uncertainty in pore-scale structures and processes that is propagated by the homogenization map. Specifically, the effective rate constant $\gamma_{\mathrm{eff}}$ ($\text{/}\mathrm{L}$) is computed as

[TABLE]

i.e., is defined solely by the pore geometry; and the effective diffusion coefficient $\mathbf{D}_{\mathrm{eff}}$ (L2/T), a second rank tensor, depends on both the pore geometry and the pore-scale processes. It is computed in terms of a closure variable $\boldsymbol{\chi}$ as

[TABLE]

where $\boldsymbol{I}$ is the identity matrix. The closure variable $\boldsymbol{\chi}(\boldsymbol{\xi})$ is a ${\breve{\mathcal{V}}}$-periodic vector defined on ${\breve{\mathcal{P}}}$, which satisfies the Laplace equation

[TABLE]

subject to the normalizing condition

[TABLE]

the boundary condition along the fluid-solid segments ${\breve{\Gamma}}$,

[TABLE]

and ${\breve{\mathcal{V}}}$-periodic boundary condition on the remaining fluid segments of the boundary of ${\breve{\mathcal{P}}}$. For the hierarchical nanoporous material in Figure 1, these auxiliary conditions reduce to

[TABLE]

and

[TABLE]

where

[TABLE]

are, respectively, the width and height of the unit cell ${\breve{\mathcal{V}}}$. In the next section, we introduce a representation for uncertainty in the pore- and Darcy-scale models that will enable us to compartmentalize various stochastic modeling and statistical inference tasks.

3. Bayesian network formulation for random PDE models of

multiscale porous media

A Bayesian network is a directed graph structure, in which each node represents a random variable with an associated PDF and each edge encodes a causal relationship [19]. A Bayesian network PDE incorporates these structured probabilistic models into a forward physical model, and allows one to capture causal relationships in a systematic way. The key components of the full statistical model include:

1. (i)

inputs $\boldsymbol{\Theta}=\{R,\theta,d,l\}$, a random vector related to the parameters describing pore-scale features in Figure 1; 2. (ii)

upscaling variable $\boldsymbol{X}$, a random vector related to the closure equations (6)–(8); and 3. (iii)

QoI $U$, a random macroscopic quantity, such as macroscopic parameters $\mathbf{D}_{\mathrm{eff}}$ and $\gamma_{\mathrm{eff}}$ in 4 and 5 or a functional of $u$ in 3.

Causal relationships arise naturally among the problem scales and hence these components: the PDF of $U$ depends on the PDF of $\boldsymbol{X}$ since, for example, $u$ and $\mathbf{D}_{\mathrm{eff}}$ depend on the closure variable $\boldsymbol{\chi}$. In turn, the PDF of $\boldsymbol{X}$ depends on the PDF of $\boldsymbol{\Theta}$ since $\boldsymbol{\chi}$ depends on the pore-scale parameter values, e.g., via 9 and 10. The Bayesian network in Figure 2 describes these causal relationships for the full statistical model. The directed network structure indicates that there is only one-way communication between the different scales.

More precisely, the Bayesian network in Figure 2 encodes conditional relationships among the PDFs for $\boldsymbol{\Theta}$, $\boldsymbol{X}$, and $U$. We denote the joint PDF of the input parameters by $P(\boldsymbol{\Theta})$. For simplicity we assume the statistical models for $\boldsymbol{X}$ and $U$ to be known, that is, uncertainty enters only through the parameters since we have fixed both the form of equations for $\boldsymbol{\chi}$ and $u$ in (3)–(8) and the numerical methods for approximating them. Under this assumption the conditional PDF of $\boldsymbol{X}$ given a sample $\boldsymbol{\Theta}$ is the Dirac delta function

[TABLE]

Similarly, the conditional PDF of $U$ given a sample $\boldsymbol{X}=\boldsymbol{\chi}$ is

[TABLE]

Then the full statistical model $P$ is given by

[TABLE]

In the remainder of this section, we discuss expanded Bayesian networks for representing causal relationships between parameters that allow us to encode correlations among pore scale features. In Section 3.1, we consider the assumption of independent priors for the pore-scale parameters and contrast this in section Section 3.2 with causality arising from natural structural constraints encountered in engineering design.

Remark 1:

Widely adopted approaches to uncertainty quantification, e.g., Monte Carlo methods, involve strategies for generating surrogates of the model $P$ in 13. The direct application of these traditional UQ methods require the form of $P$ to be known. In the sequel we introduce an approach relying on Rosenblatt transformations, generalized polynomial chaos expansions, and the popular UQ software package DAKOTA that makes sampling $P$ feasible.

3.1. Independent uniform priors

As a first simple case we revisit the analysis [47] of the hierarchical nanopore geometry in Figure 1, but in the context of the Bayesian network perspective. A naive model for representing uncertainty in each of the pore-scale parameters assumes pairwise independent priors. Recall that for $\boldsymbol{\Theta}=(\Theta_{1},\dots,\Theta_{n})$, the variables $\Theta_{1},\dots,\Theta_{n}$ are pairwise independent, $\Theta_{i}\mathrel{\perp}\Theta_{j}$ for all $i,j=1,\dots,n$ such that $i\neq j$, if and only if

[TABLE]

where $P(\Theta_{i})$ denotes the (marginal) PDF of $\Theta_{i}$ and $P(\Theta_{i},\Theta_{j})$ denotes the joint PDF of $(\Theta_{i},\Theta_{j})$ ([19]). Under the assumption of independent priors $P(\Theta_{i})$, the joint PDF $P(\boldsymbol{\Theta})$ factors into the product of the priors,

[TABLE]

The full statistical model for $P$ in 13 is then given by

[TABLE]

In the case of 15, the Bayesian network has the special form in Figure 3 where the independent priors assumption leads to an overall flatness in the graph structure for $P(\boldsymbol{\Theta})$. For PDF $P(\Theta_{i})$ with finite variance, the pairwise independence assumption implies the variables $\Theta_{i}$ are uncorrelated. In particular, 15 holds with $|\boldsymbol{\Theta}|=n=4$ for the parameters $\{R,\theta,l,d\}$ that describe pore-scale features in Figure 1. A model or PDF can then be specified for each $\Theta_{i}$, for example, the uniform priors considered in [47].

The Bayesian network in Figure 3 does not take into account geometrical structural constraints between the pore scale parameters that naturally arise when considering a periodic arrangement of unit cells ${\breve{\mathcal{V}}}$ as in the hierarchical nanoporous structure in Figure 1. For example, assuming independent priors $P(\Theta_{i})$ each uniformly distributed according to the hyperparameter ranges in Tables 1 and 2, it is possible to sample geometries that are inconsistent with Figure 1. Next we build Bayesian networks based on causal relationships that encode such natural structural constraints.

3.2. Correlations arising from pore-scale structural

constraints

Inclusion of causal relationships imposed by geometrical constraints on the pore-scale parameters adds more complexity to the structure of $\boldsymbol{\Theta}$ in the Bayesian network in Figure 3. Here and below we use $\Theta_{p}$ with labels $p\in\{R,\theta,l,d\}$ in place of indices where no confusion arises. Moreover, for each $\Theta_{p}$ we fix hyperparameters $\{p_{+},p_{-}\}$ corresponding to upper and lower bounds on $\Theta_{p}$ in Tables 1 and 2; performing inference over the Bayesian network to infer hyperparameters from relevant data is beyond the scope of this work. In particular, there are a plurality of Bayesian networks that describe structural constraints consistent with Figure 1. Different causal relationships mirror distinct engineering or design processes (see Section 6) and yield different correlation structures among pore-scale features (cf. Figure 6 and Figure 11(b)).

Figure 4 presents one possible Bayesian network capturing causal relationships that encode structural constraints among the hierarchical pore-scale parameters $\{R,\theta,l,d\}$. Choosing $\{\Theta_{R},\Theta_{\theta}\}$ as independent parameters to be consistent with our design goal and assuming uniform priors,

[TABLE]

constrains both the distribution of $\Theta_{d}$ and of $\Theta_{l}$. Specifically, in order for the sample $\Theta_{d}=d$ to be consistent with the features of the hierarchical nanopore structure in Figure 1, the nanotube diameter cannot exceed the width of the unit cell ${\breve{\mathcal{V}}}$, i.e., $d<2R\cos\theta$ (see Figure 5). That is, we are naturally able to specify the conditional PDF $\Theta_{d}$ given samples $\Theta_{R}=R$ and $\Theta_{\theta}=\theta$. Assuming an uninformed or uniform model for this conditional PDF we have

[TABLE]

The length of the nanotube is bounded below by $l>2R-\sqrt{4R^{2}-d^{2}}$ when the vertical gap between the mesopores is zero (see Figure 5) and then the PDF of $\Theta_{l}$ given $\Theta_{R}$ and $\Theta_{d}$ is

[TABLE]

where again we assume a uniform model for the conditional PDF of $\Theta_{l}$.

The joint distribution for the pore-scale input parameters is then given by

[TABLE]

where each of the conditional and prior PDFs is specified in (16)–(18) (cf. to the assumption of independent priors in 14). Once a range of hyperparameters is fixed, one can sample from the PDFs (16)–(18) and hence the correlation structure of $P(\boldsymbol{\Theta})$ can be computed empirically. In Figure 6, we compare the empirical correlation structure obtained from 19 over both a narrow hyperparameter range in Table 1 and a physical hyperparameter range in Table 2 to the correlation structure obtained from 14. The model with independent priors is not valid over the physical range in Table 2.

Incorporating 19 into the full statistical model 13, we obtain a new probabilistic model $P_{1}$,

[TABLE]

that describes the full statistical model for the multiscale system with the causal relationships in Figure 4 assuming that uncertainty only enters through the parameters. When the models for $U$ and $\boldsymbol{X}$ are known and trivial (e.g., have a known analytic form), then sampling 20 is a straightforward task. In the next section, we review tools that will enable us to feasibly sample the statistical models for $U$ and $\boldsymbol{X}$ in order to carry out UQ and GSA.

3.3. Constructing a physics-informed probabilistic model for

macroscopic QoIs

We are interested in functionals $g=g(\boldsymbol{\Theta})=g(\mathbf{D}_{\mathrm{eff}}(\boldsymbol{\Theta}),\gamma_{\mathrm{eff}}(\boldsymbol{\Theta}))$ such as the projections,

[TABLE]

corresponding to the longitudinal and the transverse components of the effective diffusion coefficient tensor. In the hierarchical nanoporous media in Figure 1, ${D_{L}}$ is related to diffusion through nanotunnels/mesopores and ${D_{T}}$ through nanotubes. In the numerical experiments presented in the sequel, we report uncertainty in these macroscopic quantities by giving an estimate of the PDF $f_{g}\sim P(g\mid\boldsymbol{\Theta})$ for the PDF of $g$ given knowledge of pore-scale inputs. In the context of the homogenized pore- to Darcy-scale model in Figure 2, we have the representation

[TABLE]

For example, the univariate PDF for ${D_{L}}=g(\mathbf{D}_{\mathrm{eff}}(\boldsymbol{\Theta}),\gamma_{\mathrm{eff}}(\boldsymbol{\Theta}))$ based on the Bayesian network in Figure 4 is

[TABLE]

In general, the form of the statistical model or PDF for $P(g\mid\boldsymbol{X})$ and hence $P(g\mid\boldsymbol{\Theta})$ in 22 is unknown. However, the PDF in 22 can be estimated empirically using simulations of the forward model. Sampling 22 based strictly on input-output pairs may be computationally expensive due to the high cost involved in simulating the forward model.

The computation of estimates for QoIs of the form 22 is made feasible using a two-step method that relies on first finding a truncated generalized polynomial chaos expansion (gPCE) for the variable $g$ and second using this surrogate to build an appropriate kernel density estimator (KDE) for the desired distribution. A surrogate $\hat{g}$ for $g(\boldsymbol{\Theta})$ is given by the gPCE ([22, 50, 15]),

[TABLE]

where $\Psi_{i}(\boldsymbol{\Theta})$ are an orthogonal multivariate polynomial basis, $G_{i}$ are the expansion coefficients, and the expansion is truncated after $N_{PC}$ terms such that

[TABLE]

where $\kappa_{i}$ is the polynomial order bound for the $i^{\mathrm{th}}$ dimension and $N_{p}$ is the number of parameters.

The second step involves producing $N$ samples $\hat{g}^{k}$ using the gPCE 24 corresponding to realizations $\boldsymbol{\Theta}^{k}$ for $k=1,\dots,N$, for a fixed number $N$. These samples are then used to construct a KDE $\bar{f}_{\hat{g}}$ for the desired density $f_{g}$, for example, using a Gaussian-kernel,

[TABLE]

where $h$ is the kernel bandwidth; for the multivariate density $\bar{f}_{\hat{g}_{1},\hat{g}_{2}}$ we similarly employ Gaussian-kernels

[TABLE]

with bandwidths $h_{1}$ and $h_{2}$. Software DAKOTA [1] was used to automatize the process of computing the coefficients, basis functions, and truncation parameters appearing in 24. This approach is taken in [47] in the context of independent priors with the aim of computing Sobol’ indices for global sensitivity analysis. Although generalizations of gPCE that handle correlated inputs exist (e.g., [31, 36]), we instead present a recipe for obtaining the desired gPCE in 24 that non-intrusively utilize existing codes and software packages such as DAKOTA that implement methods assuming uncorrelated inputs in the next section.

4. Quantifying uncertainty in Darcy-scale flow variables

The present section deals with uncertainty quantification for the Bayesian network PDE model for multiscale porous media outlined in Sections 2 and 3. The causal relationships encoded by the Bayesian network for the full statistical model in Figure 2 propagate uncertainty from pore-scale parameters to Darcy-scale variables via the homogenization map in 4, 5, 6, 7, 8, 9 and 10. Together, these allow one to study systematically the impact of microscopic structural uncertainty on macroscopic flow variables. At present, we incorporate the Bayesian networks constructed in Section 3 for informed priors into the random PDE homogenization framework thereby examining through simulations and numerical experiments the role of causality in predictions of Darcy-scale flow variable QoIs. Specifically, we report on numerical experiments concerning the joint and marginal distributions of Darcy-scale flow variables where uncertainty stems from pore-scale features with correlations arising from structural constraints encoded by the Bayesian network in Figure 4. As these numerical experiments utilize DAKOTA ([1]) to compute gPCE surrogates for Darcy-scale variables, we first illustrate in Section 4.1 a technique for decorrelating inputs to allow the non-intrusive use of existing codes and software packages. This work-flow, which relies on Rosenblatt transformations, gPCEs, and DAKOTA, enables uncertainty quantification and global sensitivity analysis by making it feasible to sample from the desired QoI with respect to a given statistical model $P$ in 13.

4.1. Non-intrusive input decorrelation using Rosenblatt

transforms

Many variance-based methods for uncertainty quantification and sensitivity analysis, such as Sobol’ indices, and hence the popular software packages that implement these methods, require models that assume statistically independent inputs. Bayesian networks, recall Figures 3 and 4, encode correlations through the specification of causal relationships (see Figure 6). Presently, we highlight how the Rosenblatt transform can be used to decorrelate inputs by mapping a vector of random variables with a specified joint distribution onto a vector of independent uniform random variables when the conditional distributions are known. This procedure, in Algorithm 1 below, enables the non-intrusive use of DAKOTA and existing codes for solving the forward model for our application of interest when the conditional dependencies are represented using Bayesian networks.

The Rosenblatt transform ([40]) turns the problem of sampling a general joint distribution into the problem of sampling a vector of independent $\mathrm{unif}(0,1)$ random variables. Let $\boldsymbol{X}=(X_{1},\dots,X_{k})$ be a random vector with a continuous joint cumulative distribution function $F(x_{1},\dots,x_{k})$. Define a transform, $\mathcal{T}(x)=\mathcal{T}(x_{1},\dots x_{k})=(z_{1},\dots,z_{k})=z$, given by

[TABLE]

where $F_{i\mid j}$ is the conditional cumulative distribution function of $X_{i}$ given $X_{j}$, i.e. $F_{i\mid j}(x_{i}|x_{j})=P(X_{i}<x_{i}\mid X_{j}=x_{j})$. The Rosenblatt transform, $\boldsymbol{Z}\mathrel{\mathop{:}}=\mathcal{T}(\boldsymbol{X})$, yields $\boldsymbol{Z}=(Z_{1},\dots,Z_{k})$ uniformly distributed on the $k$-dimensional hypercube, that is, $Z_{1},\dots,Z_{k}$ are independent and identically distributed (iid) $\mathrm{unif}(0,1)$ random variables. Note that this transform depends on the ordering of the elements in the vector $\boldsymbol{X}$ due to the serial nature of the conditioning; we denote the Rosenblatt transform and the inverse, when it exists, associated with the ordering of a particular vector $\boldsymbol{X}$ with a subscript, e.g. $\mathcal{T}_{\boldsymbol{X}}$.

For our application of interest, a target vector $\boldsymbol{\Theta}$ with a causal structure encoded by a Bayesian network can be obtained by applying the inverse Rosenblatt transform to a vector of independent uniform variables. For example, given the random vector of parameters $\boldsymbol{\Theta}=(\Theta_{R},\Theta_{\theta},\Theta_{l},\Theta_{d})$ with joint distribution in 19, the transform 28 simplifies to

[TABLE]

due to the independence and conditional independence of the variables ([19]), as indicated in Figure 4 by the absence of a causal relationships between $\Theta_{R}$ and $\Theta_{\theta}$ and between $\Theta_{l}$ and $\Theta_{d}$, respectively. Thus, the Rosenblatt transform is given by

[TABLE]

where, using the statistical models indicated in 16, 17 and 18, the components of $\boldsymbol{Z}$ are

[TABLE]

The corresponding inverse Rosenblatt transform $\boldsymbol{\Theta}=\mathcal{T}_{\boldsymbol{\Theta}}^{-1}(\boldsymbol{Z})$ is, component-wise,

[TABLE]

where the target vector $\boldsymbol{\Theta}$ has the distribution $P(\boldsymbol{\Theta})$ given in 19. Thus, the inverse transform maps a random vector $\boldsymbol{Z}\sim\mathrm{unif}(0,1)^{k}$ into a target distribution $P(\boldsymbol{\Theta})$ using knowledge of the conditional dependencies associated with the Bayesian network for $P(\boldsymbol{\Theta})$, in particular, where the conditional distribution $\Theta_{j}$ given $\Theta_{j-1},\dots,\Theta_{1}$ is known for each $j=1,\dots,k$.

Together, Bayesian networks and Rosenblatt transforms provide a strategy for non-intrusively incorporating constraints and correlations into existing computational frameworks. Algorithm 1 describes the use of DAKOTA for computing surrogates for Darcy-scale QoIs based on correlated inputs given an inverse Rosenblatt transform and an existing solver for the forward problem. A surrogate $\hat{g}$ for a QoI $g(Y(\boldsymbol{\Theta}))$, e.g. the gPCE coefficients $G$ in 24 and truncation parameter $N_{PC}$ in 25, can be computed with DAKOTA using several forward simulations of the response $Y=\mathcal{M}(\boldsymbol{\Theta})$ where $\mathcal{M}$ denotes the portion of the forward model solver that maps a random input $\boldsymbol{\Theta}$ to $Y$. For example, if the QoI is $g(Y(\boldsymbol{\Theta}))={D_{L}}$ then $\mathcal{M}$ would correspond to the projection $\mathbf{D}_{\mathrm{eff}}^{11}$ in 21 of the solution to the coupled system 5 with 6, 7 and 8. The compositional model $Y=(\mathcal{M}\circ\mathcal{T}^{-1})(\boldsymbol{Z})$ that first employs the inverse Rosenblatt transform provides a non-intrusive means of computing with DAKOTA since the statistics of the output of $\mathcal{M}\circ\mathcal{T}^{-1}$ in response to $\boldsymbol{Z}$ are identical to the statistics of the output of $\mathcal{M}$ in response to $\boldsymbol{\Theta}$ ([46]). An important observation is that Algorithm 1 returns a surrogate $\hat{g}$ for $g(Y(\boldsymbol{\Theta}))$ with respect to the input variables $\boldsymbol{Z}$, not for $\boldsymbol{\Theta}$ directly, and therefore care must be taken if reporting Sobol’ indices with respect to the input parameters ([30]). In the next section, we observe that the density functions for effective Darcy-scale QoIs exhibit non-Gaussian behavior. Together, these two challenges—correlated inputs and non-Gaussian QoIs—motivate the investigation of moment-independent global sensitivity indices in Section 5.

4.2. Numerical experiments: incorporating causality in

predictions of Darcy-scale flow

The numerical experiments that follow employ the following common setup that makes sampling the distribution of QoIs computationally feasible for uncertainty quantification and global sensitivity analysis. DAKOTA is used as a wrapper to map random inputs on pore-scale parameters to Darcy-scale responses yielding gPCE surrogates for Darcy-scale variables, i.e. effective longitudinal diffusion ${D_{L}}$, effective transverse diffusion ${D_{T}}$, and effective sorption rate constant $\gamma_{\mathrm{eff}}$.

With regard to sampling input-output pairs for the generation of gPCE surrogates, we follow the decorrelation procedure described in Algorithm 1 for each sample $\boldsymbol{\Theta}$ that relates to a possible configuration of pore-scale features consistent with the hierarchical nanoporous material in Figure 1. In particular, Algorithm 1 allows seamless, non-intrusive integration with existing codes for the numerical solution of the multiscale forward model presented in Section 2. For the numerical solutions of the forward model, we first solve the closure equations 6, 7, 8, 9 and 10 using a finite element code written in COMSOL and then compute the rate constant $\gamma_{\mathrm{eff}}$ in 4 and effective diffusions ${D_{L}}$ and ${D_{T}}$ in 21 by numerically evaluating the quadrature in 5. For the required gPCE surrogates, given in 24, we select for $\Psi_{i}$ the Askey scheme of hypergeometric orthogonal polynomials ([50]) with $N_{PC}=$626 (i.e., for $N_{p}=4$ parameters we consider polynomials of degree $\kappa_{i}=4$, for each $i=1,\cdots,N_{p}$, in [25](#S3.E25)). These gPCE surrogates are then used to construct KDEs for the desired Darcy-scale flow variables. The KDE approximations for densities below employ the Gaussian-kernels for univariate and multivariate densities, described in [26](#S3.E26) and [27](#S3.E27), respectively, with $N=${10}^{8} where the kernel bandwidths are estimated using a modified Sheather-Jones method ([4]).

As a first experiment, we compare simulations of Darcy-scale flow variables based on the causal relationships encoded by model $P_{1}$ in 20 to simulations based on the independent priors model $P_{0}$ in 15 (cf. Bayesian networks in Figures 4 and 3). Over the narrow range of hyperparameter values given in Table 1, we anticipate qualitative similarities in the resulting Darcy-scale outputs as both $P_{0}$ and $P_{1}$ exhibit statistically uncorrelated parameter distributions over this hyperparameter range as demonstrated by the empirical correlations in Figures 6(a) and 6(b). In Figure 7, we observe that the marginal distributions for Darcy-scale flow variables based on model $P_{1}$ in Figure 7(b) are qualitatively consistent with the marginals based on model $P_{0}$ in Figure 7(a). Likewise, in Figure 8 the joint distributions $({D_{L}},{D_{T}})$, $({D_{T}},\gamma_{\mathrm{eff}})$, and $(\gamma_{\mathrm{eff}},{D_{L}})$ for simulations based on model $P_{1}$ in Figure 8(b) and model $P_{0}$ in Figure 8(a) exhibit qualitative similarities. Although the qualitative nature of the simulated distributions suggest no remarkable difference in the physics over the narrow hyperparameter range, in general the simulations based on model $P_{1}$ follow a different sampling procedure than the simulations based on model $P_{0}$; in contrast to $P_{0}$, the causal relationships embedded in model $P_{1}$ ensure the pore-scale geometries sampled for the numerical experiment are always consistent with hierarchical nanoporous material in Figure 1 even over extended hyperparameter ranges such as Table 2.

A quantitative comparison suggests that the difference in sampling procedures leads to distinct distributions with statistical significance. To compare the densities in Figures 7(a) and 7(b) and in Figures 8(a) and 8(b) quantitatively, we work directly with samples from the gPCEs, employing a two-way statistical test on the equality of distributions. Since the data in Figures 7(b) and 7(a) and in Figures 8(b) and 8(a) appear commensurable, we select a nonparametric Cramér test ([4]) indicated to be sensitive against location alternatives that is applicable to both univariate and multivariate distributions so as to have a consistent presentation. Although the estimated densities in Figure 7(b) and in (Figure 8(b)) are superficially similar to the densities in Figure 7(a) (respectively, Figure 8(a)), the statistical tests each based on $2000$ sample values, summarized in Table 3, reject the hypothesis on equality of distributions with high statistical significance for all but one comparison.

Using model $P_{1}$, one is able to consider numerical experiments over the physical range of hyperparameters in Table 2 that lead to non-trivial correlations as demonstrated in Figure 6(c). Over the physical range of hyperparameters, the marginal and joint densities in Figures 7(c) and 8(c), respectively, for the Darcy-scale flow variable are markedly different from the marginal and joint densities observed in Figures 7(a), 7(b), 8(a) and 8(b). In comparing Figure 7(c) to Figure 7(b), we observe that density for the effective rate constant $\gamma_{\mathrm{eff}}$ becomes more positively skewed and more peaked suggesting less variance in the estimate of $\gamma_{\mathrm{eff}}$ over the physical hyperparameter range. In contrast, we observe in comparing Figure 7(c) to Figure 7(b) that the density for ${D_{T}}$ becomes more uniform in distribution (and more like the distribution of ${D_{L}}$) suggesting that the simulations over the physical hyperparamter range realize a richer variety of transverse diffusions. Similarly, in Figure 8(c) we observe that joint distributions for the Darcy-scale flow variables employing model $P_{1}$ over the physical hyperparameter range demonstrate similar qualitative changes, with the joint densities involving $\gamma_{\mathrm{eff}}$ narrowing in variability and realizing a more variety in the observed transverse diffusion ${D_{T}}$. Altogether, the difference between the densities for the Darcy-scale flow QoIs for model $P_{1}$ suggest that different physics is observed over the narrow versus the physical hyperparameter ranges. In the present context, these differences in physics can have a profound impact on decision support tasks such as experimental design thereby underscoring the importance of including causal relationships in mathematical and statistical modeling processes. In the section that follows, we continue to investigate the impact of causality on uncertainty in Darcy flow variables and QoIs by investigating methodologies for global sensitivity analysis.

Remark 2:

In any of the Bayesian network representations 15, 20, 38 and 23, uncertainty and error in the homogenization can be included by putting a distribution on $\boldsymbol{X}$ that is not trivial, i.e. by replacing 11 with a distribution that captures error in the homogenization map or that compares distributions resulting from different upscaling techniques. Uncertainty and error in other relevant processes, such as the choice of the KDE, might also be incorporated into the Bayesian network for the full statistical model and analyzed similarly. Thus, the Bayesian network PDE formulation provides a systematic framework for building complete predictive models including forward physical models, transitions between scales, and uncertainties in parameters, mechanisms and parameter or model constraints.

5. Global sensitivity analysis and effect ranking

Recall that we are interested in simulating porous media to inform general decision tasks concerning the design of materials with targeted macroscopic properties, such as effective diffusions and sorption constants, through engineered microscopic pore-scale structures. In this context, it is important to analyze the sensitivity of macroscopic QoIs with respect to uncertainties in pore-scale properties. Simulations of macroscopic material performance that are highly sensitivity to distributional changes in microscopic pore-scale processes and structures would undermine the generality of such investigations.

In general, sensitivity analysis is a key component of uncertainty quantification and informs decision tasks such optimal experimental design and the analysis of model robustness, identifiability, and reliability ([20, 43, 41]). Local sensitivity analysis is suited to situations where the (hyper)parameters are known with some confidence and small perturbations are relevant. In contrast, in our application of interest the mapping from pore-scale input distributions to Darcy-scale flow variables is nonlinear, includes several computational steps, and we have no a priori information on the form of the model for the Darcy-scale variables (cf. Section 3.3). Therefore, the present application demands global sensitivity analysis methods that explore the whole space of uncertain input factors ([42]).

Variance-based sensitivity analysis methods, such as Sobol’ indices [44] and total sensitivity indices [17], rank input factors and higher order interactions of input factors in terms of their contributions to the variance of a QoI. In particular, Sobol’ sensitivity indices are used in [47] to analyze the global sensitivity of Darcy-scale QoIs to first and second order interactions among input parameters for the multiscale porous media model considered here under the assumption of independent uninformed priors on the pore-scale distributions. However, such variance-based methods for assessing global sensitivity are not easily applied or interpreted in the case of the dependent input parameters introduced through causal relationships outlined in Section 3. Moreover, we observe that the distributions for macroscopic Darcy-scale flow variables exhibit non-Gaussian behavior (cf. marginal and joint densities in Figures 7 and 8 in Section 4). Methods that rely on moment information alone may be insufficient to capture the full complexity of these interactions.

5.1. Mutual information for global sensitivity analysis

Presently, we employ global sensitivity indices based on information theoretic concepts ([12, 45]) that rely on empirical distributions as opposed to moments. There is a rich literature on moment-independent indices for local sensitivity analysis ([34, 33, 20, 29, 28]) as well as global sensitivity analysis ([11, 7, 25, 8, 26, 10, 39]). In particular, global sensitivity indices based on various information theoretic notions are well established in the literature ([35, 25, 27, 48]). In [27], discrete mutual information based sensitivity indices are demonstrated as an effective tool for discrete probability distributions arising in biochemical reaction networks in systems biology. For our application of interest, sensitivity indices and rankings based on the differential mutual information provide a suitable measure of effect that overcomes the twin challenges of causally related inputs and non-Gaussian QoI. The differential mutual information has explicit connections to more general information theoretic concepts and we provide interpretations of these in the context of uncertainty quantification and global sensitivity analysis.

The differential mutual information between continuous random variables $V$ and $W$,

[TABLE]

quantifies the statistical dependence, that is, the amount of shared information, between $V$ and $W$ provided that all of the densities above exist and the marginals are non-zero. Importantly, 29 applies to dependent random variables $V$ and $W$, such as random variables that share a causal relationship, and we observe that $I(V,W)=0$ if $V$ and $W$ are independent, i.e., if $P(V,W)=P(V)P(W)$. Further, 29 applies to random variables with very general marginal and joint distributions including distributions that are non-Guassian. The differential mutual information is precisely the relative entropy $\operatorname{\mathcal{R}}$ (or Kullback–Leibler divergence) between the joint and marginal densities,

[TABLE]

a well known pseudo-distance used in variational inference ([49, 6]), machine learning ([5]), and model selection ([9]) as well as other areas. Lastly, the differential mutual information $I(V,W)$ is the limiting value of the discrete mutual information (the supremum over all partitions of $V$ and $W$) and therefore shares all of the same properties as its discrete counterpart ([12]); in particular, we will compute the differential mutual information using empirical distributions.

5.2. Estimators for mutual information global sensitivity

indices

For each QoI $g$, we consider a global sensitivity index,

[TABLE]

based on the mutual information between a Darcy-scale QoI $g$ and each uncertain pore-sclae input parameter $\Theta\in\boldsymbol{\Theta}=(\Theta_{1},\dots,\Theta_{N_{p}})$. Intuitively, the index 31 measures the predictability of $g$ given knowledge of $\Theta$ through the discrepancy between the joint density and the product of the marginal densities appearing in 29.

We estimate the sensitivity index 31 using a MC approximation that relies on empirically estimated distributions, i.e. KDEs for the density functions of Darcy-scale QoIs. Indeed, a benefit to using the mutual information is the availability of methods relying on plug-in estimators for 29 with corresponding numerical analysis (including [21, 32, 3]). Recall that we obtain approximate densitities $\bar{f}$ using univariate and multivariate KDEs, in 26 and 27, respectively, that are in turn obtained using the gPCE surrogates $\hat{g}$ defined in 24. Due to the availability of these surrogates, we do not sample input-output pairs as suggested by the form of 29 but instead consider the equivalent representation,

[TABLE]

Thus we compute the statistical estimator $\widehat{S}_{\Theta}(g)\approx S_{\Theta}(g)$,

[TABLE]

based on a MC approximation of 32 using KDEs as plug-in estimates where appropriate.

In the left-hand side of Table 4, we report the value of the statistical estimator $\widehat{S}_{\Theta}(g)$ related to the numerical experiments presented in Section 4.2 concerning the probabilistic model $P_{1}$ in 20 with causal inputs (see Figure 4) over both the narrow and physical hyperparameter ranges in Tables 1 and 2. We are interested in the global sensitivity of the Darcy-scale QoIs $g={D_{L}}$, $g={D_{T}}$, and $g=\gamma_{\mathrm{eff}}$ with respect to each of input parameters $\Theta\in\boldsymbol{\Theta}=(\Theta_{R},\Theta_{\theta},\Theta_{d},\Theta_{l})$ representing the pore-scale features identified in the hierarchical nanoporous material in Figure 1. The computed estimators 33 are based on first generating random variables $\Theta^{k}$ and $\hat{g}^{k}$ for $k=1,\dots,N=$1\text{\times}{10}^{7}$$ using the respective gPCE surrogate obtained by the workflow outlined in Section 4. These samples are used to form the respective KDE on an $\eta$-grid of size 128 by 128 which in turn are then used to form the plug-in quantity

[TABLE]

where we define $\log(\frac{0}{0})\cdot\frac{0}{0}=0$ (re-using $1\text{\times}{10}^{5}$ of the samples $(\Theta^{j},\hat{g}^{j})$ generated during density estimation). The sensitivity index is then estimated using the MC estimator, $\widehat{S}_{\Theta}(\hat{g})=\frac{1}{M}\sum_{j=1}^{M}X^{j}$ with $M=$1\text{\times}{10}^{5} samples. We observe that the $M=$1\text{\times}{10}^{5} samples utilized for the MC estimator is suggested to be sufficient by the convergence demonstrated in Figure 9 for the experiment with model $P_{1}$ over the range Table 1; similar convergence observations were made for the other experiments.

Remark 3:

The differential mutual information 29 can be expressed as an expected value of the relative entropy between conditional and marginal distributions,

[TABLE]

for the input-output pair $(g(Y(\boldsymbol{\Theta})),\Theta)$ provided $P(g,\Theta)=P(g,\mid\Theta)P(\Theta)$, i.e., all the densities exist. In [39], a family of sensitivity measures is given by replacing $\operatorname{\mathcal{R}}$ in 34 with a general class of Csiszár $\varphi$-divergences ([2, 13, 23]); here we demonstrate an implementation using gPCE surrogates that requires the representation 32 and restrict our attention to the differential mutual information owing to the clear interpretation as a measure of effect in terms of statistical dependence and shared information.

Remark 4:

One can also consider higher order effects that include interactions between a subset of parameters using the conditional differential mutual information that takes the form of conditional expectations of the relative entropy between joints and marginals, as in [27]. These higher order interactions can be interpreted similarly in terms of dependence and shared information.

5.3. Ranking impact of uncertainty in correlated pore-scale

inputs on Darcy-scale QoIs

Using the global sensitivity index 31 we form the ranking,

[TABLE]

of the relative contribution of each pore-scale parameter $\Theta\in\boldsymbol{\Theta}$ to the global sensitivity of each Darcy-scale variable $g$. We then obtain the estimate $\widehat{r}_{\Theta}(g)$ reported in the right-hand side of Table 4, by using the estimator $\widehat{S}$ in 33 in place of $S$. The values $\widehat{r}$ are also displayed graphically in Figure 10 with error bars that indicate the relative error associated with plus or minus two standard deviations of the computed sensitivity index to provide an indication of confidence in the ranking.

In Figure 10(b), we observe that the rankings suggested by the model $P_{1}$ in 20 with causal inputs in Figure 4 over the narrow hyperparameter range in Table 1 are consistent with the Sobol’ index rankings in Figure 10(a) that were obtained in [47] for the model $P_{0}$ with independent priors. As observed in Section 4.2 this is to be expected due to the similarity in the correlation structures in Figures 6(a) and 6(b) over the narrow hyperparameter range. The $\theta$, which is related to the mesopore radius, is the most influential parameter for both the longitudinal diffusion ${D_{T}}$ and the transverse diffusion ${D_{T}}$. The mesopore radius $R$, which is related to the size of the fluid-solid interface in Figure 1, is the most influential parameter for the sorption rate constant $\gamma_{\mathrm{eff}}$. In general, it is important to interpret these sensitivity rankings in the context of the hyperparameter range; the ranking $\widehat{r}_{\Theta_{d}}$ in Figure 10(b) is likely to be small as the range for $d$ in Table 1 is very narrow (i.e. the Darcy-scale QoIs are insensitive over the narrow range of admissible $d$ values).

In contrast, the rankings in Figure 10(c) for the model $P_{1}$ over the physical hyperparameter range in Table 2 demonstrates considerably different rankings to Figures 10(b) and 10(a) thereby highlighting once again the impact of causal relationships and the Bayesian network PDE modeling approach. The rankings in Figure 10(c) indicate that while the $\theta$ and the $R$ are still the most influential parameters for, respectively, the longitudinal diffusion, ${D_{L}}$, and sorption rate constant, $\gamma_{\mathrm{eff}}$, it is $d$, which is related to the diameter of the nanotube, that is the most influential parameter for the transverse diffusion ${D_{T}}$ through nanotubes. Thus, for the simple hierarchical nanoporous material in Figure 1, the rankings in Figure 10(c) using model $P_{1}$ over the physical hyperparameter range reflect our expectations of the physics better than the experimental observations with respect to models over the narrow range. Importantly, the rankings in Figure 10 are moment-independent and therefore suitable for the non-Gaussian behavior of the Darcy-scale QoIs observed in Figures 7 and 8. Moreover, as the rankings are based on the mutual information, they can be interpreted as ranking the impact of the pore-scale parameter (whether correlated or not) on the total uncertainty in the Darcy-scale quantity of interest.

6. Alternative probabilistic models reflect different

engineering and design causal relationships

We view the construction of the full statistical model for the multiscale porous media system as mirroring engineering processes and design work-flows. From this perspective, it may be desirable to build models containing different causal relationships among the pore-scale parameters than those previously considered. Recall that the model $P_{1}$ in 20 with causal priors given by the Bayesian network in Figure 4 is related to the design of nanotunnels/mesopores, i.e. the pore-scale features $R$ and $\theta$ in the hierarchical nanoporous media in Figure 1. In contrast, if it is desirable or possible to choose three aspects, such as the features $R$, $\theta$, and $l$, independently while constraining the only remaining parameter, $d$, then one obtains a second model represented by the Bayesian network in Figure 11(a).

The Bayesian network in Figure 11(a) corresponds to placing all of the constraints on $\Theta_{d}$. That is, choosing a subset of independent parameters $\{\Theta_{R},\Theta_{\theta},\Theta_{l}\}$ and assuming uniform priors,

[TABLE]

constrains the distribution of $\Theta_{d}$. Specifically, in order for the choice of $d$ to be consistent with the features depicted in Figure 1, then (i) the nanotube diameter $d<2R\cos\theta$ and (ii) if the gap between vertically mesopores is zero than the nanotube diameter must be less than $d<\sqrt{4lR-l^{2}}$ (i.e. to avoid nanotube “goiters”, see further Figure 5). The corresponding conditional distribution is then given by

[TABLE]

where again we assume a uniform model. The form of the full statistical model that includes the causal relationships in Figure 11(a) is then given by

[TABLE]

where the distributions for conditionals and priors above are given in 36 and 37.

We observe that the Bayesian networks in Figure 4 and Figure 11(a) are not equivalent. Over the physical range of hyperparameters in Table 2, we observe that correlation structure in Figure 11(d) is distinct from Figure 8(c) and thus the joint density $P(\boldsymbol{\Theta})$ for the Bayesian networks in Figure 11(a) and Figure 4 are different. Although it may be possible to select priors to ensure equality between the distributions $P_{1}$ and $P_{2}$, this is not in general the goal. Further, we can compare the Darcy-scale QoIs for the model $P_{2}$ to the model $P_{1}$ in Figure 11(c) over the physical hyperparameter ranges in Table 2, i.e. by comparing Figure 11(c) to Figure 7(c) and Figure 11(d) to Figure 8(c). Importantly, incorporating constraints directly into the models $P_{1}$ and $P_{2}$ allows us to ensure that we sample from a joint distribution that ensures realistic geometries that are consistent with the hierarchical nanoporous material in Figure 1 over the physical range of hyperparameters. On the one hand, the densities in Figure 7(c) (and Figure 8(c)) are qualitatively similar to those in Figure 11(c) (respectively, Figure 11(d)) owing to the closeness of the input densities $P(\boldsymbol{\Theta})$, see the empirical correlations in Figures 6(c) and 11(b). On the other hand, we observe that the models $P_{1}$ and $P_{2}$ lead to distinct distributions on Darcy-scale flow variables as the Cramér tests reported in Table 5 reject the hypothesis of equality of the empirical distributions with high statistical significance. We emphasize that we do not give here a method for selecting among various models but instead a collection of tools for breaking the stochastic modeling task into smaller, manageable components that enables a systematic way of building a full statistical model that incorporates engineering design constraints.

7. Conclusions

The sensitivity of Darcy-scale observables to changes in pore-scale properties and rigorous quantification of the uncertainty in predictions are some of the least studied aspects of multiscale models of flow and transport in porous materials. Our analysis leads to the following major conclusions.

• •

Causal relationships are natural and stem from physical or chemical constraints, engineering design, and expert knowledge. These relationships exist between model parameters, scales, and model components in multi-scale and multi-physics models. In the context of hierarchical nanoporous materials, we observe that causal relationships suggested by geometrical structural constraints among microscopic features yield non-trivial correlations among pore-scale model parameters over physical ranges.

• •

We incorporate causal relationships, and thereby correlations, in a unified random, multiscale PDE model using structured probabilistic graphical models, in this case Bayesian networks. This perspective ultimately gives rise to Bayesian network (random) PDEs.

• •

Due to causal relationships and resulting correlations between model parameters, global sensitivity analysis is not straightforward. Furthermore, predictive PDFs of QoIs are not necessarily Gaussian or otherwise of known analytical form. For these two reasons, it is necessary to depart from moment-based sensitivity analysis, e.g. ANOVA methods, and deploy PDF-based methods such as mutual information.

• •

The proposed mutual information global sensitivity indices for Bayesian network PDE and corresponding Darcy-scale QoIs yield parameter rankings that are consistent with our understanding of the physics of a hierarchical nanoporous material. We demonstrate that correlations stemming from causal relationships turn out to be important in determining the most influential model parameters/mechanisms and impact predictions of QoIs.

The structure of Bayesian networks or PGMs for the pore-scale parameter correlations can in general be learned from experimental or simulated data rather than based on natural structural constraints as in this work and moreover they need not adopt causal relationships (Bayesian networks) that mirror engineering design processes. One future direction includes learning the structure of the Bayesian networks from available data and inferring the distributions of pore-scale features directly. The flexibility of the Bayesian network PDE approach for incorporating uncertainty into physical models will also enable explorations of model-form uncertainty that advance traditional uncertainty quantification techniques. Another future direction for research is to apply hybrid information divergences (e.g. [16]) to bound families of model predictions based on various pore-scale parameter representations and based on different upscaling techniques (i.e. different probabilistic models) to address questions of model-form uncertainty and tracking different levels of fidelity attached to system components.

Acknowledgments

This research was supported in part by the Department of Energy (DOE) Office of Science under grant 0000241868, by the Air Force Office of Scientific Research under grants FA9550-17-1-0417 and FA9550-18-1-0214, and by National Science Foundation under grant 1606192.