Use of Deep-Learning-Accelerated Gradient Approximation for Reservoir Geological Parameter Estimation

Abstract

The estimation of space-varying geological parameters is often not computationally affordable for high-dimensional subsurface reservoir modeling systems. The adjoint method is generally regarded as an efficient approach for obtaining analytical gradient and, thus, proceeding with the gradient-based iteration algorithm; however, the infeasible memory requirement and computational demands strictly prohibit its generic implementation, especially for high-dimensional problems. The autoregressive neural network (aNN) model, as a nonlinear surrogate approximation, has gradually received increasing popularity due to significant reduction of computational cost, but one prominent limitation is that the generic application of aNN to large-scale reservoir models inevitably poses challenges in the training procedure, which remains unresolved. To address this issue, model-order reduction could be a promising strategy, which enables us to train the neural network in a very efficient manner. A very popular projection-based linear reduction method, i.e., propel orthogonal decomposition (POD), is adopted to achieve dimensionality reduction. This paper presents an architecture of a projection-based autoregressive neural network that efficiently derives an easy-to-use adjoint model by the use of an auto-differentiation module inside the popular deep learning frameworks. This hybrid neural network proxy, referred to as POD-aNN, is capable of speeding up derivation of reduced-order adjoint models. The performance of POD-aNN is validated through a synthetic 2D subsurface transport model. The use of POD-aNN significantly reduces the computation cost while the accuracy remains. In addition, our proposed POD-aNN can easily obtain multiple posterior realizations for uncertainty evaluation. The developed POD-aNN emulator is a data-driven approach for reduced-order modeling of nonlinear dynamic systems and, thus, should be a very efficient modeling tool to address many engineering applications related to intensive simulation-based optimization.

1. Introduction

In the field of reservoir engineering, the geological parameters, e.g., rock permeability and porosity, cannot be directly measured and, thus, are often uncertain. The accurate predictions of well production dynamics, e.g., oil production and oil distribution underground, largely depend on whether we can obtain these geological parameters. Thus, it is necessary to estimate these geological parameters to reduce the uncertainty of the oil well production dynamic predictions. Due to the necessity of many high-fidelity model simulations, solutions of inverse modeling or parameter estimation for computational fluid modeling system with many spatially variable parameters are almost intractable tasks. In general, gradient-based optimization algorithms are very common choices to address parameter estimation problem with the aid of adjoint models [1]. Although the conventional adjoint modeling methodology has been extensively demonstrated to be very computationally cheap for achieving analytical gradient, the realization of adjoint models is generally very difficult, especially when the simulation models are not code-intrusive.

In general, emulating high-fidelity models with computationally efficient proxy models are gaining popularity to ease the development of the adjoint models. Proxy models are able to deliver a rapid simulation of high-fidelity models and, therefore, are capable of speeding up adjoint-based optimization algorithms. To our current best knowledge, the popular proxy modeling approaches roughly consist of reduced-order modeling (ROM) and data-driven proxy modeling [2,3]. The main idea behind ROM is to obtain a set of basis functions/vectors that can be used to compress the high-dimensional model into a low-dimensional subspace. We often achieve the set of orthogonal basis functions/vectors by the use of projection-based proper orthogonal decomposition (POD) [4,5]. Such projection-based ROM has been successfully applied to various research areas [6,7,8], and shows very promising efficiency increments in the process of model simulations.

Data-driven proxy modeling methods use the simulated or experimental data to accurately emulate the high-fidelity models. The commonly used data-driven proxy modeling methods include Gaussian process regression [9,10], support vector regression [11], and artificial neural network [12]. These aforementioned data-driven proxy modeling methods have extensively proved their applicability; however, they are usually restricted to handle small-scale models. The deep learning (DL) method has recently increased popularity due to its successful applications in various research domains, such as computer vision [13,14] and image processing [15,16]. Recent applications of DL to emulate subsurface transport dynamic have been widely presented in the literature. For instance, Jin and Durlofsky [17] developed a data-driven DL proxy model to emulate the temporal evolution in reservoir simulation problem. Tang and Durlofsky [18] established a deep recursive neural network architecture to speed up the time-serial predictions. Overall, recent advances in DL architectures and easy-to-use deep learning packages  [19,20], have substantially prompted applications of DL proxy modeling to high-dimensional subsurface transport models [21,22,23,24,25], multiphysics conjugate problems [26], fast prediction of film cooling performance [27], and super-resolution reconstruction of turbulent flows [28].

The use of proxy models that substantially accelerate physic-model simulation to address parameter estimation problems has been widely investigated in various research fields. Cardoso MA [29] integrated POD and trajectory-piecewise-linearization (TPWL), i.e., referred to as POD-TPWL, and then applied it to oil reservoir production optimization. Subsequently, Heijn T [8] further prompted the use of POD-TPWL to geological parameter calibration problems. Vermeulen and Heemink [4], Altaf et al. [5], Kaleta et al. [30] initially utilized the POD-TPWL to derive a model-reduced adjoint-based parameter estimation in the community of reservoir engineering. Recently, Xiao et al. [31,32] further combined the POD-TPWL and domain decomposition technique to address large-scale reservoir simulation and history matching problems. The above POD-TPWL is a reduced-order linear surrogate model. To efficiently address model nonlinearity, the combination of POD and the discrete empirical interpolation method (DEIM) was initiated [33,34,35]. Ref. [36] proposed a POD-DEIM approach for rapidly solving reservoir models through projecting the large-scale partial differential equations into low-order subspace and then solved them efficiently. DL-based solutions of inversion problems or parameter estimations have also recently been explored in a wide range of application domains, e.g., hydraulic engineering [22], computer vision [37], geophysical inversion [38], and underground transport [39,40,41].

We can easily envision that the use of DL to emulate large-scale simulation models inevitably suffers from many challenges, for example, training DL for a large-scale model needs infeasible computational demands and computer memory storage. To address this issue, the combinations of DL surrogate and POD-based model-order reduction have recently gained popularity in many research fields. That is, we can project the original high-dimensional model into a low-order subspace and then train POD-based DL surrogates correspondingly [42,43,44,45]. For example, Wu et al. [43] developed data-driven reduced-order modeling by integrating POD and a temporal convolutional neural network. Mohan and Gaitonde [42] accelerated the turbulent flow simulation by combining the DL model and POD model-order reduction. The long short-term memory, as a variant of standard recurrent neural network, was adopted to emulate time-dependent dynamics of turbulence flow. Cheung et al. [44] combined data-driven deep neural network and a POD-based reduced-order modeling strategy for rapid prediction of flow dynamics.

Inspired by the successful applications of DL-based surrogate models to accelerate subsurface flow simulations, the focus of this study is on exploring the potential of combining POD-based model-order reduction and DL surrogate to emulate the original high-fidelity models and their adjoint models simultaneously. The colorredadvances of research that have arisen in the community of DL particularly facilitate the development of adjoint models. When the DL proxy is fully trained, the adjoint models of the trained DL model are easily obtained since the popular automatic-differentiation operator currently becomes accessible and is readily maintained in neural-network-related packages. The simulations of DL proxy are extremely faster than the high-fidelity models and, thus, support their superiority to derive adjoint models.

In our previous work, we have proposed many conceptual DL frameworks to derive deep-learning-based adjoint models [39,46]. Some promising results presented in our preliminary work have shown great potential in combining DL and POD to prompt the construction of adjoint models inside a reduced subspace and then efficiently support gradient-based parameter estimation [47]. In this paper, on the basis of our previous work, we systematically derive an efficient model-reduced adjoint model by the use of a DL surrogate. Following the aforementioned reduced-order modeling method, the POD-based model reduction and a nonlinear autoregressive neural network are used to compress the high-dimension model and evolve the time-series latent dynamic features, respectively. Specifically, a physically interpreted projection method, i.e., POD, is employed to obtain low-order subspace. To learn the time-series evolution of dynamic models, we use an autoregressive structure which contains a stack of residual neural network units. Once the projection-based autoregressive neural network model is constructed, the corresponding adjoint model can be easily derived by the use of automatic-differentiation mode in many popular DL packages.

In summary, we generally employed a gradient-based iterative algorithm to estimate the space-varying geological parameters. This algorithm needs us to compute the gradients. Compared to the finite-different method that needs to perturb the parameters one-by-one and, thus, requires intensive computation cost, the adjoint method can obtain the gradients analytically and, thus, is very computationally efficient. In this study, we mainly develop a deep learning surrogate model to speed up the construction of the adjoint model. The content of this study is structured as follows: Section 2 describes the mathematical formula of the general adjoint-based variational parameter estimation method. Section 3 provides the mathematical derivation of the model-reduced adjoint model. Section 4 describes an autoregressive neural network proxy model for the purpose of emulating temporal dynamics. In Section 5, many accuracy assessments of our proposed DL proxy model in terms of time-series fluid saturation prediction and parameter estimation results are presented. This section also comprehensively verifies the applicability of model-reduced adjoint-based data assimilation on a subsurface flow modeling with the aim of geological parameters’ estimation. Finally, Section 6 concludes our work and discusses possible research extensions.

2. Definition of Adjoint-Based Parameter Estimation for Reservoir Model

Reservoir geological parameter estimation is an important step in the process of reservoir development. Through historical fitting, various static parameters of the reservoir and fluid can be further clarified, thereby providing a more accurate understanding of reservoir characteristics. Reservoir history fitting is based on reservoir numerical simulation and multisource observation data fitting to modify model parameters. Without loss of generality, we simplify an explicit formula for a single simulation step of a discretized subsurface porous media flow system as follows:

$${\mathbf{x}}^n={\mathbf{f}}^n({\mathbf{x}}^{n-1},\mathbf{m}),n=1,\cdots ,N_d$$

(1)

where ${\mathbf{f}}^n$ denotes dynamic evolution model at two consecutive timesteps, n and $n-1$. ${\mathbf{x}}^n$ represents the grid-based reservoir pressure and saturation. $\mathbf{m}$ denotes a vector of spatially varying parameters, which are the spatial permeability fields in our study. We can find descriptions from reference [48].

Uncertain parameters can be estimated by minimizing an objective function based on adjoint-based optimization algorithm. In general, an adjoint-based parameter estimation constrained by the simulation model (e.g., Equation (1)) can be formulated as follows:

$$\begin{array}{cc}\hfill & \mathbf{m}=arg\underset{\mathbf{m}}{min}J({\mathbf{x}}^1,\cdots ,{\mathbf{x}}^n\cdots ,{\mathbf{x}}^{N_d},\mathbf{m})\hfill \\ \hfill & J({\mathbf{x}}^1,\cdots ,{\mathbf{x}}^n\cdots ,{\mathbf{x}}^{N_d},\mathbf{m})={\displaystyle \frac{1}{2}}\sum _{n=1}^{N_d}{\parallel {\mathbf{d}}^n-{\mathbf{x}}^n\parallel }_2^2\hfill \\ \hfill & {\mathbf{x}}^n-{\mathbf{f}}^n({\mathbf{x}}^{n-1},\mathbf{m})=\mathbf{0},n=1,\cdots ,N_d\hfill \end{array}$$

(2)

where ${\mathbf{d}}^n$ represents the vector of observations at the timestep $\mathit{n}$.

To estimate the geological parameters, the objective function $J({\mathbf{x}}^1,\cdots ,{\mathbf{x}}^{N_d},\mathbf{m})$ needs to be minimized by gradient-based local optimization or stochastic global optimization methods. Among them, gradient-based optimization supported by the adjoint-state method is one of the most generic options. In the adjoint-based optimization approach [49,50], the objective function gradient can be derived and supported by an adjoint model, which can be presented by

$${\displaystyle \frac{dJ}{d\mathbf{m}}}=-\sum _{n=1}^{N_d}{\left[{\lambda }_n\right]}^T{\displaystyle \frac{\partial {\mathbf{x}}^n}{\partial \mathbf{m}}}-\sum _{n=1}^{N_d}{\left[{\displaystyle \frac{\partial {\mathbf{x}}^n}{\partial \mathbf{m}}}\right]}^T({\mathbf{d}}^n-{\mathbf{x}}^n)$$

(3)

and the corresponding equation of the adjoint model is as follows:

$${\left[{\displaystyle \frac{\partial {\mathbf{x}}^{n+1}}{\partial {\mathbf{x}}^n}}\right]}^T{\mathit{\lambda }}^{n+1}-{\mathit{\lambda }}^n+({\mathbf{d}}^n-{\mathbf{x}}^n)=\mathbf{0},n=N_d,\cdots ,1$$

(4)

with an additional condition ${\mathit{\lambda }}_{N_d+1}$ = 0. More details related to the description of adjoint model are given in our previous work [31,51]. When the objective function gradient ${\displaystyle \frac{dJ}{d\mathbf{m}}}$ is obtained, the uncertain parameters can be progressively adjusted with the use of quasi-Newton or Gaussian–Newton optimization algorithms [52].

The derivation of athe djoint model needs two derivatives of the dynamic equation, i.e., ${\displaystyle \frac{\partial {\mathbf{x}}^n}{\partial \mathbf{m}}}$ and ${\displaystyle \frac{\partial {\mathbf{x}}^{n+1}}{\partial {\mathbf{x}}^n}}$, which are generally not easy to obtain. In addition, the adjoint codes in the commercial simulators are generally not available, which further worsens the construction of the adjoint model. In addition, the programming of adjoint model needs requires intensive computation cost and computer memory storage, especially for high-dimension models. To mitigate these issues, we propose a model-reduced neural network surrogate to derive these two derivatives for the construction of adjoint models.

3. Formula of Model-Reduced Adjoint Modeling

3.1. Model-Order Reduction Approach

The dynamic equations, e.g., ${\mathbf{f}}^n$, are generally high-dimensional and this hinders the application of the adjoint-based optimization approach to practical models. This limitation greatly motivates the popularity of the model-order reduction approach [8]. The use of model-order reduction is to represent the high-dimensional model state $\mathbf{x}$ using a set of low-dimensional latent-space features $\mathit{\psi }\in R^{N_{\psi }}$

$$\mathbf{x}={\mathbf{\Phi }}_x\mathit{\psi }$$

(5)

where ${\Phi }_x$ denotes a projection basis matrix, and $N_{\psi }$ denotes the number of low-dimensional patterns. This paper employs the POD-based model-order reduction due to the physical interpretation and high scalability.

In this study, the space-varying parameters, e.g., geological permeability, need to be estimated. In general, the spatial parameters preserve some underlying spatial correlations [53]. The geological permeability $\mathbf{m}$ is assumed to statistically obey a Gaussian distribution, and, thus, can be reduced into a low-dimensional subspace. On the basis of this assumption, the high-dimensional parameters $\mathbf{m}$ can be represented by a set of low-dimensional variables $\mathit{\xi }\in R^{N_{\xi }}$ as follows:

$$\mathbf{m}={\mathbf{\Phi }}_{\mathbf{m}}\mathit{\xi }$$

(6)

where ${\mathbf{\Phi }}_m$ denotes a basis matrix. $N_{\xi }$ represents the dimension of low-dimensional patterns. The basis matrix is generally obtained by using principle component analysis (PCA).

Through substituting Equations (5) and (6) into Equation (2), we obtain a new reduced-order model ${\mathbf{f}}_{POD}^n\in R^{N_{\psi }+N_{\xi }}\to R^{N_{\psi }}$ as follows

$${\mathit{\psi }}^n={\mathbf{f}}_{POD}^n({\mathit{\psi }}^{n-1},\mathit{\xi }),n=1,\cdots ,N_d$$

(7)

3.2. Model-Reduced Adjoint Modeling

Equations (2)–(4) represent the adjoint equations for a fluid dynamic model upon a full-order space. For any high-dimensional dynamic system, the derivation and solution of adjoint equations are very expensive. To address this issue, a reduced-order adjoint model will be introduced in detail. In general, the parameter estimation problem upon the reduced-order subspace can be reformulated as follows:

$$\begin{array}{cc}\hfill & \mathit{\xi }=arg\underset{\mathit{\xi }}{min}{J}_{rom}({\mathit{\psi }}^1,\cdots ,{\mathit{\psi }}^{N_d},\mathit{\xi })\hfill \\ \hfill & J_{rom}({\mathit{\psi }}^1,\cdots ,{\mathit{\psi }}^{N_d},\mathit{\xi })={\displaystyle \frac{1}{2}}\sum _{n=1}^{N_d}{\parallel {\mathbf{d}}^n-{\mathbf{\Phi }}_x{\mathit{\psi }}^n\parallel }_2^2\hfill \\ \hfill & {\mathit{\psi }}^n-{\mathbf{f}}_{POD}^n({\mathit{\psi }}^{n-1},\mathit{\xi })=\mathbf{0},n=1,···,N_d\hfill \end{array}$$

(8)

The Lagrange multiplier is generally used to formulate the model-reduced adjoint-state method, where the model-constrained optimization problem is replaced by an unconstrained optimization. We modify the original objective function $J_{rom}$ through adjoining the system equation as the following formula:

$$\begin{array}{c}\hfill {\widehat{J}}_{rom}=J_{rom}+\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T[{\mathit{\psi }}^n-{\mathbf{f}}_{POD}^n({\mathit{\psi }}^{n-1},\mathit{\xi })]\end{array}$$

(9)

where ${\mathit{\lambda }}_{rom}^n$ represents the model-reduced adjoint vector of variables at time instance n. Then, we can further obtain the total variation of ${\widehat{J}}_{rom}$ with respect to $d{\mathit{\psi }}^n$, for n = 1, $\cdots$, $N_d$, and $d\mathit{\xi }$ as follows:

$$\begin{array}{c}\hfill d{\widehat{J}}_{rom}=d{J}_{rom}+\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^{T}d{\mathit{\psi }}^n-\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial {\mathit{\psi }}^{n-1}}}d{\mathit{\psi }}^{n-1}-\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial \mathit{\xi }}}d\mathit{\xi }\end{array}$$

(10)

We can rearrange the above equation by changing the terms ${\mathit{\psi }}^{n-1}$

$$\begin{array}{c}\hfill d{\widehat{J}}_{rom}=d{J}_{rom}+\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^{T}d{\mathit{\psi }}^n-\sum _{n=0}^{N_d-1}{\left[{\mathit{\lambda }}_{rom}^{n+1}\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^{n+1}}{\partial {\mathit{\psi }}^n}}d{\mathit{\psi }}^n-\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial \mathit{\xi }}}d\mathit{\xi }\end{array}$$

(11)

and further reformulating the limits of the sums involving ${\mathit{\psi }}^n$ obtains

$$\begin{array}{cc}\hfill d{\widehat{J}}_{rom}& =d{J}_{rom}+\sum _{n=1}^{N_d}\{{\left[{\mathit{\lambda }}_{rom}^n\right]}^T-{\left[{\mathit{\lambda }}_{rom}^{n+1}\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^{n+1}}{\partial {\mathit{\psi }}^n}}\}d{\mathit{\psi }}^n-{\left[{\mathit{\lambda }}_{rom}^1\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^1}{\partial {\mathit{\psi }}^0}}d{\mathit{\psi }}^0+\hfill \\ & {\left[{\mathit{\lambda }}_{rom}^{N_d+1}\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^{N_d+1}}{\partial {\mathit{\psi }}^{N_d}}}d{\mathit{\psi }}^{N_d}-\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial \mathit{\xi }}}d\mathit{\xi }\hfill \end{array}$$

(12)

The sum of the variation of objective function can be represented as the total variations regarding to ${\mathit{\psi }}^n$ and $\mathit{\xi }$

$$\begin{array}{c}\hfill d{J}_{rom}=\sum _{n=1}^{N_d}{\displaystyle \frac{\partial J_{rom}}{\partial {\mathit{\psi }}^n}}d{\mathit{\psi }}^n+{\displaystyle \frac{\partial J_{rom}}{\partial \mathit{\xi }}}d\mathit{\xi }\end{array}$$

(13)

Substituting the above equation into Equation (12) yields

$$\begin{array}{c}\hfill d{\widehat{J}}_{rom}=\sum _{n=1}^{N_d}\{{\displaystyle \frac{\partial J_{rom}}{\partial {\mathit{\psi }}^n}}+{\left[{\mathit{\lambda }}_{rom}^n\right]}^T-{\left[{\mathit{\lambda }}_{rom}^{n+1}\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^{n+1}}{\partial {\mathit{\psi }}^n}}\}d{\mathit{\psi }}^n-\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial \mathit{\xi }}}d\mathit{\xi }+{\displaystyle \frac{\partial J_{rom}}{\partial \mathit{\xi }}}d\mathit{\xi }\end{array}$$

(14)

We set all coefficients of $d{\mathit{\psi }}^{N_d}$ to zero; finally, the adjoint model is given as follows:

$$\begin{array}{c}\hfill {\left[{\displaystyle \frac{\partial J_{rom}}{\partial {\mathit{\psi }}^n}}\right]}^T+{\mathit{\lambda }}_{rom}^n-{\left[{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^{n+1}}{\partial {\mathit{\psi }}^n}}\right]}^T{\mathit{\lambda }}_{rom}^{n+1}=\mathbf{0},n=1,\cdots ,N_d\end{array}$$

(15)

Given the ending condition ${\mathit{\lambda }}_{rom}^{N_d+1}=\mathbf{0}$, the reduced-order adjoint model can be easily solved through a back-propagation manner. The total variations of ${\widehat{J}}_{rom}$ regarding the parameters $\mathit{\xi }$ are given by

$$\begin{array}{c}\hfill {\displaystyle \frac{d{\widehat{J}}_{rom}}{d\mathit{\xi }}}={\displaystyle \frac{\partial J_{rom}}{\partial \mathit{\xi }}}-\sum _{n=1}^{N_d}{\left[{\mathit{\lambda }}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial \mathit{\xi }}}\end{array}$$

(16)

The model-reduced adjoint system equation (e.g., Equation (16)) needs two sensitivity terms of the reduced-order model, e.g., ${\displaystyle \frac{\partial {\mathbf{f}}_{POD}^{n+1}}{\partial {\mathit{\psi }}^n}}$ and ${\displaystyle \frac{\partial {\mathbf{f}}_{POD}^n}{\partial \mathit{\xi }}}$. In the following part, we propose a projection-based autoregressive neural network to emulate the reduced-order model ${\mathbf{f}}_{POD}^n$. This surrogate model can easily derive the sensitivity terms by the use of automatic-differentiation module built-in deep-learning frameworks.

4. Projection-Based Autogressive Neural Network (POD-aNN)

We will describe a novel methodology that derives the model-reduced adjoint model using a deep neural network in this section. In detail, it includes neural network architecture, training samples preparation, and the neural network training procedure. The diagram of our proposed neural network structure is displayed in Figure 1.

4.1. Description of Neural Network Architecture

Our proposed POD-aNN surrogate model emulates the model-reduced flow dynamic at two consecutive timesteps, e.g., Equation (7). Generally, the autoregressive neural network is capable of capturing the temporal dynamic process as follows:

$${\widehat{\mathit{\psi }}}^n={\widehat{\mathbf{f}}}_{POD}^n({\widehat{\mathit{\psi }}}^{n-1},\mathit{\xi },\mathit{\theta }),n=1,\cdots ,N_d$$

(17)

where ${\widehat{\mathit{\psi }}}^n$ and ${\widehat{\mathit{\psi }}}^{n-1}$ are the POD pattern coefficients at the n-th and $n-1$-th timestep, and $\mathit{\xi }$ is the PCA coefficients. $\mathit{\theta }$ represents the trainable neural network parameters. The structure of the autoregressive neural network reveals that the current POD patterns ${\widehat{\mathit{\psi }}}^n$ rely only on the previous POD patterns ${\widehat{\mathit{\psi }}}^{n-1}$ and the time-independent PCA patterns $\mathit{\xi }$, indicating that the autoregressive neural network is capable of emulating the temporal dynamic for subsurface flow models [23].

The purpose of POD-aNN is to perform a vector-to-vector regression task. The selection of a DL structure cannot be generally determined beforehand since the design of the DL structure adapted by specific tasks cannot follow strict criteria. In this study, the configuration of this POD-aNN structure is guided by our previous work [47]. A diagram of the POD-aNN architecture is depicted in Figure 1. Figure 1 illustrates the main architecture of the proposed POD-aNN neural network model. The basic unit of this neural network model is composed of many convolutional blocks. Each convolutional block includes three continuous operators, including a one-dimensional convolutional layer ($Conv1$), a batch normalization layer ($BN$), and a rectified linear activation unit ($ReLU$). To further improve the prediction accuracy of the POD-aNN surrogate model, we introduce two different kinds of residual blocks [54]. The details of the residual blocks are depicted in Figure 1a,b. These two kinds of residual blocks, i.e., $ResBlock1$ and $ResBlock2$, are connected and followed by a full-connected ($FC$) neural network. The deep insights/motivations about the POD-aNN architecture can be further explained here. On the one hand, the adjoint model requires us to obtain dynamic evolution at two consecutive timesteps, and the autoregressive structure, indeed, satisfies this purpose compared to the conventional recurrent neural network. As is known, the NN with deeper neural network layers can represent the complex relationship. To ease the training process or mitigate the gradient vanish/explosion issue, we often use residual block or dense blocks. The use of dense blocks might be an alternative. The choice of a neural network architecture is flexible, but it is still relatively subjective.

In this study, we propose a model-reduced nonlinear surrogate model through combining the POD and autoregressive neural network. The aNN neural network is mainly used to represent the nonlinear evolution of the model dynamics. As noted earlier, to efficiently address model nonlinearity, the combination of POD and discrete empirical interpolation method (DEIM) has been proposed in the literature [33,34,35]. Both the DEIM and aNN can be used to address model nonlinearity; however, their implementation theory is significantly different. In terms of POD-DEIM, the reduced nonlinear terms are obtained by applying the DEIM to a POD basis of the nonlinearities. The use of the DEIM technique can project the original full-order model to the low order and, thus, reduce complexity in solving equations by using Newton-like iteration methods. Overall, the essence of POD-DEIM is to project the high-dimensional governing equations into low-order subspace and prompt the computation efficiency of numerical solvers; by contrast, our proposed POD-aNN constructs the surrogate model.

4.2. Samples Preparation and Training Procedure

A set of samples needs to be generated for training the POD-aNN model. In subsurface flow models, the training samples constitute many pairs of geological permeability parameters to time-series fluid saturation/pressure realizations. After that, the full-order states $\mathbf{x}$ and parameters $\mathbf{m}$ are projected into reduced state variables $\mathit{\psi }$ and $\mathit{\xi }$ using POD and PCA, respectively. To emulate the temporal relationship for two consecutive timesteps, the time-series saturation realizations, i.e., (${\mathit{\psi }}_i^1$,$\cdots$,${\mathit{\psi }}_i^{N_d}$), for the i-th high-fidelity model ${\mathit{\xi }}_i$ (i = 1,$\cdots$,$N_s$), can be reorganized as follows:

$${\left\{({\mathit{\xi }}_i;{\mathit{\psi }}_i^1,\cdots ,{\mathit{\psi }}_i^{N_d})\right\}}_{i=1}^{N_s}\Rightarrow {\left\{({\mathit{\xi }}_i,{\mathit{\psi }}_i^{n-1};{\mathit{\psi }}_i^n)\right\}}_{i=1,n=1}^{N_s,N_d}$$

(18)

in which $N_s$ denotes the training sample size. After reorganizing the dataset, the temporal relationship at two consecutive timesteps can be represented by the autoregressive neural network structure.

To improve the performance of the trained POD-aNN surrogate model, we introduce a regularized loss function. The first several POD patterns mainly dominate the dynamic features, and it is advisable to improve the fitting accuracy of these POD-based variables corresponding to large eigenvalues. We define the loss function by projecting the reduced states $\mathit{\psi }$ back to full-order states $\mathbf{x}$. This loss function can be considered to be an extension of the above one weighted by basis matrix ${\mathbf{\Phi }}_x$

$$\begin{array}{c}\hfill L\left(\mathit{\theta }\right)={\displaystyle \frac{1}{N_s{N}_t}}\sum _{i=1}^{N_s}\sum _{n=1}^{N_t}{\parallel {\mathbf{\Phi }}_x{\widehat{\mathit{\psi }}}_i^n-{\mathbf{\Phi }}_x{\mathit{\psi }}_i^n\parallel }_2^2\end{array}$$

(19)

We employ the stochastic-gradient optimizer, i.e., $\mathit{Adam}$, to optimize the neural network model. Once the POD-aNN is trained, the predictions for new test samples are easily realized. Given the initial model-reduced state variables ${\mathit{\psi }}^0$ and the reduced parameter $\mathit{\xi }$, Equation (17) can be repeatedly performed to obtain time-series reduced state variables ${\widehat{\mathit{\psi }}}^n$ at all $N_d$ time instances. Specifically speaking, each output vector ${\widehat{\mathit{\psi }}}^n$ is consecutively predicted through regarding the reduced state ${\widehat{\mathit{\psi }}}^{n-1}$ and reduced parameter $\mathit{\xi }$ as inputs. This procedure is explicitly displayed in Figure 2. After that, a set of predicted variables [${\widehat{\mathit{\psi }}}^1$, ${\widehat{\mathit{\psi }}}^2$, …,${\widehat{\mathit{\psi }}}^{N_d}$] can be transformed to their corresponding full-order states [${\mathbf{x}}^1$, ${\mathbf{x}}^2$, …,${\mathbf{x}}^{N_d}$] using Equation (5).

4.3. Model-Reduced Adjoint-Based Data Assimilation Using POD-aNN

The objective function $J_{POD-aNN}$ using our trained POD-aNN surrogate model can be reformulated as follows:

$$\begin{array}{cc}\hfill & \mathit{\xi }=arg\underset{\mathit{\xi }}{min}{J}_{POD-aNN}({\widehat{\mathit{\psi }}}^1,\cdots ,{\widehat{\mathit{\psi }}}^{N_d},\mathit{\xi })\hfill \\ \hfill & J_{POD-aNN}({\mathit{\psi }}^1,\cdots ,{\mathit{\psi }}^{N_d},\mathit{\xi })={\displaystyle \frac{1}{2}}\sum _{n=1}^{N_d}{\parallel {\mathbf{d}}^n-{\mathbf{\Phi }}_x{\widehat{\mathit{\psi }}}^n\parallel }_2^2\hfill \\ \hfill & {\widehat{\mathit{\psi }}}^n-{\widehat{\mathbf{f}}}_{POD}^n({\widehat{\mathit{\psi }}}^{n-1},\mathit{\xi },\mathit{\theta })=\mathbf{0},n=1,\cdots ,N_d\hfill \end{array}$$

(20)

and the gradient of the approximated objective function $J_{POD-aNN}$

$$\begin{array}{c}\hfill {\displaystyle \frac{d{J}_{POD-aNN}}{d\mathit{\xi }}}={\displaystyle \frac{\partial J_{POD-aNN}}{\partial \mathit{\xi }}}-\sum _{n=1}^{N_t}{\left[{\widehat{\mathit{\lambda }}}_{rom}^n\right]}^T{\displaystyle \frac{\partial {\widehat{\mathbf{f}}}_{POD}^n}{\partial \mathit{\xi }}}\end{array}$$

(21)

and the model-reduced adjoint equation with the model-reduced Lagrange multipliers ${\widehat{\mathit{\lambda }}}_{rom}^n$

$$\begin{array}{c}\hfill {\left[{\displaystyle \frac{\partial J_{POD-aNN}}{\partial {\widehat{\mathit{\psi }}}^n}}\right]}^T+{\widehat{\mathit{\lambda }}}_{rom}^n-{\left[{\displaystyle \frac{\partial {\widehat{\mathbf{f}}}_{POD}^{n+1}}{\partial {\widehat{\mathit{\psi }}}^n}}\right]}^T{\widehat{\mathit{\lambda }}}_{rom}^{n+1}=\mathbf{0}\end{array}$$

(22)

Since the training of the POD-aNN surrogate model and adjoint-based parameter estimation are implemented in the PyTorch framework, the auto-differentiation module can be employed to derive a set of derivative terms, e.g., ${\displaystyle \frac{\partial {\widehat{\mathbf{f}}}_{POD}^n}{\partial \mathit{\xi }}}$ and ${\displaystyle \frac{\partial {\widehat{\mathbf{f}}}_{POD}^{n+1}}{\partial {\widehat{\mathit{\psi }}}^n}}$ at $n=1,\cdots ,N_d$, which are required by the reduced-order adjoint model. In addition, another two derivatives, i.e., ${\displaystyle \frac{\partial J_{POD-aNN}}{\partial \mathit{\xi }}}$ and ${\displaystyle \frac{\partial J_{POD-aNN}}{\partial {\widehat{\mathit{\psi }}}^n}}$, of the objective function also can be computed in an efficient manner. We can efficiently compute the least-square function, e.g., Equation (20), and the derivatives using the auto-differentiation module in $PyTorch$.

When we obtain the gradient ${\left[{\displaystyle \frac{d{J}_{POD-aNN}}{d\mathit{\xi }}}\right]}^k$ at the k-th iteration, the optimized parameters at this iteration can be obtained as follows:

$$\begin{array}{c}\hfill {\mathit{\xi }}^{k+1}={\mathit{\xi }}^k-{\epsilon }^k{\left[{\displaystyle \frac{d{J}_{POD-aNN}}{d\mathit{\xi }}}\right]}^k\end{array}$$

(23)

where ${\epsilon }^k$ denotes a step-length at the k-th iteration. The gradient descent algorithm is considered to have sufficiently converged when the predefined metrics satisfy the following:

• The changes of objective functions at two consecutive iterations are very small, i.e.,

  $${\displaystyle \frac{|J_{POD-aNN}\left({\mathit{\xi }}^{k+1}\right)-J_{POD-aNN}\left({\mathit{\xi }}^k\right)|}{max\left\{\right|J_{POD-aNN}\left({\mathit{\xi }}^{k+1}\right)\left)\right|,1\}}}<{\eta }_{J_{POD-aNN}}.$$

  (24)
• The changes of updated parameters at two consecutive iterations are very small, i.e.,

  $${\displaystyle \frac{|{\mathit{\xi }}^{k+1}-{\mathit{\xi }}^k|}{max\left\{\right|\left({\mathit{\xi }}^{k+1}\right|,1\}}}<{\eta }_{\mathit{\xi }}$$

  (25)

where ${\eta }_{J_{POD-aNN}}$ and ${\eta }_{\mathit{\xi }}$ denote the predefined convergent metrics. In our experiments, ${\eta }_{J_{POD-aNN}}={10}^{-4}$ and ${\eta }_{\mathit{\xi }}={10}^{-3}$. Algorithm 1 describes the detailed implementation procedure of the entire workflow.

Algorithm 1: Model-reduced adjoint-based parameter estimation using POD-aNN surrogate model.

5. Illustrative Example and Results

This section analyzes some numerical results based on a synthetic two-dimensional subsurface transport model with a dimension size of 20 × 20. The 2D subsurface transport model contains one pumping well and one injection well located at the two corners, respectively. In this numerical experiment, we use an open-source simulation tool, i.e., Matlab Reservoir Simulation Toolbox (MRST) [55], to run the high-fidelity model (HFM) simulations.

5.1. Description of Subsurface Transport Model

Our research mainly focuses on estimating geological parameters of a subsurface transport model with the oil/water two-phase reservoir model. In the 2D transport model, the rock permeability field is assumed to be spatially heterogeneous by following a log-Gaussian random distribution. We generate Gaussian-distributed realizations of logarithmic permeability using the Stanford Geostatistical Modeling Toolbox (SGeMS) [56]. As we mentioned previously, the log-permeability field parameters are reparameterized by PCA and, therefore, can be represented by low-dimensional PCA patterns. Specifically, we preserve the $N_{\mathit{\xi }}$ = 72 PCA pattern coefficients to accurately characterize the spatial variability and correlation of permeability fields. As illustrated in Figure 3, one realization is regarded as the true log-permeability field. The triangle and circle denote the injector and producer, respectively. The 25 white diamonds represent the data measured locations.

Table 1. Experimental model properties, including the model size, fluid properties, underground formation pressure, initial fluid saturation and well pumping schedules, and hyperparameters for training POD-aNN.

Description of Parameters	Values
Grid size, $N_x$ × $N_y$	20 × 20
Fluid density	1014 kg/m3, 859 kg/m3
Fluid viscosity	0.4 mP·s, 2 mP·s
Initial oil, water saturation	${\mathrm{S}}_o$ = 0.80, ${\mathrm{S}}_w$ = 0.20
Initial formation pressure	30 MPa
Bottom-hole pressure for pumping wells	25 MPa
Pumping rate for injection well	150 m3/day
Pumping time	1800 days
Hyperparameters
Training samples, $N_s$	400, 600, 800, 1000
Testing samples, $N_{test}$	200
Learning rate	0.0003
Optimizer	$\mathit{Adam}$
Batch size	200
Epochs number	300

lists some other model settings, including the model size along two dimensions, fluid properties, underground formation pressure, initial fluid saturation, and well pumping schedules.

5.2. Construction of Projection-Based Autogressive Neural Network

Implementation of Model-Order Reduction

As noted earlier, $N_{\mathit{\xi }}$ = 72 PCA pattern coefficients are retained to accurately characterize the geological permeability fields. A set of PCA pattern coefficients are randomly generated to reconstruct the original permeability fields. Through running the MRST simulator, these permeability fields are regarded as the inputs of HFM simulations, and the simulated time-series fluid saturation snapshots are preserved as the outputs. These pairs of original permeability fields and fluid saturation constitute the training samples. In this experiment, many high-fidelity model simulations are used to collect fluid saturation snapshots. In this case study, 50 simulation models are run, which results in 1000 saturation snapshots. The basis matrix of POD patterns based on these 1000 snapshots is obtained. As illustrated in Figure 4, increasing the number of POD patterns will preserve the dynamic features with increasingly high accuracy. For example, $N_{\psi }$ = 10, 50, and 100, POD patterns retain 87.3%, 97.2%, and 98.8% accuracy, respectively. In this case study, 50 POD modes are almost sufficient to represent the original full-order flow dynamic.

The details about the POD-aNN neural network model corresponding to this 2D synthetic model are shown in

Table 2. Illustration of POD-aNN architecture in this case study.

Layer	Input Size	Output Size
Input	($N_{\xi }$ + $N_{\psi }$, 1)	(140, 2)
$ResBlock1$	(140, 2)	(140, 2)
$ResBlock2$	(140, 2)	(130, 2)
$ResBlock1$	(130, 2)	(130, 2)
$ResBlock2$	(130, 2)	(120, 2)
$ResBlock1$	(120, 2)	(120, 2)
$ResBlock2$	(120, 2)	(240, 1)
$FC$	(240, 1)	($N_{\psi }$, 1)

. As illustrated in Table 1, the dimension size of this 2D model are $N_x$ = 20 and $N_y$ = 20. In this case study, we set the numbers of channels as follows: ${\mathrm{C}}_1$ = $N_{\psi }$ + $N_{\xi }$, ${\mathrm{C}}_2$ = 140, ${\mathrm{C}}_3$ = 130, ${\mathrm{C}}_4$ = 120, and ${\mathrm{C}}_5$ = $N_{\psi }$, respectively.

5.3. Accuracy Evaluation of POD-aNN Neural Network Model

In terms of the POD-aNN neural network surrogate model, there are two main sources of errors (SOEs), including (SOE1) approximation errors of the neural network surrogate, and (SOE2) the loss of hidden physics dynamics due to an insufficient number of preserved POD patterns, e.g., $N_{\psi }$ = 10. To evaluate the quality of the POD-aNN neural network model with respect to the number of training samples $N_s$, $N_{test}$ high-fidelity models are run to generate the test samples. The mean relative error over $N_d$ timesteps among $N_{test}$ test samples, denoted as ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$, respectively, are given by

$${\gamma }_{SOE1}={\displaystyle \frac{1}{N_{test}{N}_d}}\sum _{i=1}^{N_{test}}{\displaystyle \frac{\parallel {\widehat{\mathit{\psi }}}_i^n-{\mathit{\psi }}_i^n\parallel }{{\mathit{\psi }}_i^n}},,n=1,\cdots ,N_d$$

(26)

and

$${\gamma }_{SOE2}={\displaystyle \frac{1}{N_{test}{N}_d}}\sum _{i=1}^{N_{test}}{\displaystyle \frac{\parallel {\mathbf{\Phi }}_x{\widehat{\mathit{\psi }}}_i^n-{\mathbf{x}}_i^n\parallel }{{\mathbf{x}}_i^n}},,n=1,\cdots ,N_d$$

(27)

where ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ correspond to reduced state variables and full-order variables, respectively. The ${\gamma }_{SOE1}$ only reveals the errors from the aNN neural network, while ${\gamma }_{SOE2}$ includes the total errors from both model-order reduction and the aNN neural network. The performance of surrogate model depends on the selection of hyperparameters, e.g., learning rate (lr) and batch size (bt). The evolutions of loss function corresponding to different learning rate and batch size are shown in Figure 5. The optimal hyperparameters involving the training of the POD-aNN neural network model can be found in Table 1.

We progressively increase the number of training samples, namely, $N_s$ = 400, 600, 800, and 1000, to assess the evolution of neural network accuracy. The simulation of this synthetic subsurface flow model is at 20 timesteps. In this study, we investigate three different data collection strategies through measuring the training samples at (1) $\mathbf{S}\mathbf{1}$: n = [4, 8, 12, 16, 20]-th timestep with $N_t$ = 5; (2) $\mathbf{S}\mathbf{2}$: n = [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]-th timestep with $N_t$ = 10; (3) $\mathbf{S}\mathbf{3}$: n = [1, 2, 3, 4, …, 18, 19, 20]-th timestep with $N_t$ = 20 timsesteps. For instance, the data collection strategy $\mathbf{S}\mathbf{1}$ yields 2000, 3000, 4000, and 5000 training samples, correspondingly.

Figure 6 plots the evolution of field-average relative error ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ scores over all 200 test models with respect to the time instances. Due to the negative error accumulation effects induced by the autoregessive network structure, both ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ scores progressively increase as the time goes on. This phenomenon reveals that the prediction accuracy of the POD-aNN network model decreases over time. Although the negative error accumulation problem occurs, the POD-aNN neural network model trained by the $N_s$ = 1000 dataset is still capable of achieving relatively small ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ values. In general, increasing the number of training samples has a high possibility to mitigate the negative error-accumulation effects to a certain degree. In our methodology, this can be realized by using more model simulations (i.e., 1000 high-fidelity model runs) or collecting data more frequently (i.e., $N_t$ = 20 at $\mathbf{S}\mathbf{3}$ scheme in this case study). For example, using data measured at all 20 timesteps allows us to obtain a POD-aNN network model trained by sufficient samples, which helps to mitigate the overfitting problem and, hence, generalizes the trained POD-aNN in the test samples with low relative errors.

Figure 7 displays the evolution of ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ values corresponding to $N_{test}$ = 200 test samples. In this numerical experiment, the results for $N_{\psi }$ = 100 and 50 POD modes are analyzed. It clearly reveals that the POD-aNN neural network model could produce relatively small errors, e.g., 0.02–0.10, corresponding to solely 400 training samples. When we increase the number of training samples $N_s$ to 1000, both the ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ scores can reduce to 0.06. In this POD-aNN neural network model, the current output relies on state information at a previous timestep. This dependency might cause error accumulation, that is, the approximation error gradually increases as the time evolves. By retaining more POD patterns, the difference between ${\gamma }_{SOE1}$ and ${\gamma }_{SOE2}$ becomes smaller. The main source of approximation error stems from the POD-aNN surrogate model.

Figure 8 shows the time-series saturation distributions predicted from POD-aNN and HFM at five, i.e., n = [4, 8,12, 16, 20]-th, timesteps. We should emphasize that the continuous three rows represent the saturation predictions from the high-fidelity model (x), POD-aNN ($\widehat{\mathbf{x}}$), and their differences (x-$\widehat{\mathbf{x}}$). In this numerical experiment, we generate the POD-aNN models by $N_s$ = 1000 training samples. The POD-aNN neural network model is able to accurately predict the time-series saturation indicated by small difference values (x-$\widehat{\mathbf{x}}$). For example, the mean relative errors ${\gamma }_{SOE2}$ are 6.22% and 5.34%, corresponding to $N_{\psi }$ = 50 and 100 POD modes, respectively.

Figure 9 shows the time-series prediction of fluid saturation at three spatial points. As shown in Figure 3, these three positions are the injector $I1$, producer $P1$, and a random point D. We can clearly observe that the predictions of fluid saturation from the POD-aNN are very close to the predictions from the high fidelity model. Increasing the number of POD modes generally results in better results.

5.4. Results of Geological Parameter Estimation

The accuracy of the POD-aNN surrogate model was demonstrated in the previous section. The POD-aNN is further applied in this section to estimate geological parameter estimation for the 2D synthetic subsurface transport model. Two scenarios are designed to demonstrate the performance of the POD-aNN surrogate model. Scenario 1 is mainly used to verify the proposed methodology with spatially sparse measurements, e.g., the 25 white diamonds shown in Figure 4. The purpose of Scenario 2 is to demonstrate the effectiveness and robustness of the POD-aNN surrogate model in high-dimensional inverse modeling by assimilating spatially grid-based measurements.

The root-mean-square error of the estimated geological permeability field is used to assess the results of parameter estimation, which can be generally defined by

$$\begin{array}{c}\hfill e_m=\sqrt{{\displaystyle \frac{1}{N_{\mathbf{m}}}}{\sum }_{i=1}^{N_{\mathbf{m}}}{({\mathbf{m}}_i-{\mathbf{m}}_{true,i})}^2}\end{array}$$

(28)

where $N_{\mathbf{m}}$ denotes the number of geological permeability parameters, ${\mathbf{m}}_i$ is the feature at i-th gridblock, and ${\mathbf{m}}_{true}$ is the true geological permeability field.

5.4.1. Scenario 1—Sparse Measurements

In Scenario 1, we generate the synthetic observations by adding some noise to the simulation results of the high-fidelity model parameterized by the true geological permeability field. Both the true permeability field and the measured locations can be found in Figure 4. Based on the simulated results at the n = [4, 8, 12, 16, 20]-th timestep, we collect the fluid saturation at the measured locations, which results in 125 measurements.

The sensitivity analysis of the POD-aNN surrogate to the number of training samples is performed first. The comparison between the finite-difference method (FD) and the POD-aNN surrogate model is also conducted. The dynamic evolution of the logarithmic objective function values with respect to the number of training samples is shown in Figure 10. The objective function value generated by the adjoint-based optimization method based on POD-aNN is equivalent to or even smaller than that of the FD method. Figure 11 displays the time-series fluid saturation field predicted by the POD-aNN surrogate model. Although the estimated geological permeability fields are significantly different, the prediction results of the saturation distribution are very identical to the saturation predictions from the true model. This result further reveals that multiple local optimal solutions have been obtained. In terms of computational cost, the FD method requires more than 3000 forward model simulations, while our proposed POD-aNN assisted framework can produce comparable parameter estimation results with much smaller amount of computation cost.

The sensitivity analysis of POD-aNN surrogate to the number of POD patterns is further investigated. The number of $N_{\psi }$ = 10, 50, and 100 POD modes separately corresponding to 87.3%, 97.2%, and 98.8% model-order reduction accuracy was designed to analyze the performance of parameter estimations. Figure 12a displays the dynamic changes of logarithmic objective-function values with respect to the number of POD modes. Although the use of fewer POD modes yields a faster convergence performance, we can obtain smaller objective function values as the number of POD modes increases. The main reason behind these results can be explained by the fact that the parameter estimation results from our proposed POD-aNN approach are substantially influenced by the selection and representation of the POD patterns. In this case study, the POD-aNN surrogate model with $N_{\psi }$ = 50 POD modes obtained almost comparable objective function values with the FD method. Figure 12b–d show the updated geological permeability fields and the parameter errors $e_m$ from our proposed POD-aNN surrogate model. The main parameter features, e.g., high permeable zones, were almost reconstructed. The use of a small number of POD modes also has the potential of capturing the geological features of the parameter field. Additional features of the geological permeability field would be gradually achieved with an increasing number of POD modes. In addition, increasing the number of POD modes also produces small parameter errors.

5.4.2. Scenario 2—Assimilation of Spatially Dense Measurements

In this scenario, the synthetic observations are generated by measuring the fluid saturation at all gridblocks. This type of data can nearly mimic a complex situation, e.g., spatially dense observations interpreted from a seismic survey. We assume that the spatially fluid saturation at all gridblocks are measured at the n = [4, 8, 12, 16, 20]-th timestep, which results in 2000 measurements.

Both the logarithmic objective functions and updated log-permeability fields are shown in Figure 13. As shown in Figure 13a, increasing the number of training samples gradually results in small objective functions. In this scenario, although the FD method achieves relatively smaller objective functions than that of our proposed POD-aNN surrogate model, the POD-aNN surrogate model needs much fewer high-fidelity model simulations than that of the FD method. In addition, compared with Scenario 1, the computational cost almost does not increase with the number of measurements. This result further demonstrates the scalability of the POD-aNN approach. Figure 13b–f separately show the true, prior, and posterior log-permeability fields. Comparing to FD approach, our proposed POD-aNN obtains satisfactory posterior parameter fields, which are very identical to the true model. It is also, surprisingly, found that the conventional FD method fails to accurately estimate the true model in this case study. This result further demonstrates that our proposed method has the ability to find the global optimal solution at some degree. Figure 14 displays the temporal–spatial distribution of fluid saturation predicted from POD-aNN, HFM, and their differences. It clearly can be observed that spatially dense measurements produce highly accurate predictions of the saturation.

Figure 15a displays the dynamic changes of the logarithmic objective-function with respect to optimization steps corresponding to $N_{\psi }$ = 100, 50, and 10 POD modes, respectively. Preserving more POD modes reasonably generates smaller objective function values. In this case study, preserving $N_{\psi }$ = 100 and 50 POD modes are able to obtain comparable results. For example, retaining 50 POD modes is almost able to recover the true logarithmic permeability field, as indicated in Figure 15b–d. That is to say, 50 POD modes can contain enough information extracted from the 2000 measurements. In terms of parameter estimation in this case study, ignoring some unimportant POD patterns corresponding to small eigenvalues might introduce some positive effects, i.e., efficiently searching the optimal solutions in a reduced subspace.

Figure 16 displays the time-dependent fluid saturation predicted from POD-aNN, HFM, and their differences. The differences between HFM and POD-aNN are relatively large when we only retain 10 POD modes. However, the flow patterns, e.g., the position of fluid fronts, can be approximately captured. In the inverse modeling, these flow patterns are able to provide sufficient information for updating the parameters. These results partially explain that retaining only 10 POD modes still generates satisfactory parameter estimation.

To further prove the effectiveness of our proposed POD-aNN proxy model, by using 200 different prior models, 200 posterior realizations are obtained, correspondingly. Figure 17a,b show the box plots of the parameter errors $e_m$ in these two cases, respectively. Compared with the prior models, the parameter errors $e_m$ are greatly reduced. Since a large number of measurements is assimilated in Scenario 2, the parameter error $e_m$ reduction in Scenario 2 is greater than that in Scenario 1. The average and standard deviation (STD) of ensemble log-permeability fields are shown in Figure 17c,d. The comparison of the ensemble average parameter field clearly reveals that the POD-aNN surrogate model is able to effectively identify the geological properties of the true model. In comparison to the prior models, the posterior models are very close to the true model. In addition, the geological permeability uncertainty is significantly decreased, which can be indicated by small posterior model STD. The gradient-based parameter estimation algorithm inevitably falls into local minimums. Our proposed surrogate-based parameter estimation framework does not require many other HFM simulations. We can easily obtain multiple posterior realizations. All these results confirm that the POD-aNN method can generate many posterior parameters for the uncertainty evaluation while maintaining high scalability and efficiency.

6. Discussion

6.1. Comments on Computational Cost

In terms of computational effort, the proposed POD-aNN surrogate model needs about 0.01 s running on an NVIDIA GeForce GTX 745 GPU, which is nearly 600 times faster than that of the original high-fidelity model, e.g., about 6.0 s for this 2D synthetic case study. However, the offline training phase of the neural network model has to consider additional computational costs involving training dataset generation and network parameters optimization. It is meaningless to train the POD-aNN surrogate model unless it is definitely employed in the specific applications where many simulations are repeatedly performed. Because numerous high-fidelity model simulations are usually needed in the application of simulation-based parameter estimation, POD-aNN is undoubtedly beneficial in this situation.

In this study, since the running time of the POD-aNN surrogate model can be almost ignorable, the number of HFM simulations should be a suitable metric to quantify the computational expense. For example, training the POD-aNN surrogate model by 1000 samples requires about 18 min, which equals 200 HFM simulations. In summary, our proposed model-reduced adjoint-based inverse modeling requires about 1200 model runs; the conventional FD method, by contrast, takes 13,250 HFM simulations. It can be concluded that the use of POD-aNN significantly reduces the running times in the process of solving parameter estimation problem in this case study. Both effectiveness and efficiency of our proposed method were successfully validated through a synthetic 2D subsurface transport model. We can envision that the computational efficiency of our proposed POD-aNN becomes very substantial in the large-scale simulation-based engineering applications.

6.2. Comments on the Scalability of POD-aNN

In this study, a conceptual 2D subsurface transport model is used to demonstrate the performance, e.g., accuracy and efficiency, of our proposed POD-aNN surrogate model. It is worthy of investigating whether the proposed method can scale with higher grid resolutions and larger time horizons. Both the retained POD pattern coefficients and PCA pattern coefficients are the input of the POD-aNN surrogate model. The complexity of the neural network model, e.g., the number of trainable parameters, depends on the size of input variables. We generally should provide more training samples to train complex neural network models with satisfactory accuracy. The application models with higher grid resolutions and larger time horizons usually need to preserve more POD pattern coefficients and PCA pattern coefficients, which increases the complexity of neural network models. This limitation seemingly hinders further applications to large-scale models.

Here, we provide a possible option to scale our proposed method to practical models with higher grid resolutions than those in the provided experiments. To train the deep neural network, numerous training samples obtained from high-fidelity model simulations are needed, especially for the large-scale models. To reduce the computation burden, we can further reduce the requirement of high-fidelity model simulations by training subdomain-surrogates instead of one surrogate for the global domain. As conducted in our previous work [39,40,41], we propose a subdomain-based reduced-order modeling algorithm for the estimation of space-varying parameter patterns in numerical reservoir models. Domain decomposition is applied to construct separate approximations to the numerical model in every subdomain. Motivated by this idea, we also can combine the subdomain-based model order reduction with a deep neural network surrogate and its parallel implementation. Such an approach will have the advantage that the subdomain surrogate models can be naturally trained in parallel. For the global domain deep learning surrogate, a single GPU is used to train the deep learning surrogate model. The combination of domain decomposition and deep learning enables us to train multiple subdomain deep learning surrogate models using a large number of CPUs and multiple GPUs in parallel. Motivated by this idea, we can train the POD-aNN surrogate model at each subdomain individually. Increasing the number of subdomains enables us to train subdomain deep learning surrogate models with low complexity.

7. Conclusions

To address complex parameter estimation problems with a short timeframe, we combine reduced-order modeling and deep-learning surrogate to emulate model-reduced adjoint effectively and efficiently. In detail, a structure of projection-based autoregressive neural network, referred to as POD-aNN, is proposed. The main innovation around this type of POD-aNN is to approximate objective function gradients supported by our proposed model-reduced adjoint approach. This new approach does not require us to code an adjoint model and, thus, results in very efficient implementations of gradient-based inverse modeling compared to the conventional finite-difference method. In terms of the model sensitivity approximation, our proposed POD-aNN is capable of making full use of the auto-differentiation module inside the deep learning frameworks, e.g., $PyTorch$. As a result, both the forward simulation and backward integration with the model-reduced adjoint result in ignorable computational cost. Some numerical experiments were conducted on a synthetic two-dimensional subsurface flow model. The results show that the POD-aNN approach could achieve accurate parameter estimations while the computation cost remains low. The combination of model-order reduction procedure and deep-learning surrogate prompts the simulation-based sensitivity analysis and optimization in various engineering applications, especially for large-scale and computation-intensive physical models. The proposed deep learning hybrid workflow offers an alternative avenue to rapidly achieve comparable simulations and predictions from original high-fidelity physics-based models. The promising application community includes hydrology engineering, underground energy resources (geothermal recovery), and other geoscience domains.