Model order reduction for Bayesian approach to inverse problems

Abstract

This work presents an approach to solve inverse problems in the application of water quality management in reservoir systems. One such application is contaminant cleanup, which is challenging because tasks such as inferring the contaminant location and its distribution require large computational efforts and data storage requirements. In addition, real systems contain uncertain parameters such as wind velocity; these uncertainties must be accounted for in the inference problem. The approach developed here uses the combination of a reduced-order model and a Bayesian inference formulation to rapidly determine contaminant locations given sparse measurements of contaminant concentration. The system is modelled by the coupled Navier-Stokes equations and convection-diffusion transport equations. The Galerkin finite element method provides an approximate numerical solution-the ’full model’, which cannot be solved in real-time. The proper orthogonal decomposition and Galerkin projection technique are applied to obtain a reduced-order model that approximates the full model. The Bayesian formulation of the inverse problem is solved using a Markov chain Monte Carlo method for a variety of source locations in the domain. Numerical results show that applying the reduced-order model to the source inversion problem yields a speed-up in computational time by a factor of approximately 32 with acceptable accuracy in comparison with the full model. Application of the inference strategy shows the potential effectiveness of this computational modeling approach for managing water quality.

Background

Hydrodynamic processes such as contaminant transport in lakes and reservoirs have a direct impact on water quality. The contaminants will appear, spread out, and decrease in concentration, etc. because of some processes such as convection, diffusion, time rate release of contaminants, and distance of travel. To simulate such processes, a coupled system of partial differential equations (PDEs) including the Navier-Stokes equations (NSEs) and contaminant transport equations needs to be solved. A better understanding of these processes is important in managing water resources effectively.

The direct or forward problems compute the distribution of contaminant directly from given input information such as contaminant location, contaminant properties, fluid flow properties, boundary conditions, initial conditions, etc. On the contrary, the inverse problems infer the unknown physical parameters, boundary conditions, initial conditions, or geometry given a set of measured data. These known data can be obtained experimental or computational. In realistic applications, data are not perfect because of error due to sensor noise. In addition, the model may contain some uncertain parameters such as wind velocity. Thus, inverse problems are generally stochastic problems.

A Bayesian inference approach to inverse problems is one way to deal with the stochastic problem. This approach takes into account both the information of the model parameters with uncertainty and the inaccuracy of data in terms of a probability distribution [1, 2]. The Bayesian approach provides a general framework for the formulation of wide variety of problems such as climate modeling [3], contaminant transport modeling [4–6], and heat transfer [7]. Under the Bayesian framework, simulated solutions need to be evaluated repeatedly over different samples of the input parameters. There are available sampling strategies associated with Bayesian computation such as Markov chain Monte Carlo (MCMC) methods [8–11]. However, if we use traditional PDE discretization methods, such as finite element or finite volume methods, the resulting numerical models describing the system will be very large and expensive to solve.

This paper presents an efficient computational approach to solve the statistical inverse problem. The approach uses the combination of a reduced-order model and a Bayesian inference formulation. The Galerkin finite element method [12] provides an approximate numerical solution - the ‘full model’ or ‘forward model’. The MCMC method is applied to solve the Bayesian inverse problems. In particular, we are interested in inferring an arbitrary source location in a time-dependent convection-diffusion transport equation, given a velocity field and a set of measured data of the concentration field at sparse spatial and temporal locations. We obtain the velocity field by solving the Navier-Stokes equations.

Solution of large-scale inverse problems can be accelerated by applying model order reduction [13–15]. The idea is to project the large-scale governing equations onto the subspace spanned by a reduced-space basis, yielding a low-order dynamical system. Proper orthogonal decomposition (POD) is the most popular method to find the reduced basis for a given set of data. The snapshot POD method, which was proposed by Sirovich [16], provide an efficient way for determining the POD basis vectors. Snapshots are solutions of a numerical simulation at selected times or choices of the parameters. The choice of snapshots should ensure that the resulting POD basis captures the most important characteristics of the system.

Since the convection term in contaminant transport equations is a velocity-dependent term, we cannot apply the standard POD reduction framework. The evaluation of the velocity-dependent convection term in the reduced system still depends on the full finite element dimension and has the same complexity as the full-order system. Thus, a system with velocity-dependent convection term needs an additional treatment to obtain an efficient reduced-order model.

This paper is outlined as follows. In the section ‘Methods’, we introduce the mathematical model, the finite element approximate technique, and the Bayesian inference formulation and its solution using an MCMC method. In section ‘Model order reduction’, we describe the model order reduction technique and the treatment of the velocity-dependent convection term. In the section ‘Results and discussion’, we use the numerical example to demonstrate the solution of the statistical inverse problem and the reduced-order model performance. We provide some concluding remarks in the final section.

Methods

In this section, the mathematical model, finite element approximation technique, Bayesian inference formulation, and MCMC method are briefly introduced.

Transport model

Consider the fluid flow through a physical domain $\mathcal{D}\subset {\mathbb{R}}^d,d=1,2,3$ with boundary Γ. The contaminant concentration c:≡c(x,t;θ) satisfies the dimensionless parabolic PDEs, boundary conditions, and initial conditions as follows:

$$\begin{array}{lcr}\frac{\mathrm{\partial c}}{\mathrm{\partial t}}+\mathbf{u}·\nabla c-\frac{1}{\mathit{\text{Pe}}}\mathrm{\Delta c}& =& f\text{in}\mathcal{D}\times [t_0,t_f],\end{array}$$

(1)

$$\begin{array}{l}c=g\text{on}{\Gamma }_D\times [t_0,t_f],\end{array}$$

(2)

$$\begin{array}{l}\frac{\mathrm{\partial c}}{\partial \mathbf{n}}=0\text{on}{\Gamma }_N\times [t_0,t_f],\end{array}$$

(3)

$$\begin{array}{l}c(\mathbf{x},t_0)=c_0\left(\mathbf{x}\right)\text{in}\mathcal{D},\end{array}$$

(4)

where $\mathbf{x}\in \mathcal{D}$ denotes the spatial coordinates, t∈ [ t₀,t _f ] denotes time, and c₀ is the given initial condition. The inlet boundary Γ _D is subjected to a Dirichlet condition, while the remainder of the boundary Γ _N =Γ∖Γ _D satisfies Neumann conditions. The external source f(x,t;θ) in Equation 1 is used here as a input parameter, where $\theta \in \mathcal{D}$ is the source location. In this study, we suppose that the contaminant source drives the transport process, given a velocity field $\mathbf{u}\in {\mathbb{R}}^d$. Thus, the Péclet number $\text{Pe}=\frac{\parallel \mathbf{u}\parallel \mathrm{L}}{\kappa }$, where κ is the diffusivity and L is the characteristic length. The velocity field is a function of x and t, u(x,t), which can be obtained by solving the Navier-Stokes equations, as follows:

$$\begin{array}{lcr}\frac{\partial \mathbf{u}}{\mathrm{\partial t}}+\left(\right.\mathbf{u}·\nabla \left)\right.\mathbf{u}-\frac{1}{\text{Re}}{\nabla }^2\mathbf{u}+\nabla p& =& {\mathbf{f}}_c\text{in}\mathcal{D}\times [t_0,t_f],\end{array}$$

(5)

$$\begin{array}{l}\nabla ·\mathbf{u}=0\text{in}\mathcal{D}\times [t_0,t_f],\end{array}$$

(6)

$$\begin{array}{l}\mathbf{u}={\mathbf{u}}_D\text{on}\Gamma \times [t_0,t_f],\end{array}$$

(7)

$$\begin{array}{l}\mathbf{u}(\mathbf{x},t_0)={\mathbf{u}}_0\left(\mathbf{x}\right)\text{in}\mathcal{D},\end{array}$$

(8)

where Re is the Reynolds number, p is the pressure field, and f _c is the body force. So, for any given f and u, one can solve for the solution c(x,t;θ).

Finite element approximations

The finite element method [12] associated with the stabilized second-order fractional-step method is applied to solve the incompressible Navier-Stokes equations (5 to 8). The detailed discussion and derivation of this method is beyond the scope of this paper. For more details, please refer to [17, 18].

We also use the finite element method for spatial discretization of the contaminant transport equations in Equations 1 to 4. This leads to the semi-discrete equations

$$\begin{array}{l}\mathbf{M}\dot{\mathbf{c}}+\left(\right.\mathbf{C}\left(u\right)+K\mathbf{c}=\mathbf{F}(t;\theta ),\end{array}$$

(9)

$$\begin{array}{l}\mathbf{c}\left(t_0\right)={\mathbf{c}}_0,\end{array}$$

(10)

with the output of contaminant solution at sensor locations in the domain

$$\begin{array}{l}\mathbf{y}=\mathit{\text{Bc}}.\end{array}$$

(11)

Here, $\mathbf{c}\left(t\right)\in {\mathbb{R}}^N$ is the discretized approximation of c(x,t) and contains N state unknowns where N is the number of grid points in the spatial discretization. $\dot{\mathbf{c}}$ is the derivative of c with respect to time. The vector $\mathbf{y}\in {\mathbb{R}}^{N_o}$ contains the N _o outputs of the system related to the state through the matrix $\mathbf{B}\in {\mathbb{R}}^{N_o\times N},N_o\ll N$. The matrices $\mathbf{M},K,C\left(u\right)\in {ℝ}^{N\times N}$ and the vector $\mathbf{F}(t;\theta )\in {\mathbb{R}}^N$ are defined in the finite element space as follows:

$$\begin{array}{l}{\mathbf{M}}_{i,j}=\underset{\mathcal{D}}{\int }{\varphi }_i\left(\mathbf{x}\right){\varphi }_j\left(\mathbf{x}\right)d\mathbf{x},\end{array}$$

(12)

$$\begin{array}{l}{\mathbf{K}}_{i,j}=\underset{\mathcal{D}}{\int }\frac{1}{\mathit{\text{Pe}}}\nabla {\varphi }_i\left(\mathbf{x}\right)·\nabla {\varphi }_j\left(\mathbf{x}\right)d\mathbf{x},\end{array}$$

(13)

$$\begin{array}{lcr}\left(\right.\mathbf{C}\left(u\right)i,j& =& \underset{\mathcal{D}}{\int }\mathbf{u}·\nabla {\varphi }_i\left(\mathbf{x}\right){\varphi }_j\left(\mathbf{x}\right)d\mathbf{x},\end{array}$$

(14)

$$\begin{array}{lcr}\left(\right.\mathbf{F}(t;\theta )j& =& \underset{\mathcal{D}}{\int }f(t;\theta ){\varphi }_j\left(\mathbf{x}\right)d\mathbf{x}.\end{array}$$

(15)

Here ϕ _i with i=1,⋯,N is the finite element basis and $\theta \in {\mathbb{R}}^d$ is the coordinate vector of the source location.

Bayesian inference to inverse problems

The Bayesian formulation of an inference problem is derived from Bayes’ theorem [1]:

$$p\left(\theta \right|\mathbf{y})=\frac{p\left(\mathbf{y}\right|\theta \left)p\right(\theta )}{p\left(\mathbf{y}\right)},$$

(16)

where $\mathbf{y}\in {\mathbb{R}}^{N_o}$ is the observed data and $\theta \in {\mathbb{R}}^d\subseteq \mathcal{D}$ is the vector of input parameters. p(θ) denotes the prior distribution of the parameters, p(y|θ) is the likelihood function, and p(θ|y) is the posterior distribution of parameters given the observed data. p(y) is the prior predictive distribution of y and may be set to an unknown constant that is not needed in the MCMC solution method. In this case, we write Bayes’s formula as

$$p\left(\theta \right|\mathbf{y})\propto p(\mathbf{y}\left|\theta \right)p\left(\theta \right).$$

(17)

A simple model for the likelihood function can be obtained from the following relationship,

$$\mathbf{y}=G\left(\theta \right)+\eta ,$$

(18)

where η is the noise which we assume to be white Gaussian noise η∼N(0,σ²I). The input-output relative in Equations 1 to 4 is denoted as the forward model G(θ). Subsequently, the likelihood can be written as

$$p\left(\mathbf{y}\right|\theta )=\frac{1}{\sigma \sqrt{2\pi }}exp\left(-\frac{1}{2{\sigma }^2}\parallel \mathbf{y}-G\left(\theta \right){\parallel }^2\right).$$

(19)

There are many approaches to define the prior information, such as Gaussian Markov random fields or uniform distribution models, etc. In this work, we assume that our only prior information on the source location is given by the bounds on the domain. Thus, using the principle of maximum entropy [19], we take our prior to be a uniform distribution. Therefore, Equation 17 becomes

$$p\left(\theta \right|\mathbf{y})\propto \prod _{i=1}^{K}exp\left[-\frac{1}{2{\sigma }^2}{\left({\mathbf{y}}_i-G{\left(\theta \right)}_i\right)}^T\left({\mathbf{y}}_i-G{\left(\theta \right)}_i\right)\right].$$

(20)

Here, K is the number of time steps over which we collect the output data.

Markov chain Monte Carlo for posterior sampling

MCMC simulation provides a sampling strategy from a proposal distribution to the target distribution using the Markov chain mechanism. Hastings introduced the Metropolis-Hastings (MH) algorithm [10], which is a generalized version of the Metropolis algorithm [8, 9], in which we can apply both symmetric and asymmetric proposal distributions. In this work, the MH algorithm is used to solve the Bayesian inverse problem. The MH algorithm is summarized as follows:

• 1. Initialize the chain θ ⁰ and set n=0

• 2. Repeat

• n=n+1

• Generate a proposal point θ^∗∼q(θ^∗|θ)

• Generate U from a uniform U(0,1) distribution

• Update the state to θⁿ⁺¹ as

  $$\begin{array}{l}{\theta }^{n+1}=\left\{\begin{array}{cc}{\theta }^{\ast },& \text{if}\beta <\mathrm{U}\\ {\theta }^n,& \text{otherwise}\end{array}\right.\end{array}$$

• 3. Until n=N _mcmc→ stop.

Here, β is the acceptance-rejection ratio, given by

$$\beta =min\left(1,\frac{\pi \left({\theta }^{\ast }\right)q\left({\theta }^{n-1}|{\theta }^{\ast }\right)}{\pi \left({\theta }^{n-1}\right)q\left({\theta }^{\ast }\right)|{\theta }^{n-1}}\right).$$

(21)

N_mcmc is the total number of samples, π(θ) is the target distribution, q(θ^∗|θ) is a proposal distribution, and U is a random number.

The discretized model in the form of ordinary differential equations (ODEs) (Equations 9 to 11) may have very large dimensions and be expensive to solve. The MCMC method requires evaluating repeatedly the solution of the forward model (many thousands or even millions of times). Hence, these simulations can be a computationally expensive undertaking. In such situations, the reduced-order model is needed to approximate the large-scale model, which allows efficient simulations.

Model order reduction

This section presents the model order reduction framework. This includes the reduction via projection, the proper orthogonal decomposition method, the additional treatment for the velocity-dependent convection term, and the error estimation.

Reduction via projection

A reduced-order system approximating the ODEs (9 to 11) can be obtained by approximating the full state vector c as a linear combination of m basis vectors as follows:

$$\mathbf{c}\approx V{\mathbf{c}}_r,$$

(22)

where ${\mathbf{c}}_r\in {\mathbb{R}}^m$ is the reduced order state and $V=[v_1{v}_2\cdots v_m]\in {\mathbb{R}}^{N\times m}$ is an orthonormal basis, i.e., V^TV=I. The low-order model can be derived by projecting the governing system (9 to 11) onto the reduced space formed by the column span of basis V. This yields

$$\begin{array}{lcr}{\mathbf{M}}_r{\dot{\mathbf{c}}}_r+\left(\right.{\mathbf{C}}_r\left(\mathbf{u}\right)+{\mathbf{K}}_r\left)\right.{\mathbf{c}}_r& =& {\mathbf{F}}_r(t;\theta ),\end{array}$$

(23)

$$\begin{array}{l}{\mathbf{c}}_r\left(t_0\right)={\mathbf{c}}_{0r},\end{array}$$

(24)

$$\begin{array}{l}{\mathbf{y}}_r={\mathbf{B}}_r{\mathbf{c}}_r,\end{array}$$

(25)

where

$$\begin{array}{l}{\mathbf{M}}_r=V^T\mathbf{M}V,\end{array}$$

(26)

$$\begin{array}{l}{\mathbf{K}}_r=V^T\mathbf{K}V,\end{array}$$

(27)

$$\begin{array}{l}{\mathbf{C}}_r\left(\mathbf{u}\right)=V^T\mathbf{C}\left(\mathbf{u}\right)V,\end{array}$$

(28)

$$\begin{array}{lcr}{\mathbf{F}}_r(t;\theta )& =& V^T\mathbf{F}(t;\theta ),\end{array}$$

(29)

$$\begin{array}{l}{\mathbf{B}}_r=\mathbf{B}V,\end{array}$$

(30)

$$\begin{array}{l}{\mathbf{c}}_{0r}=V^T{\mathbf{c}}_0.\end{array}$$

(31)

The model reduction task is to find a suitable basis V so that m≪N and the reduced-order model yields accurate results. This study will consider POD as the method to compute the basis.

Proper orthogonal decomposition

POD provides a method to compute the reduced-order basis V and construct the low-order system. Here, we briefly describe the general POD method (more details may be found in [16]).

Supposed that we have the snapshot matrix $X\in {\mathbb{R}}^{N\times Q}$, where Q is the number of snapshot solutions, which are built from the instantaneous solution c^s(t _k ) of the forward model, corresponding to s=1,⋯,S random input source locations, picked up at k=1,⋯,T different time steps (so that Q=S T). The POD basis vectors are the m left singular vectors of X corresponding to the largest singular values (m≤Q). Let σ _i ,i=1,2,⋯,Q be the singular values of X in decreasing order. We choose the number of POD basis vectors to retain in the reduced-order model as m≤Q so that

$$\sum _{i=1}^m{\sigma }_i^2/\sum _{j=1}^Q{\sigma }_j^2\ge {\epsilon }_E,$$

(32)

where ε _E (%) is the required amount of energy, typically taken to be 99% or higher.

After obtaining the POD basis vectors for the contaminant, we have defined our reduced-order system. However, as mentioned earlier, the dimension of the system of Equations 24 to 26 is reduced only in state (concentration). The reduced convection matrix (C _r (u)) in Equation 29 still depends on the full dimension of the velocity field as in Equation 14. This means that we have to re-compute the full convection matrix (C(u)) at each time-step before projecting it onto the reduced-space basis to obtain the reduced convection matrix (C _r (u)). Hence, we need an additional treatment to avoid this computational cost.

Linear expansion techniques for velocity-dependent term

We proceed by representing the velocity field as a linear combination of POD velocity basis vectors. The velocity-dependent convection matrix can then be expanded based on this velocity expansion. In particular, let u(x,t) be a given velocity field and $\left\{\right.\mathbf{u}(\mathbf{x},t_k)k=1T$ be a set of ‘snapshots’ of velocities which are obtained from Equations 5 to 8 and taken over T time steps. The velocity field is decomposed as

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}_m\left(\mathbf{x}\right)+{\mathbf{u}}^{\prime }(\mathbf{x},t),$$

(33)

where ${\mathbf{u}}_m\left(\mathbf{x}\right)=\frac{1}{T}\sum _{k=1}^T\mathbf{u}(\mathbf{x},t^k)$ is the mean velocity field and u^′(x,t) is the fluctuating velocity field. The fluctuating velocity field is then represented by

$${\mathbf{u}}^{\prime }(\mathbf{x},t)=\sum _{k=1}^T{\alpha }_k\left(t\right){\Psi }_k\left(\mathbf{x}\right),$$

(34)

where Ψ _k (x) is the k th POD velocity basis and α _k (t) is the corresponding time dependent amplitude. We now consider the expansion of the velocity field as

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}_m\left(\mathbf{x}\right)+\sum _{k=1}^{N_u}{\alpha }_k\left(t\right){\Psi }_k\left(\mathbf{x}\right),$$

(35)

where N _u ≪T is the number of POD basis vectors use to represent the velocity.

Galerkin projection of the time-dependent velocity term in Equation 1 has the form as

$$\begin{array}{l}(\mathbf{u}·\nabla c,V_i)=\left(\left[{\mathbf{u}}_m+\sum _{k=1}^{N_u}{\alpha }_k{\Psi }_k\right]·\nabla c,V_i\right)\\ =\left[{\mathbf{u}}_m·\nabla c,V_i\right]+\left(\sum _{k=1}^{N_u}{\alpha }_k{\Psi }_k·\nabla c,V_i\right).\end{array}$$

(36)

The discrete form of Equation 37 then becomes

$$\begin{array}{l}{\mathbf{C}}_r\left(\mathbf{u}\right)=\left[V_i,{\mathbf{u}}_m·\nabla V_j\right]+\sum _{k=1}^{N_u}{\alpha }_k\left[V_i,{\Psi }_k·\nabla V_j\right]\\ ={\mathbf{C}}_{\mathit{\text{rm}}}\left(\right.{\mathbf{u}}_m\left)\right.+\sum _{k=1}^{N_u}{\alpha }_k{\mathbf{C}}_{\mathit{\text{rk}}}\left(\right.{\Psi }_k\left)\right..\end{array}$$

(37)

Here, C _{r m} (u _m (x)) and C _{r k} (Ψ _k (x)),k=1,⋯,N _u are the reduced-order forms of the full-order convection matrices C(u _m (x)) and C(Ψ _k (x)), respectively. They are computed only once in the offline step. Thus, the first case of the reduced-order model (with Equation 29) is only reduced in state (concentration), but the second case (with Equation 38) is now reduced in both state (concentration) and parameter (velocity).

Error estimation

In order to estimate the efficiency of the reduced model relative to the full model, we use the relative errors of state solutions and relative errors of outputs. These errors are defined as follows:

$$\begin{array}{lcr}{\epsilon }_F\left(t^k,\theta \right)& =& \frac{\parallel \mathbf{c}\left(t^k;\theta \right)-V{\mathbf{c}}_r\left(t^k;\theta \right){\parallel }_{L_2\left(\mathcal{D}\right)}}{\parallel \mathbf{c}\left(t^k;\theta \right){\parallel }_{L_2\left(\mathcal{D}\right)}},\end{array}$$

(38)

$$\begin{array}{lcr}{\epsilon }_s\left(t^k,\theta \right)& =& \frac{\parallel \mathbf{y}\left(t^k;\theta \right)-{\mathbf{y}}_r\left(t^k;\theta \right){\parallel }_{L_2\left(\mathcal{D}\right)}}{\parallel \mathbf{y}\left(t^k;\theta \right){\parallel }_{L_2\left(\mathcal{D}\right)}}.\end{array}$$

(39)

Here, c(t^k;θ),c _r (t^k;θ),1≤k≤T are the full and reduced solutions, while y(t^k;θ) and y _r (t^k;θ),1≤k≤T are the full and reduced outputs of interest, respectively.

Results and discussion

We consider a numerical example based on a 2D mathematical model. The velocity field in the reservoir is obtained by solving a 2D laterally averaged Navier-Stokes model. The Bayesian formulation of the inverse problem is then solved to determine an uncertain contaminant source location. We compare the effectiveness of solving the inverse problem using both the full model and reduced model techniques.

Model setup

The physical domain is illustrated in Figure 1a, which represents a simplified model of a 2D laterally averaged reservoir system. The reservoir system includes a main reservoir section and the river connections. We assume that the contaminant is present within the main reservoir section and that the contaminant transport processes are mainly affected by the inflow and wind velocity. We assume that other factors (such as heat exchange, etc.) have little influence on the processes. The spatial domain is discretized with a finite element mesh, as shown in Figure 1b. The total number of grid points is N=1,377 and the total number of elements is N _E =2,560. The simulation is run from initial time t₀=0 to final time t _f =40, with time-step size Δ t=0.08. Thus, the number of time steps is T=500.

Figure 1

The domain of interest. (a) Physical domain. (b) Computational domain.

A time-dependent velocity field is obtained from the 2D laterally averaged system, which is given in Equations 5 to 8, where the body force f _c = [f _{c x} ; f _{c z} ]^T= [ 0; g]^T with g the gravitational acceleration.

The boundary conditions are set up as follows (refer to Figure 1a for the boundary definitions):

$$\begin{array}{lcr}(u,w)& =& (1,0)\text{on}{\Gamma }_{12},\end{array}$$

(40)

$$\begin{array}{lcr}(u,w)& =& \left(-16\ast (2.0+z)\ast (1.5+z),0\right)\text{on}{\Gamma }_4,\end{array}$$

(41)

$$\begin{array}{lcr}(u,w)& =& \left(-16\ast (1.0+z)\ast (0.5+z),0\right)\text{on}{\Gamma }_8,\end{array}$$

(42)

$$\begin{array}{lcr}(u,w)& =& (0.03V_a,0)\text{on}{\Gamma }_{11},\end{array}$$

(43)

$$\begin{array}{l}p=0\text{on}{\Gamma }_{11}.\end{array}$$

(44)

The velocity on remaining boundaries are set to zero. Here, V _a is the wind speed at 10 m above the water surface. In this example, we assumed that V _a =2 m/s for the entire simulation time. The Reynolds number is set to Re=5.0e⁵, and a mixing length turbulence model is used [20]. The full model system is given in Equations 9 to 10. The Crank-Nicolson method [21] is used to discretize the system in time.

The setup for the contaminant transport part of the problem is as follows. We use a contaminant source that is defined as the superposition of Gaussian sources, each one active on the time interval t_0k∈ [ t₀,t_off] and centered at ${\theta }_k\in \mathcal{D}$, with strength h _k and width σ _{s k} :

$$f(\mathbf{x},t;\theta )=\sum _{k=1}^{n_s}\frac{h_k}{2\pi {\sigma }_{\mathit{\text{sk}}}^2}exp\left(-\frac{|{\theta }_k-\mathbf{x}{|}^2}{2\underset{\mathit{\text{sk}}}{\overset{2}{\sigma }}}\right)\delta (t-t_{0k}).$$

(45)

Here, to simplify the problem, we choose the number of sources to be n _s =1, located at θ₁=(x _c ,z _c ), with strength h₁=0.2 and width σ_{s 1}=0.05. The active time of the source is t₀₁∈ [ 0,t_off] with t_off=10. The inflow boundary and other solid boundaries satisfy a homogeneous Dirichlet condition on the contaminant concentration; the outflow boundaries and free surface boundary satisfy a homogeneous Neumann condition. The diffusivity coefficient is assumed to be constant, κ=0.005. Thus, the Péclet number $\text{Pe}=\frac{\parallel \mathbf{u}\parallel \mathrm{L}}{\kappa }=100$, where the length of the inflow section is used as the characteristic length L=0.5. The contaminant concentration is assumed to be zero at initial time t₀=0.

With this setup, we have specified our forward problem with the input parameter $\theta \in {\mathbb{R}}^2$ in the range Θ:= [ 0, 2]× [ −2, 0]. Figure 2 shows the contaminant solution c(x,t;θ) of the forward model with θ₁=(0.5,−0.5) at specific times. The contaminant field grows while the source is active. After the shutoff time of the source, the contaminant moves away, spreads out, and decreases in concentration due to convection and diffusion, until it flows out of the domain.

Figure 2

Contaminant field of forward model at specific times. Figures are taken at (a) t=4.0s, (b) t=8.0s, (c) t=16.0s and (d) t=32.0s.

The outputs of interest are the values of contaminant solution c(x,t;θ) at selected sensor locations in the computational domain. These sensors are located on a 4×4 uniform grid covering the reservoir domain as shown in Figure 1b.

Model order reduction

To perform the model order reduction, we need to obtain the POD expansion of velocity field.

Figure 3a shows all eigenvalues of the velocity field. As the number of eigenvalues increases, the significant of eigenvalues decreases, respectively. Figure 3b shows the relative errors of representing the velocity field with various numbers of POD velocity modes. The goal of this step is to find a suitable N _u that can provide a fast computation of the velocity-dependent convection term and to ensure that the relative errors are acceptable. In this example, we choose the number of terms in the velocity expansion to be N _u =12 yields the relative error around 5×10⁻³.

Figure 3

Properties of POD velocity basis vectors: (a) eigenvalues and (b) relative error analysis.

Figure 4 illustrates the first four POD basis vectors of the velocity field. Clearly, they represent the most important characteristics of the velocity field.

Figure 4

The first four POD basis vectors of the velocity field. (a) First POD velocity basis vector. (b) Second POD velocity basis vector. (c) Third POD velocity basis vector. (d) Fourth POD velocity basis vector.

A comparison between the original velocity and the approximate velocity representation at final simulating time is given in Figure 5. We observe that the approximate velocity profiles at x=0.5 (Figures 5a,b) and x=1.0 (Figures 5c,d) are able to represent well the characteristics of the original velocity profiles.

Figure 5

Velocity profile taken at t=t _f at specific locations. Figures are captured (a) U profile at x=0.5, (b) V profile at x=0.5, (c) U profile at x=1 and (d) V profile at x=1. Here U(full) and V(full) are the original velocity while Ur(mor) and Vr(mor) are the approximate velocity using a POD expansion with N _u =12 modes.

To generate the snapshots needed for the POD basis to represent the contaminant concentration, we choose S location samples in the computational domain. In this example, we generate random input locations $\theta \in \mathcal{D}\subset {\mathbb{R}}^2$ for source term f(x,t;θ).

Table 1 shows the relative errors and computational time ratio of the reduced-order models (ROM) for different sample sizes of snapshot matrices. Here, we use the notation ROM-1 to denote the reduced model mentioned in the section ‘Proper orthogonal decomposition’ and ROM-2 for the one mentioned in the section ‘Linear expansion techniques for velocity-dependent term’. The POD basis is chosen based on the energy as in Equation 33, with ε _E =99.999%. We observe that while the relative errors in the outputs are similar for both cases, the ratio of computational time differs significantly. The case ROM-2 shows speedups from 15 to 111 times, while the case ROM-1 is in fact more expensive to solve than the full forward model (t _{f u l l} ≈285 seconds)¹. This is because the first case consumes a lot of time in re-computing the velocity-dependent convection matrix and then projecting it onto the reduced space at each time step. We choose as our reduced-order model the case ROM-2 with sample size S=30.

Table 1 The properties of various reduced-order models

Figure 6a presents the eigenvalues of contaminant snapshot matrix for the sample test case S=30. We observe that POD eigenvalues of the contaminant field are decayed very fast, and they are significantly important only at the first few hundred values. The first four POD contaminant basis vectors are shown in Figure 7.

Figure 6

Properties of POD contaminant basis vectors: (a) eigenvalues and (b) relative error analysis.

Figure 7

The first 4 POD basis vectors of the contaminant field. (a) First POD contaminant basis vector. (b) Second POD contaminant basis vector. (c) Third POD contaminant basis vector. (d) Fourth POD contaminant basis vector.

Figure 6b shows relative errors of ROMs for different numbers of POD basis vectors. With the energy taken to 99.99999%, the ROM resulted in 434 POD basis vectors with relative error around 4×10⁻³. In the trade-off between accuracy results and computational efforts (please refer to Table 1), we can choose our ROM with 99.999% of energy captured which resulting in the size of m=257 POD basis vectors and an acceptably small error.

Figure 8a shows the outputs of the full model (y) and ROM-2 (y _r ) with size m=257 at selected sensor locations. The selected sensors (sensor 6, 7, 10, and 11) capture more information of the contaminant than the others for the case studied. We observe that ROM-2 is able to capture well the behavior of the full model at the sensor locations. This reduced-order model will be utilized as the forward solver in solving the inverse problem.

Figure 8

A comparison of the full model (N=1377) and ROM (m=257) output of interest. (a) Output comparison and (b) a set of generated data at sensors 6, 7, 10 and 11.

Inverse problems

MCMC is now used to solve the inverse problem for a variety of source locations using the ROM-2 solver above. By accounting for both measured noise and uncertain information in model parameters, Bayesian inference to inverse problem leads to a well-posed problem resulting in a posterior distribution of the unknown [1, 2]. Once obtaining the posterior distribution sample, all statistical questions related to the unknown could be answered with sample averages.

To represent the behavior of uncertain variables such as wind velocity, we generate noisy data at sensor locations. The ‘real’ data (y_real) is generated synthetically by adding noise to the ideal data (y_ideal) as in Equation 18. The noise is assumed to be Gaussian η∼N(0,σ²I). Figure 8b shows the noisy data at sensor locations with σ=0.3.

We use the snapshots with S=30 samples above to generate several reduced-order models with different energies. To estimate the relative error of the inverse problem solutions (ε _θ ), we choose the solution corresponding to the largest order as a ‘truth’ solution. We do the MCMC simulation with the starting point θ_ini=(0.821; −0.955)^T. The total number of MCMC samples is set to N_mcmc=5000. The initial burn-in period is set to N_burnin=500. After this stage, data is saved to compute summary statistics of source locations.

Table 2 shows the results of the inverse source problems with different numbers of POD basis. For the case with m=257, the computational time is around 24 h to obtain the solution. If we use the full forward solver, with the speedup factor as given in Table 1, the computational time is estimated around 773 h or approximately 32 days.

Table 2 Estimated inverse solutions for different numbers of POD basis

Figure 9a shows the credible regions of the source locations θ. In this figure, both a pairwise scatter plot and 1D marginal distributions are displayed. The dash line (in θ₁ axis) and dash-dot line (in θ₂ axis) show the probability density function of each parameter while the solid blue regions show the 95% and 99% credible regions in which the parameters appeared the most (the blue dots).

Figure 9

Credible regions and historical trajectory plots of the MCMC convergence. (a) A parameter plot showing estimated source location with credible regions. (b) A contour plot of the posterior pdf.

To assess the convergence of the MCMC approach, we generate four different initial points for the MH sampler. Figure 9b shows the contour plot of the posterior probability density function (the contour plot corresponds to the result of Init1) and the historical trajectories of MCMC samplers for each starting point. We observe that for all four initial points, the MCMC simulations converge to the same final solution. However, depending on the position of the initial points, the computational time varies, as given in Table 3.

Table 3 Estimated inverse solutions for different MCMC starting points

Conclusion

This study has applied successfully the combination of a model order reduction technique based on the POD and a Bayesian inference approach to solve an inverse problem that seeks to identify an uncertain contaminant location. Applying an additional POD expansion to approximate the velocity results in a reduced-order model in both state (concentration) and parameter (velocity). The resulting reduced-order model is efficient for solution of the forward problem. Solution of the Bayesian formulation of the inverse problem using the reduced-order solver is much more rapid than using the full model, yielding the probability density of the source location in reasonable computational times. The computational time of using a reduced-order model with m=257 degrees of freedom is about a factor of 32 times lower than using the full model with size N=1,377. This reduction is important in real-time water quality management applications because it reduces time cost and storage requirements.

Endnote

¹ The simulations were performed on a personal computer (PC) with processor Intel(R) Core(TM)2 Duo CPU E8200 @2.66GHz 2.66GHz, RAM 3.25GB, 32-bit Operating System.