Numerical approximation of poroelasticity with random coefficients using Polynomial Chaos and Hybrid High-Order methods

TL;DR

This paper develops a numerical method combining Polynomial Chaos and Hybrid High-Order discretizations to efficiently approximate the behavior of poroelastic materials with uncertain properties, validated through convergence and uncertainty propagation studies.

Contribution

It introduces a non-intrusive Polynomial Chaos approach integrated with Hybrid High-Order methods for stochastic poroelasticity problems, supporting complex meshes and arbitrary approximation orders.

Findings

Convergence of Polynomial Chaos approximation with sparse spectral projection.

Effective propagation of input uncertainty to the solution.

Validation through numerical experiments on an injection-extraction problem.

Abstract

In this work, we consider the Biot problem with uncertain poroelastic coefficients. The uncertainty is modelled using a finite set of parameters with prescribed probability distribution. We present the variational formulation of the stochastic partial differential system and establish its well-posedness. We then discuss the approximation of the parameter-dependent problem by non-intrusive techniques based on Polynomial Chaos decompositions. We specifically focus on sparse spectral projection methods, which essentially amount to performing an ensemble of deterministic model simulations to estimate the expansion coefficients. The deterministic solver is based on a Hybrid High-Order discretization supporting general polyhedral meshes and arbitrary approximation orders. We numerically investigate the convergence of the probability error of the Polynomial Chaos approximation with respect to the level of the sparse grid. Finally, we assess the propagation of the input uncertainty onto the solution considering an injection-extraction problem.

Equations

$${\mathrm{d}}_t({\rho }_{\mathrm{f}}\phi )+\vec{\nabla }\cdot ({\rho }_{\mathrm{f}}\phi \vec{v})={\rho }_{\mathrm{f}}g,$$

$${\mathrm{d}}_t{\rho }_{\mathrm{f}}=\frac{{\rho }_{\mathrm{f}}}{K_{\mathrm{f}}}{\mathrm{d}}_{t}p,$$

$${\mathrm{d}}_t\phi =\frac{1}{M}{\mathrm{d}}_{t}p+\alpha {\mathrm{d}}_t(\vec{\nabla }\cdot \vec{u}),$$

$$\alpha :-1-\frac{K_{\mathrm{d}}}{K_{\mathrm{m}}},M:-\frac{K_{\mathrm{m}}}{\alpha -\phi },$$

$${\mathrm{d}}_t({\rho }_{\mathrm{f}}\phi )=\phi {\mathrm{d}}_t{\rho }_{\mathrm{f}}+{\rho }_{\mathrm{f}}{\mathrm{d}}_t\phi ={\rho }_{\mathrm{f}}\left(c_0{\mathrm{d}}_{t}p+\alpha {\mathrm{d}}_t(\vec{\nabla }\cdot \vec{u})\right),$$

$$c_0:-\frac{\alpha -\phi }{K_{\mathrm{m}}}+\frac{\phi }{K_{\mathrm{f}}}.$$

$$c_0{\mathrm{d}}_{t}p+\alpha {\mathrm{d}}_t(\vec{\nabla }\cdot \vec{u})+\vec{\nabla }\cdot (\phi \vec{v})=g.$$

$$\phi \vec{v}=-\frac{\mathbb{K}}{{\mu }_{\mathrm{f}}}\vec{\nabla }p,$$

$$c_0{\mathrm{d}}_{t}p+\alpha {\mathrm{d}}_t(\vec{\nabla }\cdot \vec{u})-\vec{\nabla }\cdot (\vec{\kappa }\vec{\nabla }p)=g.$$

$$-\vec{\nabla }\cdot \widetilde{\text{matr}}\sigma =\vec{f},$$

$$\widetilde{\text{matr}}\sigma =\text{matr}\sigma ({\vec{\nabla }}_{\mathrm{s}}\vec{u})-\alpha p{\vec{I}}_d.$$

$$\text{matr}\sigma (\overrightarrow{ϵ})=2\mu \overrightarrow{ϵ}+\lambda \mathrm{tr}(\overrightarrow{ϵ}){\vec{I}}_d=2\mu \overrightarrow{\mathrm{dev}}(\overrightarrow{ϵ})+K\mathrm{tr}(\overrightarrow{ϵ}){\vec{I}}_d,$$

$$\mathrm{tr}(\overrightarrow{ϵ}):-\sum _{i=1}^d{ϵ}_{ii},\overrightarrow{\mathrm{dev}}(\overrightarrow{ϵ}):-\overrightarrow{ϵ}-\frac{\mathrm{tr}(\overrightarrow{ϵ}){\vec{I}}_d}{d},$$

$$-\vec{\nabla }\cdot \left(2\mu {\vec{\nabla }}_{\mathrm{s}}\vec{u}+(\lambda \vec{\nabla }\cdot \vec{u}-\alpha p){\vec{I}}_d\right)=\vec{f}.$$

$${\alpha }^2=c_0(K-K_{\mathrm{d}}).$$

$$K_{\mathrm{m}}=\frac{K{c}_0-{\alpha }^2}{c_0(1-\alpha )}.$$

$$K{c}_0^2-\left(\alpha +\alpha \phi +\frac{\phi K}{K_{\mathrm{f}}}-\phi \right)c_0+\frac{{\alpha }^2\phi }{K_{\mathrm{f}}}=0.$$

$$\frac{3\phi }{2+\phi }\le \alpha .$$

$$\mathcal{D}$$

$$={\phi }^2{\left(\text{nicefrac}K{K}_{\mathrm{f}}-1\right)}^2+2\alpha \phi \left(1+\phi -2\alpha \right)\left(\text{nicefrac}K{K}_{\mathrm{f}}-1\right)+{\alpha }^2(1-\phi {)}^2$$

$$={\left[\phi \left(\text{nicefrac}K{K}_{\mathrm{f}}-1\right)+\alpha (1+\phi -2\alpha )\right]}^2+{\alpha }^2(1-\phi {)}^2-{\alpha }^2(1+\phi -2\alpha {)}^2$$

$$={\left[\phi \left(\text{nicefrac}K{K}_{\mathrm{f}}-1\right)+\alpha (1+\phi -2\alpha )\right]}^2+4{\alpha }^2(1-\alpha )(\alpha -\phi ).$$

$$\mathbb{E}(h):-{\int }_{\Theta }h(\theta )\mathrm{d}\mathcal{P}(\theta ).$$

$$L_{\mathcal{P}}^2(\Theta ):-\left\{h:\Theta \to \mathbb{R};\mathbb{E}(h^2)<+\infty \right\}.$$

$$-\vec{\nabla }\cdot \text{matr}\sigma (\vec{x},\theta ,t)+\vec{\nabla }(\alpha (\vec{x},\theta )p(\vec{x},\theta ,t))$$

$${\mathrm{d}}_t\varphi (\vec{x},\theta ,t)-\vec{\nabla }\cdot (\vec{\kappa }(\vec{x},\theta )\vec{\nabla }p(\vec{x},\theta ,t))$$

$$\vec{u}(\vec{x},\theta ,t)$$

$$\text{matr}\sigma (\vec{x},\theta ,t)\vec{n}+\alpha (\vec{x},\theta )p(\vec{x},\theta ,t))\vec{n},$$

$$p(\vec{x},\theta ,t)$$

$$\vec{\kappa }(\vec{x},\theta )\vec{\nabla }p(\vec{x},\theta ,t)\cdot \vec{n}$$

$$\varphi (\vec{x},\theta ,0)$$

$$\text{matr}\sigma (\vec{x},\theta ,t)$$

$$\varphi (\vec{x},\theta ,t)$$

$${\int }_D{\varphi }_0(\cdot )=0,{\int }_{D}g(\cdot ,t)=0,\text{ and }{\int }_{D}p(\cdot ,\cdot ,t)=0\forall t\in (0,t_F).$$

$$L_{\mathcal{P}}^2(\Theta ;V(D)):-\left\{v:D\times \Theta \to \mathbb{R};\Vert v{\Vert }_{L_{\mathcal{P}}^2(\Theta ;V(D))}:-\mathbb{E}{\left(\Vert v{\Vert }_{V(D)}^2\right)}^{\frac{1}{2}}<+\infty \right\}.$$

$$\vec{U}$$

Full text

Key words. Biot problem, poroelasticity, Uncertainty Quantification, Polynomial Chaos expansions, Pseudo-Spectral Projection methods, Hybrid High-Order methods

AMS subject classification. 65C30, 65M60, 65M70, 35R60, 76S05

1 Introduction

This work aims at the study poroelasticity problems where the model coefficients are uncertain. The interest of this type of hydro-mechanical coupled models is particularly manifest in geosciences applications [39, 41, 51, 50], where subsurface fluid flows can induce a deformation of the rock matrix. We rely here on the linear Biot model [6, 61] that describes Darcean fluid flows in saturated porous media under the assumptions of small deformations and small variations of the porosity and the fluid density. This model depends on physical parameters that are often poorly known, justifying a stochastic description, due spatial heterogeneities, measurement inaccuracies, and sometimes the ill-posedness of inverse problems inherent to parameter estimation techniques. Although there is an extensive literature on poroelasticity models (we mention, in particular, the comprehensive textbooks [19, 23, 64]) and their numerical approximation, to our knowledge few works have addressed the impact of uncertainty in Biot’s formulation of hydromechanical coupling. In [13, 38] the authors include uncertainty in one-dimensional consolidation analysis by incorporating heterogeneity in the consolidation coefficient. The stochastic consolidation model was extended in [53] to nonlinear uncertain soil parameters. Variabilities in the initial pore pressure and the heterogeneous hydraulic mobility have been considered respectively in [21] and [28]. We also mention [22], where a stochastic Galerkin approach is proposed to solve the poroelasticity equations with randomness in all material parameters and tested on a one-dimensional problem.

Uncertainty Quantification (UQ) methods have been developed in the last decades to take into account the effect of random model inputs on the model predictions. Compared to a simple deterministic simulation, advanced UQ methods can characterize an uncertain model prediction, computing its statistical moments (e.g., mean and variance), probability distributions, and performing global sensitivity analyses. Among the techniques designed for UQ in numerical models, the stochastic spectral methods have received considerable attention. The principle of these methods is to decompose random quantities on suitable stochastic approximation bases. In particular, Polynomial Chaos (PC) expansions represent quantities as finites series in stochastic orthogonal polynomial bases. PC were initially introduced by Wiener [65] and applied by Ghanem and Spanos [35] to solid mechanics and by Le Maître and Knio [44] to fluid mechanics. PC expansion have been used to treat a large variety of problems, including elliptic models (see, e.g., [1, 49]), flow and transport in porous media (see, e.g., [34, 33]), thermal problems (see, e.g., [37, 46]), and hyperbolic systems (see, e.g., [62]). In this work, we consider a non-intrusive spectral projection method, that only requires the resolution of several deterministic problems to construct the stochastic spectral expansion of quantities of interest (“black box” approach). The numerical solution of the deterministic poroelasticity coupled system requires a discretization method able to

(i) treat complex geometry with polyhedral meshes and nonconforming interfaces,

(ii) handle possible heterogeneities of the poromechanical parameters, and

(iii) prevent localized pressure oscillations arising in the case of low-permeable and low-compressible porous media.

We choose to discretize the poroelasticity problem using the fully coupled method developed in [8], which meets all of these requirements. Therein, the elasticity operator is discretized using the Hybrid High-Order method of [25] (c.f. also [10, 24]), while the Darcy operator relies on the Symmetric Weighted Interior Penalty method of [26].

The contribution of this work is threefold. First, a probabilistic framework is introduced to study the Biot model with uncertain coefficients. Special attention is given to defining an almost surely physically admissible set of stochastic poroelastic parameters by taking into account the dependences between the poroelastic coefficients. Second, the well-posedness of the stochastic Biot model is proven at the continuous level. Third, a non-intrusive PC approach is implemented in order to investigate the effect of the random poroelastic coefficients on the displacement and the pressure fields.

The material is organized as follows. In Section 2 we present the deterministic poroelasticity model and identify the relations among the poromechanical coefficients allowing to express the constrained specific storage coefficient as a function of the other parameters and primary variables. In Section 3 we present the linear poroelasticity problem with random coefficients in strong and weak form. We then give the probabilistic assumptions on the random model coefficients yielding the well-posedness of the stochastic variational problem. In Section 4 we propose an uncertainty model for the poroelastic coefficients. The sparse Pseudo-Spectral Projection (PSP) method is used to compute the PC expansions of the stochastic solutions from a finite set of deterministic problems solved using the method of [8]. Finally, we present a complete panel of numerical results in Section 5 to illustrate the performance of the method.

2 The Biot model

In this section we introduce the model problem and present some preliminary relations among the coefficients characterizing the porous medium. This linear poroelasticity model, usually referred to as the Biot model, consists of two coupled governing equations, one describing the mass balance of the fluid and the other expressing the mechanical equilibrium of the porous medium.

2.1 Fluid mass balance

We consider a spatial domain in dimension $d\geq 1$. The fluid mass conservation law in a fully saturated porous medium reads

[TABLE]

where ${\rm d}_{t}$ denotes the time derivative, $\rho_{\rm f}\ [{\rm kg}/\rm{m}^{3}]$ is the fluid density, $\varphi\ [-]$ the soil porosity, $\vec{v}\ [\rm{m}\rm{s}^{-1}]$ the velocity field, and $g\ [\rm{s}^{-1}]$ the fluid source. The equation of state for slightly compressible fluids (cf. [19, Section 3.1] and [67]) is

[TABLE]

with $p\ [{\rm Pa}]$ the pore pressure and $K_{\rm f}>0\ [{\rm Pa}]$ the fluid bulk modulus. Under the assumptions of isotropic and isothermal conditions, infinitesimal strains, and small relative variations of porosity, following [19, Section 4.1], the change in the porosity caused by the fluid pressure and the mechanical displacement field $\vec{u}\ [\rm{m}]$ is

[TABLE]

where $\vec{\nabla}{\cdot}\vec{u}$ is the divergence of the displacement field. We rely on the definitions of the Biot–Willis coefficient $\alpha\in(0,1]\ [-]$ and the Biot tangent modulus $M>0\ [{\rm Pa}^{-1}]$ given in [19, 23]:

[TABLE]

with $K_{\rm d},K_{\rm m}>0\ [{\rm Pa}]$ denoting the bulk moduli of the drained medium and the solid matrix, respectively. The coefficient $\alpha$ quantifies the amount of fluid that can be forced into the medium by a variation of the pore volume at a constant fluid pressure, while $M$ measures the amount of fluid that can be forced into the medium by pressure increments due to the compressibility of the structure. The case of a solid matrix with incompressible grains ($K_{\rm m}\to+\infty$) corresponds to the limit value $\nicefrac{{1}}{{M}}=0$. From (2) and (3), it follows that the variation of fluid content in the medium is given by

[TABLE]

where, according to (2), (3), and (4) the constrained specific storage coefficient is defined by

[TABLE]

Using (1) and assuming that the fluid density $\rho_{\rm f}$ is uniform in the medium, we obtain

[TABLE]

The fluid velocity $\vec{v}$ is related to the pore pressure through the well-known Darcy law (see for instance [23, Section 4.1.2]). Consistent with the small perturbations hypothesis adopted here, Darcy’s law can be considered in its simplest form

[TABLE]

where $\mathbb{K}\ [\rm{m}^{2}]$ is the tensor-valued intrinsic permeability of the medium and $\mu_{\rm f}\ [{\rm Pa}\cdot\rm{s}]$ is the fluid viscosity. Thus, introducing the hydraulic mobility $\vec{\kappa}\coloneq\nicefrac{{\mathbb{K}}}{{\mu_{\rm f}}}$, the mass conservation equation (6) becomes

[TABLE]

For the sake of simplicity, we also assume in what follows that the mobility field is isotropic and strictly positive, namely $\vec{\kappa}=\kappa\vec{I}_{d}$, with $0<\kappa<+\infty$ and $\vec{I}_{d}$ the identity matrix.

2.2 Momentum balance

The momentum conservation equation under the quasi-static assumption, namely when the inertia effects in the elastic structure are negligible, reads

[TABLE]

where $\widetilde{\matr}{\sigma}\ [{\rm Pa}]$ is the total stress tensor and $\vec{f}\ [{\rm Pa}/\rm{m}]$ is the loading force (e.g., gravity). According to Terzaghi’s decomposition [61], the stress tensor is modeled as the sum of an effective term (mechanical effect) and a pressure term (fluid effect), i.e.,

[TABLE]

The symmetric part of the gradient of the displacement field $\vec{\nabla}_{\rm s}\vec{u}$ measures the strain accordingly to the small deformation hypothesis. In the context of linear isotropic poroelasticity, the mechanical behavior of the soil is described through the Cauchy strain-stress relation defined, for all deformation tensors $\vec{\epsilon}\in\mathbb{R}^{d\times d}_{\rm sym}$, by

[TABLE]

where $\mu>0\ [{\rm Pa}]$ and $\lambda>0\ [{\rm Pa}]$ are Lamé’s coefficients, $K\coloneq\nicefrac{{2\mu}}{{d}}+\lambda\ [{\rm Pa}]$ is the bulk modulus, and

[TABLE]

are the trace and deviator operator, respectively. Incidentally, as noted in [4, 7], physical and experimental investigations suggest that the mechanical behavior of porous solids may be nonlinear. More general stress-strain relations could be used in place of the linear law (10), as done in [10, 11]. Inserting (10) and (9) into (8), leads to

[TABLE]

2.3 Constrained specific storage coefficient

In order to perturb the Biot model coefficients, we need to investigate the dependences between the poroelastic coefficients. We propose here to express the specific storage coefficient $c_{0}$ as a function of other physical parameters. For that purpose, we introduce the Gassmann equation (cf. [30, 50]), that relates the bulk moduli $K$ and $K_{\rm d}$ to the Biot–Willis and storage coefficients $\alpha,c_{0}$:

[TABLE]

Using (4) to express the drained bulk modulus $K_{\rm d}$ as a function of $\alpha$ and $K_{\rm m}$, equation (12) gives

[TABLE]

Plugging the previous equation into (5) and rearranging yields the following quadratic equation for $c_{0}$:

[TABLE]

The parameters $\varphi,K_{\rm f},K$, and $\alpha$ can then be used to define $c_{0}$ through the previous equation. Some conditions need to be prescribed in order to avoid non-physical solutions $c_{0}$. Owing to the definition of the Biot–Willis coefficient, it holds that $0\leq\varphi\leq\alpha\leq 1$. As observed in [67, Section 3], a stricter lower bound can be imposed, namely

[TABLE]

Lemma 1** (Existence).**

Let $K_{\rm f},K>0$, and $\varphi,\alpha$ satisfying condition (15). Then, there exists $c_{0}\in\mathbb{R}^{+}$ solution of (14).

Proof.

We prove existence by assessing the positivity of the discriminant $\mathcal{D}$ associated to (14). Computing $\mathcal{D}$ and rearranging, leads to

[TABLE]

Since $0\leq\varphi<\alpha\leq 1$ owing to (15), the second term in the previous sum is positive and, as a result, $\mathcal{D}\geq 0$. ∎

Remark 2*.*

The previous lemma yields the existence of two real solutions $c_{0}^{-}\leq c_{0}^{+}$. We consider $c_{0}^{+}$ as the unique solution to (14) because, for admissible values of $(K,K_{\rm f},\varphi,\alpha)$, we might have $c_{0}^{-}<\nicefrac{{\varphi}}{{K_{\rm f}}}$, violating (5). For instance, this is the case if $\alpha=\nicefrac{{2\varphi}}{{(1+\varphi)}}$ for any $\phi\leq\nicefrac{{3}}{{4}}$ and $K,K_{\rm f}>0$.

3 The Biot problem with random coefficients

In this section, we introduce the extension of the linear poroelasticity problem in the case of uncertain poroelastic coefficients. We discuss the assumptions on the random coefficients, derive a weak formulation of the resulting stochastic problem, and prove its well-posedness.

3.1 Strong and weak formulations

Let $n\in\mathbb{N}$ and denote by $X\subset\mathbb{R}^{n}$ a measurable set. Spaces of functions, vector fields, and tensor fields defined over $X$ are respectively denoted by italic capital, boldface Roman capital, and special Roman capital letters. Thus, for example, $L^{2}(X)$, $\textbf{L}^{2}(X)$, $\mathbb{L}^{2}(X)$ denote the spaces of square integrable functions, vector fields, and tensor fields over $X$, respectively. For any $m\in\mathbb{N}$, we denote by $H^{m}(X)$ the usual Sobolev space of functions that have weak partial derivatives of order up to $m$ in $L^{2}(X)$, with the convention that $H^{0}(X)\coloneq L^{2}(X)$, while $C^{m}(X)$ and $C_{\rm c}^{\infty}(X)$ denote, respectively, the usual spaces of $m$-times continuously differentiable functions and infinitely continuously differentiable functions with compact support on $X$.

We introduce an abstract probability space $(\Theta,\mathcal{B},\mathcal{P})$, where $\Theta$ is the set of possible outcomes, $\mathcal{B}$ a $\sigma$-algebra of events, and $\mathcal{P}:\mathcal{B}\to[0,1]$ a probability measure. For any random variable $h:\Theta\to\mathbb{R}$ defined on the abstract probability space, the expectation of $h$ is

[TABLE]

We assume hereafter that all random quantities are second-order ones, namely they belong to

[TABLE]

The Biot problem consists in finding a vector-valued displacement field $\vec{u}$ and a scalar-valued pressure field $p$ solutions of (7) and (11). We now consider the stochastic version corresponding to the case of random poroelastic coefficients. We assume that $D\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, is a (deterministic) bounded connected polyhedral domain with boundary $\partial D$ and outward normal $\vec{n}$. In order to close the problem, we enforce boundary conditions corresponding to a medium that is clamped on $\Gamma_{D}\subset\partial D$, traction-free on $\Gamma_{N}\coloneq\partial{D}\setminus\Gamma_{D}$, permeable with free drainage on $\Gamma_{d}\subset\partial D$, and impermeable on $\Gamma_{n}\coloneq\partial{D}\setminus\Gamma_{d}$, as well as a deterministic initial condition prescribing the initial fluid content $\phi_{0}$. For a given finite time $t_{F}>0$, the resulting problem is given by

[TABLE]

For the sake of simplicity, we do not consider the case of $\Gamma_{D}$ with zero $(d-1)$-dimensional Hausdorff measure. This situation simply requires an additional compatibility condition on $\vec{f}$ and the prescription of the rigid-body motions of the medium. We remark that inhomogeneous and possibly random boundary conditions, as well as random loading, can also be considered up to minor modifications.

Before giving the variational formulation of problem (16), we introduce some notations. For a given vector space $V(D)$ of real-valued functions on $D$, equipped with the norm $\|\cdot\|_{V(D)}$, we define the Bochner space of second-order random fields by

[TABLE]

Since the Biot problem is unsteady, we also need to introduce, for a given vector space $W$, the Bochner space $L^{2}((0,t_{\rm F});W)$ of square integrable $W$-valued functions of the time interval $(0,t_{\rm F})$. Similarly, the Sobolev space $H^{1}((0,t_{\rm F});W)\subset L^{2}((0,t_{\rm F});W)$ is spanned by $W$-valued functions having square integrable first-order weak derivative. In what follows, for all $(\vec{x},\theta,t)\in D\times\Theta\times(0,t_{F}]$, $v\in L^{2}_{\mathcal{P}}(\Theta;V(D))$, and $w\in L^{2}((0,t_{\rm F});L^{2}_{\mathcal{P}}(\Theta;V(D)))$, we use $v(\vec{x},\theta)$ and $w(\vec{x},\theta,t)$ as shorthand notations for $(v(\theta))(\vec{x})$ and $(w(t))(\vec{x},\theta)$, respectively.

At each time $t\in(0,t_{\rm F}]$, the natural functional spaces for the random displacement field $\vec{u}(t):D\times\Theta\to\mathbb{R}^{d}$ and pressure field $p(t):D\times\Theta\to\mathbb{R}$ taking into account the boundary conditions (16c)–(16f) are, respectively,

[TABLE]

with $\,\textbf{H}_{0,\Gamma_{D}}^{1}(D)\coloneq\left\{\vec{v}\in\textbf{H}^{1}(D):{\vec{v}}_{|\Gamma_{D}}=\vec{0}\right\}$, $\,{H}^{1}_{0,\Gamma_{d}}(D)\coloneq\left\{q\in H^{1}(D):q_{|\Gamma_{d}}=0\right\}\,$, and

[TABLE]

We consider the following weak formulation of problem (16): For a loading $\vec{f}\in L^{2}((0,t_{\rm F});\textbf{L}^{2}(D))$, a fluid source $g\in L^{2}((0,t_{\rm F});L^{2}(D))$, and an initial datum $\phi_{0}\in L^{2}(D)$ that verify (16j), find $\vec{u}\in L^{2}((0,t_{\rm F});\vec{U})$ and $p\in L^{2}((0,t_{\rm F});P)$ such that, for all $\vec{v}\in\vec{U}$, all $q\in P$, and all $\psi\in C_{\rm c}^{\infty}((0,t_{\rm F}))$

[TABLE]

where $(\cdot,\cdot)_{D}$ denotes the usual inner product in $L^{2}(D)$ and the bilinear forms $a:\vec{U}\times\vec{U}\to\mathbb{R}$, $b:\vec{U}\times P\to\mathbb{R}$, and $c:P\times P\to\mathbb{R}$ are defined such that, for all $\vec{v},\vec{w}\in\vec{U}$ and all $q,r\in P$,

[TABLE]

Above, we have introduced the Frobenius product such that, for all $\matr{\tau},\matr{\eta}\in\mathbb{R}^{d\times d}$, $\matr{\tau}:\matr{\eta}\coloneq\sum_{1\leq i,j\leq d}\tau_{ij}\eta_{ij}$ with corresponding norm such that, for all $\matr{\tau}\in\mathbb{R}^{d\times d}$, $|\matr{\tau}|_{d\times d}\coloneq(\matr{\tau}:\matr{\tau})^{\nicefrac{{1}}{{2}}}$.

3.2 Well-posedness

The aim of this section is to infer a stability estimate on the displacement and pressure $(\vec{u},p)\in L^{2}((0,t_{\rm F});\vec{U})\times L^{2}((0,t_{\rm F});P)$ solving (17) yielding, in particular, the well-posedness of the weak problem. The existence and uniqueness of a solution to the deterministic Biot problem has been studied in [58, 63]. The results therein establish the existence, for almost every (a.e.) $\theta\in\Theta$, of a solution $(\vec{u}(\theta),p(\theta))$ to (16). Following the argument of [15, Theorem 2.2], we note that the mapping $(\vec{u},p):\Theta\to L^{2}((0,t_{\rm F});\textbf{H}_{0,\Gamma_{D}}^{1}(D))\times L^{2}((0,t_{\rm F});H^{1}(D))$ is measurable because it is a continuous function of the input coefficients $\mu,\lambda,\alpha,c_{0},\kappa\in L^{2}_{\mathcal{P}}(\Theta)$. Therefore, proving that $\theta\mapsto(\vec{u}(\theta),p(\theta))$ is bounded in $L^{2}_{\mathcal{P}}(\Theta)$ gives the existence of a solution to (17). This result is established by Proposition 5.

We assume some additional conditions on the input random data that will be needed in the proofs. First, we recall that, for all $\vec{x}\in D$ and all $\theta\in\Theta$, the coupling coefficient $\alpha(\vec{x},\theta)\in(0,1]$ satisfies the lower bound (15). We assume that the reference porosity of the medium is strictly positive (otherwise the Biot problem would reduce to a decoupled one), so that we have $0<\underline{\alpha}\coloneq\frac{3\varphi}{2+\varphi}<\alpha(\vec{x},\theta)$. We also assume that Lamé’s coefficients satisfy the following:

Assumption 3** (Elastic moduli).**

The shear modulus $\mu\in L^{\infty}(D\times\Theta)$ is uniformly bounded away from zero, i.e. there exists a positive constant $\underline{\mu}$ such that

[TABLE]

The bulk modulus $K\in L^{\infty}(D\times\Theta)$ is uniformly bounded from above, i.e. there exists a constant $\overline{K}$ such that

[TABLE]

The next result establishes the coercivity and the inf-sup condition of the bilinear forms $a$ and $b$ in (17a).

Lemma 4** (Coercivity and inf-sup).**

Under Assumption 3, the following bounds hold:

[TABLE]

where $C_{\rm K}>0$ denotes the constant in Korn’s first inequality, $C_{\rm is}>0$ the inf-sup constant, and $\|\cdot\|_{\vec{U}}\coloneq\|\cdot\|_{L^{2}_{\mathcal{P}}(\Theta;{\textbf{H}_{0,\Gamma_{D}}^{1}(D)})}$. In the case $\Gamma_{D}=\partial D$, (21) holds for all $q\in L^{2}_{\mathcal{P}}(\Theta;L_{0}^{2}(D))$.

Proof.

The first bound directly follows from Assumption 3. Indeed, since $|\Gamma_{D}|_{d-1}>0$, we can apply Korn’s first inequality, yielding

[TABLE]

To obtain (21), we use [5, Lemma 7.2] establishing the existence, for any $q\in L^{2}_{\mathcal{P}}(\Theta;L^{2}(D))$ (or $q\in L^{2}_{\mathcal{P}}(\Theta;L_{0}^{2}(D))$ in the case $\Gamma_{d}=\partial D$), of ${\vec{v}}_{q}\in\vec{U}$ such that $\vec{\nabla}{\cdot}{\vec{v}}_{q}=q$ and $C_{\rm is}\|{\vec{v}}_{q}\|_{\vec{U}}\leq\|q\|_{L^{2}_{\mathcal{P}}(\Theta;L^{2}(D))}$, with $C_{\rm is}>0$ depending on $D$. Thus, we conclude

[TABLE]

∎

In order to prove the stability estimate of Proposition 5 below, we infer from (19) a lower bound for the specific storage coefficient $c_{0}$. We rewrite (13), derived from the definition of $\alpha$ in (4) and the Gassmann equation (12), in the form

[TABLE]

Since the left-hand side of the previous relation is always non-negative (according to [19, Chapter 4] we have $K_{\rm d}\leq K\leq K_{\rm m}$), we infer that $\nicefrac{{\alpha}}{{c_{0}}}\leq K$ and, as a result,

[TABLE]

Proposition 5** (A priori estimate).**

Let $(\vec{u},p)\in L^{2}((0,t_{\rm F});\vec{U}\times P)$ solve (17). Then, under Assumption 3, it holds

[TABLE]

with $G:(0,t_{\rm F})\to L^{2}(D)$ defined by

[TABLE]

Proof.

Since the loading $\vec{f}$ is deterministic, from (17a) we infer that $a(\vec{u}(t),\vec{v})+b(\vec{v},p(t))=\left(\vec{f}(t),\mathbb{E}(\vec{v})\right)_{D}$, for ${\rm a.e.}\;t\in(0,t_{\rm F})$ and all $\vec{v}\in\vec{U}$. Setting $\vec{v}(\vec{x},\theta)=\vec{u}(\vec{x},\theta,t)$ and integrating the previous relation on $(0,t_{\rm F})$, yields

[TABLE]

Adapting the argument of [9, Chapter 4] and [11, Remark 2], we infer that $b(\vec{u}(t),q)-c(p(t),q)\in H^{1}((0,t_{\rm F}))$ for all $q\in P$. Thus, we can integrate by parts (17b) and obtain

[TABLE]

Letting $t\in(0,t_{\rm F})$, integrating the previous identity on $(0,t)$, and then taking $q(\vec{x},\theta)=p(\vec{x},\theta,t)$, leads to

[TABLE]

We define $z(t)=\int_{0}^{t}p(s){\rm d}{s}$ and observe that ${\rm d}_{t}z(t)=p(t)$ and $z(0)=0$. Using the linearity of $d$ and the formula ${\rm d}_{t}z^{2}(t)=2z(t){\rm d}_{t}z(t)$ to rewrite the third term in the left hand side of the previous relation, and recalling (24), we get

[TABLE]

Then, integrating on $(0,t_{\rm F})$ and summing the resulting identity to (25), gives

[TABLE]

To establish the result, we now bound the right-hand side of the previous relation. First, we observe that, owing to the Jensen inequality, for all $\vec{v}\in\vec{U}$ it holds

[TABLE]

Using the previous relation and the Cauchy–Schwarz and Young inequalities followed by (20), it is inferred that

[TABLE]

Proceeding similarly for the second term in the right-hand side of (26) and recalling the lower bound (22), one has

[TABLE]

Plugging the two previous estimates into (26), multiplying by a factor $2$, and using $d(z(t_{\rm F}),z(t_{\rm F}))\geq 0$, yields the conclusion. ∎

Some remarks are in order.

Remark 6* (Inf-sup condition).*

We observe that it is possible to derive the previous a priori estimate without (22). Indeed, the inf-sup condition (21) together with (17a) and (19) allow to bound the second term in the right-hand side of (26) using the following: for a.e. $t\in(0,t_{\rm F})$,

[TABLE]

The resulting stability estimate would have, compared to (23), an additional dependence on $\underline{\alpha}C_{\rm is}$.

Remark 7* (Data regularity).*

In order to prove the a priori bound (23), no additional time regularity assumption on the loading term $\vec{f}$ and the mass source $g$ are needed. However, under the additional requirements $\vec{f}\in C^{1}([0,t_{\rm F}];\textbf{L}^{2}(\Omega))$ and $g\in C^{1}([0,t_{\rm F}];L^{2}(\Omega))$, a stronger version of (23) can be inferred using the argument of [48] and [54].

Remark 8* (Quasi-incompressible media).*

In order to prove the stability estimate (23), no additional assumption on the mobility $\kappa:D\times\Theta\to\mathbb{R}^{+}$ is required. Thus, Proposition 5 can handle the case of locally poorly permeable media (i.e. $\kappa\ll 1$). However, assuming (19) does not allow to consider quasi-incompressible materials for which $\lambda\gg 1$. To obtain a robust estimate in the case of $\lambda$ unbounded, we can proceed as in [43, Theorem 1]. Assuming $\vec{f}\in H^{1}((0,t_{\rm F});\textbf{L}^{2}(D))$ and $\kappa$ uniformly bounded away from zero, the Darcy term gives a $L^{2}((0,t_{\rm F});P)$ estimate of the pressure and, as a consequence, (19) and (22) are not needed. We remark that, in a medium featuring very low permeability, an incompressible fluid cannot flow unless the material is compressible, namely the two limit cases $\kappa\ll 1$ and $\lambda\gg 1$ cannot occur simultaneously. Therefore, the stochastic poroelasticity problem need to be dealt with uncertainty models preventing this situation.

Remark 9* (Lamé’s coefficients).*

The well-posedness of problem (17) holds under weaker assumption on $\mu$ and $\lambda$. More precisely, one can assume instead of (18) and (19), that $\mathcal{P}$-a.e. in $\Theta$ it holds

[TABLE]

where $\underline{\mu}(\theta)^{-1}$ and $\overline{K}(\theta)$ are second-order random variables. An assumption of this type is convenient when $\lambda$ is an unbounded (e.g. Gaussian or lognormal) random variable at $\vec{x}\in D$. See [2, Lemma 1.2] and [14] for a discussion in the context of elliptic PDEs with random data.

4 Probabilistic framework

In many applications, only limited information about the poroelastic coefficients in (7) and (11) is available. In the context of geomechanics, e.g., even when it is possible to carry out a large number of measurements, the actual knowledge of the soil properties typically suffers from inaccuracies: due to the presence of different layers, the physical properties can have strong variations that are difficult to estimate.

4.1 Uncertain poroelastic coefficients

In this section we present the probabilistic setting and we introduce the parametrization of the random coefficients using canonical random variables. Some preliminary investigations on the possible choices to model the uncertain coefficients are subsequently presented.

4.1.1 Parametrization with canonical variables

We use a set of canonical random variables, collected into a random vector $\vec{\xi}=(\xi_{1},\ldots,\xi_{N}):\Theta\to\Xi$, to describe the uncertainty of the poromechanical coefficients. To ensure the independence of the random coefficients, it suffices to use one canonical random variable per coefficient, and rely on changes of measure to define the mappings between the $\xi_{i}$ and their respective coefficients. Specifically, we deal with uniformly iid canonical variables, namely $\xi_{i}\sim\mathcal{U}([-1,1])$ and $\xi_{i}\perp\xi_{j}$ if $i\neq j$. Thus, the $\xi_{i}$’s have zero mean, are independent with product form probability law $\mathcal{P}_{\vec{\xi}}$ and joint density function $\rho(\vec{\xi})$. We denote by $\mathcal{B}_{\vec{\xi}}$ the Borel $\sigma$-algebra on $\Xi$ and by $(\Xi,\mathcal{B}_{\vec{\xi}},\mathcal{P}_{\vec{\xi}})$ the image probability space. We also define the space of second-order random variables on $(\Xi,\mathcal{B}_{\vec{\xi}},\mathcal{P}_{\vec{\xi}})$ as the weighted Lebesgue space $L^{2}_{\rho}(\Xi)$. The expectation operator on the image space is denoted using brackets and is related to the expectation on $(\Theta,\mathcal{B},\mathcal{P})$ through the identity

[TABLE]

The variance operator of $h\in(\Xi,\mathcal{B}_{\vec{\xi}},\mathcal{P}_{\vec{\xi}})$ is then defined by

[TABLE]

4.1.2 Probabilistic model

In order to account for the uncertainty of the poroelastic material, a model for its properties is needed. For the sake of simplicity, we assume that the fluid bulk modulus and the porosity are deterministic: $K_{\rm f}=2.2\ {\rm GPa}$ and $\varphi=0.2$. Our input uncertainty is then parametrized by a vector $\vec{\xi}$ of dimension $N=4$. The set of random coefficients consists in the two Lamé’s constants $\lambda(\vec{\xi})$ and $\mu(\vec{\xi})$, the Biot–Willis coefficient $\alpha(\vec{\xi})$, and the hydraulic mobility $\kappa(\vec{\xi})$. These coefficients will be considered all mutually independent, with no spatial variabilities. The case of spatial variability can be handled by introducing a Karhunen–Loève expansion of the corresponding stochastic fields [47, 57] but involving a high number of input random variables. Even if the experiments presented in this work are mainly academic model problems, we try to consider realistic sets of parameters and have a particular regard for the underlying physical phenomena. The ranges of variation of the uncertain input data are inspired from [56, Table 1], [67, Section 3], and [19, Table 4.1], and correspond to an ideal water filled sand soil with a $20\%$ reference porosity. In particular, for the numerical investigations of Section 5.1, we assume that $\mu,\lambda$, and $\kappa$ have a log-uniform distribution, while $\alpha$ is uniformly distributed. We set

[TABLE]

According to (4) and (15), we take $\alpha_{\rm min}=3\varphi(2+\varphi)^{-1}$ and $\alpha_{\rm max}=1$. This choice yields

[TABLE]

where ${\rm c_{v}}(\cdot)\coloneq\frac{\mathbb{V}{\rm ar}(\cdot)^{\nicefrac{{1}}{{2}}}}{\langle\cdot\rangle}$ denotes the coefficient of variation operator. In particular we observe that the previous average values have the same magnitude as the coefficients considered in [29, Table 2]. The case of other standard distributions can be handled similarly introducing the inverse of the cumulative distribution function of the considered coefficient.

According to the results presented in Section 2.3, the constrained specific storage coefficient $c_{0}$ can be expressed in terms of the elastic moduli $\lambda,\mu$ and the coupling Biot coefficient $\alpha$ (see also [18, 36] for theoretical and empirical investigations on the storage coefficient) as well as the porosity and the fluid bulk modulus. Therefore, we let

[TABLE]

The advantage of this approach is threefold:

(i) it produces a strategy of perturbation with independent model coefficients,

(ii) it ensures that the set of poroelastic coefficients belong to the physical admissible set, and

(iii) it allows to reduce the uncertainty dimension and then the size of the sampling (i.e. the number of sparse grid points).

We illustrate here the definition of $c_{0}(\vec{\xi})$ in (29), reporting its probability distribution for random poromechanical coefficients as in (27). For each realization of the random poromechanical coefficients, we solve (14) in order to compute the corresponding value of $c_{0}$.

In Figure 1 we plot the histogram and probability density function of $c_{0}$, estimated using $10^{5}$ realizations of the poromechanical coefficients. We observe that our choice leads to a specific storage coefficient $c_{0}$ ranging in $(2\cdot 10^{-4},3\cdot 10^{-1})\ {\rm kPa}^{-1}$ and its probability density function is far from being uniform with a highest probable value close to zero.

Different choices are possible to describe the mechanical behavior of the material and, in practice, the selection of the particular coefficients defining the elastic properties of the medium is mostly one of convenience. For instance, one might consider the Young modulus $E$ (e.g. in [42]) or the bulk modulus $K$ and the Poisson ratio $\nu$ (e.g. in [56, 18]), as the independent driving coefficients. Here we choose instead Lamé’s coefficients $\mu$ and $\lambda$ as primary independent variables, because they explicitly appear in the governing equations and are also more conveniently defined through simple transformations of independent canonical variables. The corresponding couple $(K,\nu)$ can be deduced from Lamé’s coefficients $(\mu,\lambda)$ using the relations

[TABLE]

Figure 2 shows the marginals of $K$ and $\nu$ for independent primary mechanical variables $(\mu,\lambda)$ distributed according to (27).

Starting from independent $K$ and $\nu$ results in a complex dependent distribution of the coefficients $\lambda$ and $\mu$ appearing in the model equations. This is due to the non-linear relation $(K,\nu)\mapsto(\lambda,\mu)$ that can be appreciated from Figure 3 for the case where $K$ is log-uniformly distributed and $\nu$ is uniformly distributed.

4.2 Stochastic discretization

Now that the uncertain coefficients and their parametrization have been detailed, we introduce in this section the Polynomial Chaos expansions of the random model solution and outline the Pseudo-Spectral Projection algorithm used for its determination. The PSP method is one of the so-called non-intrusive techniques, which rely on an ensemble of deterministic model simulations to estimate the expansion coefficients of the solution.

4.2.1 Polynomial Chaos expansion

PC expansions were initially proposed by Wiener [65] for the approximation of second order random quantities from standard iid Gaussian variables. It was later proposed to use these polynomial expansions for UQ [35]. Concerning the generalization to non-Gaussian measure and non-polynomial expansions, see [66] and references within [45].

Let us denote by $\{\phi_{\vec{k}}(\vec{\xi}):\vec{k}\in\mathbb{N}^{N}\}$ an Hilbertian basis of $L^{2}_{\rho}(\Xi)$, where $\phi_{\vec{k}}$ is a multivariate polynomial in $\vec{\xi}$ and $\vec{k}=(k_{1},\ldots,k_{N})$ is a multi-index indicating the polynomial degree in the $\xi_{i}$’s. The total degree of $\phi_{\vec{k}}$ is $|\vec{k}|_{1}\coloneq\sum_{i=1}^{N}k_{i}$. The basis functions are commonly chosen to be orthogonal with respect to the inner product in $L^{2}_{\rho}(\Xi)$ characterized by the probability density function $\rho:\Xi\to\mathbb{R}^{+}$,

[TABLE]

In our case, where $\Xi$ is the hypercube $[-1,1]^{N}$ with uniform density, the $\phi_{\vec{k}}$ are in fact the multivariate Legendre polynomials. If $X(\vec{\xi})\in L^{2}_{\rho}(\Xi)$ is a second-order random variable, then it admits the PC expansion [65, 12],

[TABLE]

where the deterministic coefficients of the series $\{X_{\vec{k}}:\vec{k}\in\mathbb{N}^{N}\}$ are called the spectral modes. The PC approximation $X_{\mathcal{K}}(\vec{\xi})$ of $X(\vec{\xi})$ is obtained by truncating the expansion above to a finite series,

[TABLE]

where $\mathcal{K}\subset\mathbb{N}^{N}$ is a finite set of multi-indices. Observing that the degree zero polynomial is $\phi_{\vec{0}}=1$, the expectation and variance of $X(\vec{\xi})$ simply express as

[TABLE]

Different truncation strategies, i.e. the selection of the set $\mathcal{K}$, can be considered, with different resulting error $e_{\mathcal{K}}(\vec{\xi}):=X(\vec{\xi})-X_{\mathcal{K}}(\vec{\xi})$. In particular, if the approximation is defined as the $L^{2}$-orthogonal projection of $X$, the spectral coefficients are defined by

[TABLE]

and it is clear that the $L^{2}$-error norm does not increase when $\mathcal{K}$ is enlarged. Indeed, we have in this case

[TABLE]

A classical truncation strategy is based on the degree of the expansion; for instance, one can define $\mathcal{K}$ such that the approximation basis contains all multivariate polynomials with total ($|\vec{k}|_{1}\leq p$) or partial ($|\vec{k}|_{\infty}\coloneq\max_{1\leq i\leq N}k_{i}\leq p$) degree up to a prescribed value $p\geq 0$. The convergence rate of the PC expansion of $X$ with respect to the degree $p$ depends on the regularity (in $\vec{\xi}$) of $X$. In particular, an exponential convergence is expected for analytical variables; more details about the convergence conditions can be found in [12, 27].

4.2.2 Sparse Pseudo-Spectral Projection

Several methods have been proposed to compute the spectral modes of the PC expansion [45, 40]. Among these methods, we distinguish between the intrusive and non-intrusive ones. The former involve a reformulation of the problem to derive a set of governing equations for the solution modes, while the latter rely on the availability of a deterministic solver only. As a first attempt to apply PC expansions to propagate uncertainty in the Biot model, we opted for a non-intrusive projection method which, in our experience, presents the best trade-off between numerical complexity and precision, compared to alternatives such as regression and compressed sensing.

The Non-Intrusive Spectral Projection (NISP) method is based on approximating the correlation in the right-hand-side of (34) by a deterministic numerical quadrature formula, such as

[TABLE]

where the $N_{q}$ quadrature nodes $\vec{\xi}^{(q)}$ and weights $w^{(q)}$ are classically constructed by tensorization of one-dimensional quadrature rules (e.g. Gauss, Féjèr, Clenchaw–Curtis, …). The complexity of the NISP method is governed by the number $N_{q}$ of evaluations of the model, while the accuracy depends on the basis which can be considered given the quadrature rule. Clearly, these two aspects are intertwined, as larger expansion bases require higher order quadrature rules with larger computational efforts. In practice, the basis is defined as the largest one such that the discrete projection is exact for any function belonging to the linear-span of the basis or, in other words, that the non-intrusive projection incurs no internal aliasing. Equivalently, the multi-index set $\mathcal{K}$ must satisfy

[TABLE]

Exploiting the product structure of the stochastic space, sparse tensorizations of one-dimensional quadrature rules [31, 32] have been proposed to mitigate computational complexity, especially for a larger number $N$ of canonical random variables. The sparse constructions generally rely on the Smolyak’s formula [59] and lead to the so-called sparse-grid projection methods. One drawback of using sparse quadrature rules, as opposed to fully tensorized quadrature rules, is that some of the weights $w^{(q)}$ are negative even for one-dimensional formulas with positive weights only. Therefore, the sparse quadratures are not constituting discrete inner products, with some restrictions on the possible aliasing-free bases as a result.

These limitations were recognized in [16, 17], where the authors proposed the so-called Pseudo-Spectral Projection (PSP) methods. The key-idea of PSP is to apply the Smolyak’s formula directly on the projection operator, rather than on the integration operator. The projection is seen as a sequence of nested subspace-projections, where each subspace-projection is computed using a fully-tenzorized quadrature rule. Consequently, the determination of each spectral mode uses an adapted fully-tensorized quadrature rule with improved properties compared to the unique one in (35). Specifically, for the same sparse grid, with $N_{q}$ nodes $\vec{k}^{(q)}$, the PSP method can compute without internal aliasing the projection coefficients $X_{\vec{k}}$ of a larger set $\mathcal{K}$ of basis functions $\phi_{\vec{k}}$, compared to the NISP approach.

In this work, we then use a sparse PSP method, constructed on nested one-dimensional Clenshaw–Curtis quadrature rules. As the canonical random variables account for uncertainties in different poromechanical coefficients whose respective influences on the model solution are unknown, we decided to rely on an isotropic construction. In the PSP context, it amounts to consider isotropic partial tensorization, with a unique parameter controlling the complexity and accuracy of the approximation: the maximum level $l\in\mathbb{N}$ of isotropic the sparse grid. As the level $l$ of the PSP method increases, both the number of sparse grid nodes $N_{q}(l)$ and the multi-index set $\mathcal{K}(l)$ increase.

5 Test cases

We now assess the performance of the PSP method detailed in the previous section to treat the stochastic part of the Biot problem (16) due to parametric uncertainties. Recall that the PSP method requires to solve $N_{q}$ deterministic Biot problems (discretized by the method of [8]) for different values of uncertain input parameters. Note that the PSP method can be applied with alternative solvers, possibly based on a non-monolithic formulation (i.e., splitting methods). Our numerical experiments feature a simplified model with spatially homogeneous random coefficients. More sophisticated uncertain models, requiring (possibly strongly correlated) spatially varying input coefficients with non-standard distributions will be the object of future works. We begin this section with two validation cases; the first one has a fluid forcing term while a mechanical boundary condition drives the second one. An injection-production test case at large spatial scale is then presented to assess the methodology on a more realistic problem and perform a sensitivity analysis to illustrate the interest of the method.

5.1 Validation test cases

To start, we implement two experiments to study the convergence of the PSP method on the Biot problem. These two test cases are complementary in the sense that they consider a fluid or mechanical solicitation of the porous medium. More precisely, the first case is an injection test designed to investigate the effect of a fluid source on the deformation of the domain, while the second case is meant to evaluate the mechanical effect of a traction term on the fluid pressure. These two cases are defined on the unit square domain $D=[0,1]^{2}$ and the loading term $\vec{f}$ in (16a) as well as the initial condition $\phi_{0}$ in (16g) are both equal to zero. The probability distributions of the random poroelastic coefficients are defined in Section 4.1. After describing the test cases, we numerically investigate the convergence of the stochastic error with respect to the level of the sparse grid. We also analyze the variance and covariance fields of the PC solutions.

5.1.1 Point injection

The first experiment is inspired by the Barry and Mercer test case [8, Section 6.2] (see also  [3, 55]) in which the displacement and pressure fields are driven by a stationary fluid source

[TABLE]

where $\vec{x}_{0}=(0.25,0.25)$ denotes the source location and $\delta$ the Dirac delta function. Homogeneous mechanical and fluid boundary conditions are enforced on the whole boundary of the domain,

[TABLE]

where $\vec{\tau}$ denotes the tangent unit vector on $\partial D$ such that $\vec{\tau}\perp\vec{n}$ with arbitrary orientation. The following results have been obtained using the method of [8] with polynomial degree $k=3$ on a Cartesian mesh containing $1{,}024$ elements. Using the static condensation procedure [8, Section 5], the discrete system has $27{,}648$ unknowns. We are interested in the stationary solution, namely the pressure and displacement fields manifesting after the initial transient phase, and report the solution at the final time $t_{\rm F}=1{\rm s}$ obtained after $10$ time steps of the implicit Euler scheme.

Figure 4 represents the two first statistical moments (33) of the pressure field approximated using a level $l=5$ sparse grid. We plot the solution on the deformed configuration obtained by applying the (scaled) computed mean displacement field to the reference unit (undeformed) square domain $D$. As expected, we observe a dilatation of the domain caused by the injection. The computed solutions are symmetric with respect to the diagonal according to the source location and the homogeneous boundary conditions and to the fact that the computational mesh also has the same symmetry. The mean pressure field exhibits maximum values at the source, while only small perturbations are observed near the boundary. The variance of the pressure field has the same behavior as the mean, with higher values close to the injection point. We notice that the magnitude of the standard deviation, i.e. the square root of the variance, is comparable to the mean field range. We draw attention to the fact that the variance (and covariance) fields must be carefully computed since the spatial dependence involves the products of broken polynomial functions on the mesh and its skeleton.

5.1.2 Poroelastic footing

As a second example, we examine the 2D footing problem proposed in [52, 29], where the domain is assumed to be free to drain and fixed along the base and vertical edges. A uniform load is applied on a central portion of the upper boundary $\Gamma_{N}=\{\vec{x}=(x_{1},x_{2})\in\partial D:x_{2}=1\}$ simulating a footing step compressing the medium. The boundary conditions are defined as

[TABLE]

As $\vec{f}=\vec{0}$ and $\phi_{0}=0$, the displacement and pressure fields are determined only by the stress boundary condition. We set $k=2$ in the method of [8] and discretize the domain using a triangular mesh with $3{,}584$ elements. We focus here on the pressure profile at early times and, for each sparse grid node, the linear system of dimension $54{,}720$ is solved twice to reach $t_{\rm F}=0.2{\rm s}$. Because of the free drainage boundary condition, the fluid is promptly squeezed out of the soil to reach an equilibrium configuration with no pressure and where the displacement balances the Neumann condition. We mention, in passing, that the spurious pressure oscillations observed in [52] do not occur with our discretization in agreement with results in [8, Section 6.2].

Figure 5 shows the mean and variance of the pressure field over the deformed domain. The deformation of the domain corresponds again to the mean displacement field. The domain is particularly warped on the load segment and the maximum pressure is found near the center of the domain. As in the previous test case, we observe that the magnitude of the standard deviation is comparable to the mean field range.

5.1.3 Convergence analysis

We evaluate the accuracy of the PC expansions of the stochastic solution by using a validation set independent from the sparse grid nodes. The Mean-Squared Error (MSE) is computed from $N^{*}=500$ Latin Hypercube Samples (LHS), denoted $\vec{\xi}^{(i)}$,

[TABLE]

where $X(\vec{x},\cdot)$ denotes the exact (pressure or displacement) field and $X_{\mathcal{K}}(\vec{x},\cdot)$ its PC approximation. We plot on Figures 7 and 8 the MSE fields of the PC solutions. Owing to the symmetry of each configuration, we only represent the error of the first component of the displacement field for the injection test and the second displacement component for the footing test. The results are plotted for the sparse grids of level $l=1,3,5$, and we observe a rapid decrease of the error for the two test cases and both the displacement and pressure fields. We especially notice the strong decay of the maximum error value with the increasing sparse grid level. To complete the convergence analysis, we also plot inn Figure 6 the $L^{2}(D)$-norm of the mean-squared error fields, $\|{\rm MSE}(\vec{u}-{\vec{u}}_{\mathcal{K}})\|_{\textbf{L}^{2}(D)}$ and $\|{\rm MSE}(p-{p}_{\mathcal{K}})\|_{L^{2}(D)}$ up to a sparse grid level equal to 5. The linear trends observed in the semi-log scale of the successive PC approximations suggest that the method achieves roughly an exponential convergence rate with respect to the level of the sparse grid, as predicted by the theory for a sufficiently regular quantity of interest.

5.1.4 Covariance fields

We finally analyze the correlation between the displacement and pressure fields. The covariance between two random variables $h_{1},h_{2}\in(\Xi,\mathcal{B}_{\vec{\xi}},\mathcal{P}_{\vec{\xi}})$ is defined by

[TABLE]

In the case of two random variables represeneted on by PC expansions in the same basis, $X_{\mathcal{K}}$ and $Y_{\mathcal{K}}$, the covariance reads

[TABLE]

On Figure 9 we plot the covariance fields of the injection test case. As expected, we observe a symmetry with the domain diagonal of the covariance field between the two displacement components (left plot). The sign of the covariance depends on the joint variation of the variables. Indeed the covariance field of ${u}_{\mathcal{K},1}$ and ${u}_{\mathcal{K},2}$ is in accordance with the fact that both $u_{1}$ and $u_{2}$ are negative on the bottom-left area, both $u_{1}$ and $u_{2}$ are positive on the top-right area, while either $u_{1}$ or $u_{2}$ is positive when the other component is negative on the top-left and bottom-right areas of the domain. High values of $\mathbb{C}{\rm ov}({u}_{\mathcal{K},1},p_{\mathcal{K}})$ and $\mathbb{C}{\rm ov}({u}_{\mathcal{K},2},p_{\mathcal{K}})$ (middle and right panel) are found close to the injection point, where the pressure and mechanical effects of the fluid source are maximum. We also remark that, owing to the symmetric configuration of the injection test, the covariance field $\mathbb{C}{\rm ov}({u}_{\mathcal{K},1},p_{\mathcal{K}})$ corresponds roughly to the $90^{\circ}$ clockwise rotation of $\mathbb{C}{\rm ov}({u}_{\mathcal{K},2},p_{\mathcal{K}})$.

Figure 10 represents the covariance fields for the footing test case. The covariance field of ${u}_{\mathcal{K},2}$ and $p_{\mathcal{K}}$ (right panel) has the same pattern as the pressure variance field plotted on Figure 5 with the highest value below the loading segment. The covariance is here negative since the pressure increases during the vertical compression. On the middle panel, the covariance of the horizontal displacement ${u}_{\mathcal{K},1}$ and the pressure $p_{\mathcal{K}}$ is lower by one order of magnitude consistently with the vertical loading of this test case. The sign and the location of maximum and minimum values of $\mathbb{C}{\rm ov}({u}_{\mathcal{K},1},{u}_{\mathcal{K},2})$ and $\mathbb{C}{\rm ov}({u}_{\mathcal{K},1},p_{\mathcal{K}})$ (left and middle panel) are consistent with the deformed configuration plotted on Figure 5.

5.2 Realistic problem

The last numerical experiment deals with a fluid injection-extraction process and is inspired by [43, Section 7]. Since we aim to assess the propagation of the uncertainty of input parameters on the solution, we assume that all the coefficients $\mu,\lambda,\alpha$ and $\kappa$ are log-uniformly distributed with the same coefficient of variation ${\rm c_{v}}(\mu)={\rm c_{v}}(\lambda)={\rm c_{v}}(\alpha)={\rm c_{v}}(\kappa)\sim 0.63$. We set

[TABLE]

We remark that all the possible realizations of $\alpha$ satisfy (15) if $\varphi\leq\nicefrac{{2}}{{29}}$. Thus, we consider a medium with reference porosity $\varphi=\nicefrac{{2}}{{29}}$ that is filled with water, namely $K_{\rm f}=2.2\ {\rm GPa}$. For each realization of the poromechanical coefficients, $c_{0}$ is again computed according to (14). We observe that the average values of $\mu$ and $\lambda$ in (36) have the same magnitude as the mechanical coefficients considered in [43, Table 1] and, assuming that $\mu_{\rm f}=10^{-3}\ {\rm Pa}\cdot\rm{s}$ (water), the average value of $\kappa$ corresponds to an intrinsic permeability of order $1{\rm D}=10^{-12}{\rm m^{2}}$.

5.2.1 Computational case description

The rectangular domain $D=[0,4\,{\rm Km}]\times[0,1\,{\rm Km}]$ includes two holes representing the injection and extraction wells (see Figure 11). For the sake of simplicity, we impose steady Dirichlet conditions on the holes boundaries,

[TABLE]

where $\Gamma_{{\rm d},1}$ and $\Gamma_{{\rm d},2}$ denote the left and right hole boundaries, respectively. We remark that, due to the linearity of the Biot model, any solution with the conditions above can be reconstructed from a combination of the solutions of two elementary problems: one with $p_{1}=0,p_{2}=1$ and another with $p_{1}=1,p_{2}=0$. We treat in this section the case $p_{2}=-p_{1}$ and consider the following set of boundary conditions:

[TABLE]

where $\vec{\tau}$ denotes again the tangent vector on $\partial D$ with arbitrary orientation. The loading term $\vec{f}$ in (16a), the fluid source $g$ in (16b), and the initial condition $\phi_{0}$ in (16g) are all set to zero. We use $k=1$ in the method of [8] with the Voronoi mesh having $10{,}000$ elements depicted in Figure 11. This choice yields a spatial discretization of dimension $151{,}240$. The steady state of the Biot model is considered achieved at $t_{\rm F}=1\ {\rm day}$, simulated with $10$ time steps. The PC approximations are computed using the PSP method with level $l=3$ for a total of $N_{q}(l)=209$ sparse grid nodes.

Figure 12 shows the mean of the pressure (left) and vertical displacement fields (right). We verify that the mean pressure is equal to deterministic the boundary conditions, namely $\pm 100\,{\rm kPa}$, at the holes boundary. It increases linearly between the two holes (in the $x$-direction) and remains nearly uniform away from the holes. Regarding the displacement, within vertical sections, the vertical component increases from the bottom of the domain (where it is equal to the homogeneous boundary condition) to its maximum value at the free surface. The mean displacement is negative on the left part of the domain (extraction side) and positive on the right part (injection side), consistently with the imposed pressure boundary conditions. Further, the mean displacement vanishes at the vertical centerline of the domain, owing to the symmetry of the problem.

Figure 13 shows the variance fields of the Biot solution, on the left column, and the three covariance fields, in the right column. The pressure and vertical displacement variances are small in between the two holes. Indeed, for the considered configuration, the pressure profile at steady state is linear in the $x$-direction and mostly independent of the model coefficients. Thus the pressure variance is maximal at the lateral boundaries where the uncertainty in poroelastic coefficients is the most significant. In contrast, the horizontal displacement variance is the highest between the holes (note that the range of $\mathbb{V}{\rm ar}({u}_{\mathcal{K},1})$ is half that of $\mathbb{V}{\rm ar}({u}_{\mathcal{K},2})$).

The three covariance fields display different features. Concerning the correlations $\mathbb{C}{\rm ov}({u}_{\mathcal{K},1},p_{\mathcal{K}})$ and $\mathbb{C}{\rm ov}({u}_{\mathcal{K},2},p_{\mathcal{K}})$, it is seen that the later reaches higher magnitudes at the top corners of the domain, as one could have guessed from the variance fields. Further, the pressure and vertical displacement are anti-correlated everywhere in the domain, when the horizontal displacement and pressure have a covariance that alternates sign in the domain. Finally, consistently with the symmetry of the problem, the covariance between the two displacement components is negative (resp. positive) in the injection (resp. extraction) side of the domain.

5.2.2 Sensitivity analysis

We conclude the numerical tests by evaluating the different contributions of the input parameters on the variance of the solution. We are mainly interested here in illustrating the information carried by the PC expansion to assess the effect of the random coefficients on the deformation caused by the injection-extraction. We perform the sensitivity analysis of the vertical displacement only. Similar analyses can be conducted for the other components of the solution. We compute the first-order partial variances $\mathbb{V}{\rm ar}_{i}(u_{\mathcal{K},2}(\vec{x},\vec{\xi}))$ defined by

[TABLE]

As for the variance, these partial variances can be explicitly derived from the PC expansion of the considered field [20], and we use here the level $l=3$ sparse grid PSP approximation. Since each $\xi_{i}$ in (36) parametrizes only one of the poromechanical model coefficients, the first-order variances characterize the contributions to the variance of the individual poro-mechanical parameters ($\mu(\xi_{1})$, $\lambda(\xi_{2})$, $\alpha(\xi_{3})$ and $\kappa(\xi_{4})$) on the displacement. We also report the total-order variance associated to the $\xi_{i}$, denoted $\mathbb{V}{\rm ar}_{T_{i}}(u_{\mathcal{K},2})$, which are defined by

[TABLE]

where we have denoted $\vec{\xi}_{\setminus i}=\left(\xi_{j\neq i}\right)$ the triplet containing all the random variables but $\xi_{i}$. The total-order variance $\mathbb{V}{\rm ar}_{T_{i}}$ measure the variability induced by $\xi_{i}$ (that is the poromechanical coefficient it parametrizes) including all its interactions with other uncertainty sources.

In Figure 14 we compare the total variance $\mathbb{V}{\rm ar}(u_{\mathcal{K},2})$ with the sum of the first-order partial variances $\sum_{i=1}^{4}\mathbb{V}{\rm ar}_{i}(u_{\mathcal{K},2})$ of the vertical displacement. The difference between these two fields amounts to the contribution of higher order partial variances, namely, the part of the variance incurring to interactions between random coefficients. Here, these interactions are seen to be small but not negligible. We expect coupled effects between uncertain parameters in the present uncertainty model, because of the dependence of $c_{0}$ on the three first coefficients (see (14)). In Figure 15 we plot the first-order (left half) and total-order partial variances (right half) related to each input parameter. We observe that the parameter with the main influence on the vertical displacement is the shear modulus $\mu$ (parametrized by $\xi_{1}$), while the effects of the dilatation modulus $\lambda$ (parametrized by $\xi_{2}$) and hydraulic mobility $\kappa$ (parametrized by $\xi_{4}$) are much weaker. In particular, we notice that $\kappa$ has almost no influence on the solution variability. Finally, we remark that the effect of the uncertain Biot-Willis coefficient $\alpha$ on the vertical displacement at the surface is less important than the one associated to $\mu$, but considerably larger than the one related to $\lambda$.

6 Conclusion

This paper contains a theoretical and numerical study of the poroelasticity problem with random coefficients. It contains three major contributions. The first one is a study of the dependence between the physical coefficients, leading to an uncertainty model where the specific storage coefficient is expressed in terms of the other coefficients. This construction ensures that the resulting problem remains almost surely within the set of physically possible configurations. This uncertainty model also results in a reduction in the number of independent coefficients for the stochastic problem to be solved. The second contribution concerns the well-posedness of the stochastic Biot problem; we prove the existence of a second order solution provided for random model coefficients satisfying some weak mathematical assumptions that we detail. The existence of a second order solution appears to be a new mathematical result. The third contribution is a discretization of the stochastic Biot problem using a (non-intrusive) sparse pseudo spectral projection method. This approach relies on an ensemble of deterministic simulations carried out using a Hybrid High-Order discretization. The proposed method is further tested on a set of model problems.

Acknowledgements

This work was partially funded by the Bureau de Recherches Géologiques et Minières. The work of M. Botti was additionally partially supported by Labex NUMEV (ANR-10-LABX-20) ref. 2014-2-006.