Simulation of multi-phase flow in a subsurface porous media is an essential task for a number of engineering applications including ground water management, contaminant transport, and effective extraction of hydrocarbon resources [2, 3]
. The physics governing subsurface flow simulations are mainly modeled by a system of coupled nonlinear partial differential equations (PDEs) parametrized by subsurface properties (e.g. porosity and permeability). In realistic settings, subsurface models are computationally expensive (i.e. large number of grid block is needed) as the subsurface properties are heterogeneous and the solution exhibit multiscale features [5, 2].
. Monte-Carlo methods are usually employed to propagate the uncertainties in the subsurface properties to the flow response. Monte-Carlo methods are computationally very expensive since a large number of forward simulations are necessary to estimate the statistics of the engineering quantities of interest[2, 3, 6]
. Likewise, Bayesian inference tasks require a very large number of forward simulations to sharpen our knowledge about the unknown model parameters by utilizing field observation data[5, 3]
. For example, Markov-chain Monte-Carlo (MCMC) method (and its variants) requires a large number (in millions) of reservoir simulations to reach convergence and to avoid biased posterior estimates of the model parameters.
Surrogate models can be used to overcome the computational burden of multi-query tasks (e.g. uncertainty quantification, model based optimization) governed by large scale PDEs [7, 8, 9, 10, 11, 12]. Surrogate models are computationally efficient mathematical models that can effectively approximate the main characteristics of the full-order model (full model) 
. A number of surrogate modeling techniques have been developed and could be broadly classified into three classes: simplified physics based models[13, 11], data-fit black-box models [7, 14, 15], and projection based reduced order models commonly referred to as reduced model [16, 17, 18, 19]. Physics based surrogate models are derived from high-fidelity models using approaches such as simplifying physics assumptions, using coarse grids, and/or upscaling of the model parameters [13, 7, 9, 20]. Data-fit models are generated using the detailed simulation data to regress the relation between the input and the corresponding output of interest [7, 15, 21, 22]. For a complete review of various surrogate modeling techniques, we refer the readers to the following papers by Asher et al. , Frangos et al. , Koziel and Leifsson , Razavi et al. .
In projection based reduced order models (utilized in this paper), the governing equations of the full model are projected into a low-dimensional subspace spanned by a small set of basis functions via Galerkin projection [17, 18]. Projection based ROMs relies on the assumption that most of the information and characteristics of the full model state variables can be efficiently represented by linear combinations of only a small number of basis functions. This assumption enables reduced model to accurately capture the input-output relationship of the full model with a significantly lower number of unknowns [7, 17, 18]. Projection based reduced order models are broadly categorized into: system based methods and snapshot based methods. System based methods like balanced truncation realization methods , and Krylov subspace methods  use the characteristics of the full model and have been developed mainly for linear time-invariant problems, although much progress has been done on extensions of these methods to nonlinear problems . Snapshot based methods such as reduced basis methods , proper orthogonal decomposition (POD) [29, 16] derive the projection bases from a set of full model solutions (the snapshots).
In this work, we employ POD based reduced model to accelerate Monte-Carlo simulation of subsurface flow models. The basis functions obtained from the POD is optimal in the sense that, for the same number of basis functions, no other bases can represent the given snapshot set with lower least–squares error than the POD bases [17, 29] (see section 3 for further details). Lumley  was the first to apply POD techniques in fluid flow simulations. Since then, POD procedures has successfully been applied in a number of application areas [e.g., 29, 31, 32, 33, 34, 35, 36].
In the context of fluid flow in porous media, Vermeulen et al.  introduced POD in the confined, groundwater flow problems (linear subsurface flow model). Vermeulen et al.  applied POD in gradient-based optimization problem involving groundwater flow model. McPhee and Yeh  employed POD to enhance the groundwater management optimization problem. Siade et al. 
introduced a new methodology for the optimal selection of snapshots in such a way that the resulting POD basis functions account for the maximal variance of the full model solution. Within the context of oil reservoir simulation,Heijn et al.  and Van Doren et al.  applied POD to accelerate the optimization of a waterflood process. Cardoso et al.  incorporated a new snapshot clustering procedure to enhance the standard POD for oil–water subsurface flow problems.
In the context of Monte-Carlo simulations applied to stochastic subsurface flow problems, POD based ROMs were mainly employed only when the governing equation was linear (or nearly linear) [44, 45, 46, 47]. Pasetto et al.  employed POD based reduced model to construct MC realizations of two dimensional steady state confined groundwater flow subject to a spatially distributed random recharge. Pasetto et al.  applied POD to accelerate the MC simulations of transient confined groundwater flow models with stochastic hydraulic conductivity. Baú  derived a set of POD ROMs for each MC realization of hydraulic conductivity to solve a stochastic, multi-objective, confined groundwater management problem. Boyce and Yeh  applied a single parameter-independent POD reduced model to the deterministic inverse problem and the Bayesian inverse problem involving linear groundwater flow model. In addition to the limitation of using only linear flow models, the UQ tasks in the aforementioned literature involve only low dimensional uncertain parameters.
Within the context of nonlinear subsurface flow problems, the target application of POD was mainly hydrocarbon production optimization, where POD ROMs were used mainly to optimize well control parameters (e.g., bottomhole pressure) [44, 49, 50, 51, 52]. Recently, Jansen and Durlofsky  has done an extensive review on the use of reduced-order models in well control optimization. For the well control applications, POD achieved reasonable levels of accuracy only when the well controls in test runs were relatively close to those used in training runs. In the case where the test controls substantially differ from those used in the initial training runs, additional computational steps were needed. For example refitting the POD basis functions was performed in , which impose some additional computational overhead. Although POD combined with Galerkin projection has been applied more frequently to nonlinear flow problems [33, 16, 51], the effectiveness of POD-Galerkin based model in handling nonlinear systems is limited mainly by two factors. The first factor is related to the treatment of the nonlinear terms in the POD-Galerkin reduced model [53, 54, 44] and the second factor is related to maintaining the overall stability of the resulting reduced model [44, 55, 9, 56, 57].
In relation to computing reduced non-polynomial nonlinear functions, POD based ROMs is usually dependent on the full model state variables and henceforth, the computational cost of evaluating the reduced model is still a function of full model dimension. Several techniques have been developed to reduce the computational cost of evaluating the nonlinear term in POD ROMs including trajectory piecewise linearization (TPWL) , gappy POD technique , missing point estimation (MPE) , best point interpolation method , and discrete empirical interpolation method (DEIM) [59, 53]. Among these techniques, TPWL and DEIM are widely used for efficient treatment of nonlinearities in multi-phase flow reservoir simulations [61, 55, 9].
In TPWL method , the nonlinear function is first approximated by a piecewise linear function obtained by linearizing the full-order model at a predetermined set of points in the time and the parameter space. Then the nonlinear full model is replaced by an adequately weighted sum of the selected linearized systems . Finally, the reduced model can be obtained by projecting the resultant linearized full-order system using standard techniques like POD . The TPWL method was first introduced in  for modeling nonlinear circuits and micromachined devices. In the context of subsurface flow problems, TPWL procedures were applied in [44, 49, 50, 51] to accelerate the solution of production optimization problems.
In DEIM, the nonlinear term in the full model is approximated by a linear combination of a set of basis vectors. The coefficients of expansion are determined by evaluating the nonlinear term only at a small number of selected interpolation points . DEIM was developed in 
for model reduction of general nonlinear system of ordinary differential equations (ODEs) and have been used in several areas[62, 63, 64]. Within the context of subsurface flow problems, Chaturantabut and Sorensen  applied DEIM for model reduction of viscous fingering problems of an incompressible fluid through a two dimensional homogeneous porous medium. Alghareeb and Williams  combined DEIM with POD procedures and the resultant reduced model was applied in waterflood optimization problem. Recently, Ghasemi  applied POD with DEIM to an optimal control problem governed by two-phase flow in a porous media. Next, Ghasemi 
used machine learning technique to construct a number of POD-DEIM local reduced-order models. In that work, machine learning technique was used to construct a number of POD-DEIM local reduced-order models and then a specific local reduced-order model was selected with respect to the current state of the dynamical system during the gradient based optimization task. Similarly,Yoon et al.  used multiple local DEIM approximations in POD reduced model framework to reduce the computational costs of high-fidelity reservoir simulations.
The overall convergence and stability is another issue that limits the applicability of POD-Galerkin based ROMs. POD-Galerkin projection methods manage to decrease the computational complexity by orders of magnitude as a result of state variable’s dimension reduction. However, this reduction goes hand in hand with a loss in accuracy. Moreover, slow convergence and in some cases model instabilities [57, 55, 56] are observed as the errors in the reduced state variables are propagated in time. More specifically, the performance of POD-Galerkin ROMs is directly influenced by the number of POD basis used in the POD-Galerkin projection. However, in many applications involving nonlinear conservation laws (e.g. high Reynold number fluid flow), POD-Galerkin reduced order models have shown poor performance even after retaining a sufficient number of POD basis [57, 29, 16].
Several stabilization techniques have been proposed in the recent literature to build a stabilized POD based reduced models. A notable stabilization technique relies on closing the POD reduced model using a set of closure models similar to those adopted in turbulence modeling [16, 57]. The objective of applying closure models within POD based reduced model is to include the effects of the discarded POD basis functions in the extracted reduced model [16, 57]. Wang et al.  showed that POD-Galerkin reduced model yielded inaccurate and physically implausible results when applied to the numerical simulation of a 3D turbulent flow past a cylinder at Reynold number of . Wang et al.  addressed the aforementioned accuracy and stability issues of POD reduced model by various closure models, where artificial viscosity was added to the real viscosity parameter to stabilize the POD based reduced model.
Another major approach to enhance the stability of POD-Galerkin reduced model is to compute a new set of optimal basis or to improve the POD basis vectors by solving a constrained optimization problem. Bui-Thanh et al.  determined a new set of optimal basis vectors by formulating an optimization problem constrained by the equations of the resultant reduced model and demonstrated the stability of the proposed approach on linear dynamical systems. We note that POD-Galerkin reduced model orthogonally projects the nonlinear residual into the subspace spanned by the POD basis vectors. Unlike POD-Galerkin reduced model, Petrov-Galerkin projection scheme design a different set of orthonormal basis called left reduced order basis into which the nonlinear residual is projected. Carlberg et al.  formulated stable Petrov-Galerkin reduced model in which the left reduced order basis vectors were computed from an optimization problem at every iteration of the Gauss Newton method. He  observed that poor spectral properties of the reduced Jacobian matrix could cause numerical instabilities in POD-Galerkin TPWL reduced model. Hence, He  improved the stability of the POD based reduced model by determining the optimal dimension of the reduced model through an extensive search over a range of integer numbers. We note that all the above mentioned optimization procedures involve computationally expensive procedures to maintain stability and in many cases, the stability of the extracted reduced model is still not guaranteed [55, 9].
Recently, data-fit black-box models have been combined with POD  to develop non-intrusive POD based ROMs, where the data-fit models are used to regress the relationship between the input parameter and the reduced representation of the full model state vector. Hence, non-intrusive ROMs do not require any knowledge of the full-order model and is mainly developed to circumvent the shortcomings in accessing the governing equations of the full model . However, it can also be used to address the stability and nonlinearity issues of POD based ROMs. Wang et al.  developed a non-intrusive POD reduced model using Recurrent Neural Network (RNN) as a data-fit model and presented two fluid dynamics test cases namely, flow past a cylinder and a simplified wind driven ocean gyre. RNN is a class of artificial neural network [71, 72] which has at least one feedback connection in addition to the feedforward connections [73, 71, 74]. In the context of data-fit models, RNN has been successfully applied to various sequence modeling tasks such as automatic speech recognition and system identification of time series data [75, 76, 77, 78]. Additionally, RNN has been applied to emulate the evolution of nonlinear dynamical systems in a number of applications [79, 80] and henceforth has large potential in building reduced order models. However, the applicability of non-intrusive ROMs is severely undermined in many real-world problems, where increasing the dimensionality of the input parameter space increases the complexity and training time of the data-fit model.
In summary, among many surrogate modeling techniques, POD-Galerkin reduced model is a viable option for accelerating multi-query tasks like UQ. Generally, POD-Galerkin reduced model is well established for linear systems and for nonlinear systems with parametric dependence, POD could be either combined with TPWL or with DEIM for modeling subsurface flow systems [44, 49, 50, 61]. However, POD reduced model does not preserve the stability properties of the corresponding full order model and current state of the art POD stabilization techniques [57, 55, 9] are not cost effective and ultimately do not guarantee stability of the extracted reduced order models.
In this paper, we use DR-RNN  to alleviate the potential limitations of POD-Galerkin reduced models. More specifically, we combine DR-RNN with POD-Galerkin and DEIM methods to derive an accurate and computationally effective reduced model for uncertainty quantification (UQ) tasks. The architecture of DR-RNN is inspired by the iterative line search methods where the parameters of the DR-RNN are optimized such that the residual of the numerically discretized PDEs is minimized [81, 82, 1]. Unlike the standard RNN which is very generic, DR-RNN  uses the residual of the discretized differential equation. In addition, the parameters of the DR-RNN are fitted such that the computed DR-RNN output optimally minimizes the residual of the targeted equation. In this context, DR-RNN is a physics aware RNN as it is tailored to leverage the physics embedded in the targeted dynamical system (i.e. residual of the equation or reduced residual in the current manuscript).
The resultant reduced model obtained from DR-RNN combined with POD-Galerkin and DEIM algorithm has a number of salient features. First, the dynamics of DR-RNN is explicit in time with superior convergence and stability properties for large time steps that violate the numerical stability conditions [1, 83]. Second, as the dynamics modeled in DR-RNN are explicit in time, there is a reduction in the computational complexity of the extracted reduced model from corresponding to implicit POD-DEIM reduced order models, to , where is the size of the reduced model. Third, DR-RNN requires only very few training samples (obtained by solving the full model) to optimize the parameters of the DR-RNN as it accounts for the physics of the full model within the RNN architecture (via the reduced residual). This is a major advantage when compared to pure data-driven algorithms (e.g. standard RNN architectures). Moreover, DR-RNN can effectively emulate the parameterized nonlinear dynamical system with a significantly lower number of parameters in comparison to standard RNN architectures .
In this work, we demonstrate the superior properties of DR-RNN in accelerating UQ tasks for subsurface reservoir models using Monte-Carlo method. As far as we are aware, the use of a single parameter-independent POD-Galerkin reduced model in Monte-Carlo method involving nonlinear subsurface flow with high dimensional stochastic permeability field has not been previously explored. The reason is that the resultant reduced model might require significantly more basis functions to reconstruct stable solutions [44, 49, 47, 61]. However, only a single set of small number of POD basis functions would be sufficient to reconstruct the solution with reasonable accuracy using least–squares (see section 3.2 for more details). Hence, the aim of this paper is to illustrate how DR-RNN could be used to reconstruct stable solutions emulating the full model dynamics using only a small set of POD basis functions. The proposed DR-RNN technique is validated on two forward uncertainty quantification problems involving two-phase flow in subsurface porous media. The two flow problems are commonly known within the reservoir simulation community as the quarter five spot problem and the uniform flow problem 
. In these two numerical examples, the permeability field is modeled as log-normal distribution. The obtained results demonstrate that DR-RNN combined with POD-DEIM provides an accurate and stable reduced order model with a drastic reduction in the computational cost. The reason for selecting simplified flow problems is to illustrate the potential benefit of DR-RNN to formulate an accurate and computationally effective POD-DEIM reduced model for flow problems where the standard POD-Galerkin reduced models are inaccurate and possibly unstable. We also note that DR-RNN architecture is generic and could be used to emulate any well-posed nonlinear dynamical system including subsurface flow problems while accounting for capillary pressure effects, gravity effects and compressibility.
The outline of the rest of this manuscript is as follows: In section 2, we present the formulation of multi-phase flow problem in a porous media. In section 3, we introduce POD-Galerkin method for model reduction followed by a discussion of DEIM for handling nonlinear systems. In Section 4, we describe the architecture of DR-RNN and in section 5, we evaluate the reduced model derived by combining DR-RNN with POD-DEIM on two uncertainty quantification test cases. Finally, in Section 6, we present the conclusions of this manuscript.
2 Problem Formulation
The equations governing two-phase flow of a wetting phase (water) and non-wetting phase (e.g. oil) in a porous media are the conservation of mass (continuity) equation and Darcy’s law for each phase [4, 9, 84, 85]. The continuity equation for each phase takes the form
where the subscript denotes the water phase, the subscript denotes the oil phase,
is the absolute permeability tensor,is the phase mobility, with the relative permeability to phase and the viscosity of phase , is the phase pressure, is the density of phase , is the gravitational acceleration, is the depth, is the porosity, is the saturation of the phase and is the phase source and sink terms [4, 84]. Further, the phase saturations are constrained by , since the oil and the water jointly fill the void space [4, 9].
The phase velocities are modeled by the multiphase Darcy’s law to relate the phase velocities to the phase pressures and takes the form
where is the phase velocity. The phase relative permeabilities and the capillary pressure () are usually modeled as functions of the phase saturations . Neglecting the capillary pressure, the compressibility effects, the gravitational effects, and assuming the density ratio to be equal to one, the continuity equations (Eq. (1)) can be combined with the Darcy’s law (Eq. (2)) to derive a global pressure equation and the saturation equation for water phase [4, 9, 85]. The simplified global pressure equation takes the form
where is the global pressure, is the total mobility, is the source and sink term. The saturation equation for the water phase takes the following form
where is a function of saturation termed as the fractional flow function for the water phase, is the total velocity vector and is the water saturation [4, 84]. In the rest of the paper, we write the water phase saturation as for simplicity. The coupled equations Eq. (3) and Eq. (4) could then be solved for the evolution of the saturation by providing the appropriate initial and boundary conditions. Equation (3) and Eq. (4) are continuous (in space and time) form of the full model.
The discrete form of the full model is obtained by dividing the problem domain into grid blocks and then applying the finite volume method to discretize the spatial derivatives of Eq. (3) and Eq. (4). The discretized pressure equation takes the form
where , , and is the pressure vector in which each component of represent the pressure value at the th grid block. Similarly, the spatially discretized saturation equation takes the form
where , , is the total velocity vector, and is the saturation vector in which each component of is the saturation value at the th grid block.
Eq. (5) and Eq. (6) are the discrete form of the full model for multi-phase flow problem under consideration. These two equations exhibit two way coupling from the dependence of the matrix on the mobilities in the pressure full model (Eq. (5)) and from the dependence of the matrix on the velocity vector in the saturation full model (Eq. (6)). In this paper, we adopt an implicit sequential splitting method to solve the full model (Eq. (5) and Eq. (6)). In this method, the saturation vector from the present time step is used to assemble the matrix in Eq. (5) and then the pressure full model (Eq. (5)) is solved for the pressure vector . Following that, the velocity vector (computed from the pressure vector ) is used to assemble the matrix in Eq. (6) and then the saturation full model (Eq. (6)) is solved implicitly in time for the saturation at the next time step. In the following section, we formulate a Galerkin projection based reduced model to reduce the computational effort for multi-query tasks (e.g. uncertainty quantification) involving repeated solutions of Eq. (5) and Eq. (6), when (the number of grid block) is large [53, 61].
3 Reduced Order Model Formulation
In this section, we formulate the POD-Galerkin reduced model (POD reduced model) and POD-DEIM reduced model where POD-Galerkin is combined with DEIM for handling the nonlinear terms. Both methods are introduced to reduce the computational effort associated with solving the full model (Eq. (5) and Eq. (6)).
3.1 POD basis
As stated in section 1, POD based reduced model is a projection based reduced order model in which the governing equations are projected onto an optimal low-dimensional subspace spanned by a small set of basis vectors. Galerkin projection reduced model is based on the assumption that most of the system information and characteristics can be efficiently represented by linear combinations of only a small number of basis vectors .
The optimal basis vectors
in POD are computed by singular value decomposition (SVD) of the solution snapshot matrix. The solution snapshot matrix is obtained from a set of solution vectors of size obtained by solving the full model at selected points in the input parameter space. The SVD of is expressed as
where, , is the left singular matrix and is the diagonal matrix containing the singular values of the snapshot matrix in descending order. The dominant left singular vectors corresponding to the first largest singular values represents the basis vectors to span the optimal subspace of POD based reduced model. Thus, the first step in deriving the POD based reduced model is to express the state vector of the full-order model by a linear combination of basis vectors as following
where is the reduced state vector representation of full dimensional state vector , and is the matrix that contains orthonormal basis vectors in its columns.
By following this step, for example, the optimal basis vectors for the saturation state vector are obtained from the SVD of the saturation snapshot matrix , where is the number of time steps and is the number of samples of input parameter used to build the snapshot matrix. The SVD of is expressed as
where is the left singular matrix, is the diagonal matrix containing the singular values of the snapshot matrix in descending order. The saturation state vector is optimally expressed as
where is the reduced state vector representation of , is the matrix that contains orthonormal basis vectors in its columns. Similarly, we can represent the pressure state vector from its reduced state vector representation using optimal basis matrix obtained from the SVD of the pressure snapshot matrix .
3.2 Least–squares approximation
The capacity of a set of basis functions to represent a new solution vector could be tested using least–squares fitting [86, 87]. For example, the least–squares solution for approximating a saturation state vector is defined as
The associated error termed as least–squares errors in approximating by using only basis vectors is given by
The least–squares error (Eq. (12)) is equivalent to the omitted energy [88, 16]. In practice, is commonly chosen as the smallest integer such that the relative omitted energy is less than a preset value (e.g. ), where the omitted energy is defined by the following equation
Similar expressions mentioned in Eq. (11), Eq. (12), and Eq. (13) can be obtained for the pressure state vector as well. We note that least–squares errors are not necessarily equivalent to the omitted energy for state vectors not included in the snapshot matrix or for the state vector solved at a new point in the input parameter space as these new vectors might not fall within the span of the snapshot matrix [7, 88]. The least–squares solution is the best approximation of the state variables in the sense that, for the chosen low dimensional subspace , no other low dimensional approximation can represent the given snapshot set with a lower least–squares error [17, 29, 16]. In this paper, we use the best approximation of the state variables to assess the quality of the approximation obtained from different reduced order models in the numerical examples presented in section 5.
Once the POD basis vectors are obtained, the reduced representation of the pressure vector is substituted into the pressure full model (Eq. (5)), followed by Galerkin projection of the pressure equation into the subspace spanned by . The resulting POD based reduced model for the pressure equation then takes the following form
where and . Similarly, POD based reduced model for the saturation equation (Eq. (6)) takes the form
where and .
The POD based reduced model formulated by Eq. (14) and Eq. (15) is of the reduced dimension . However, the nonlinear function in Eq. (15) is still of the order of full dimension . Moreover, the reduced Jacobian matrix needed for Newton like iterations to solve this nonlinear equation is also of order  as it relies on evaluating the full order nonlinear function . Therefore, for problems with general nonlinear functions involved in POD based reduced model, the computational cost of solving the reduced system is still a function of the full system dimension .
Discrete Empirical Interpolation Method (DEIM) was introduced in  to approximate the nonlinear terms in POD based reduced model using a limited number of points that are independent of the full system dimension . Similar to POD, the first step of DEIM is to approximate the nonlinear function in Eq. (15) using a separate set of basis vectors as
where is the coefficient of expansion of the nonlinear function in the reduced subspace spanned by , is the matrix containing the first columns of the left singular matrix obtained from the SVD of the the snapshot matrix of the nonlinear function . We note that no additional computational costs are associated with collecting the snapshot matrix of the nonlinear terms as it is already evaluated during the computation of the state snapshot vectors. The nonlinear term in Eq. (15) can then be expressed as
The matrix factor in Eq. (17) is precomputed before solving Eq. (15). The overdetermined system is approximated using the DEIM algorithm introduced in  by first computing a matrix that selects rows of the matrix to obtain as following
Using this expression of to approximate the nonlinear function in Eq. (17), we obtain a nonlinear term that is independent of that takes the form
where the matrix termed as the DEIM matrix. Similarly, the Jacobian of the nonlinear term in Eq. (15) is approximated using DEIM as following
We note that POD-DEIM formulation is independent of the full model dimension and that the DEIM procedure exploits the structure of the nonlinear function as component-wise operation at .
4 Deep Residual RNN
POD-DEIM reduced order models, as introduced in the last chapter, could be used to perform parametric UQ tasks. However, the POD-DEIM formulation is nonlinear and relies on using Newton method at each time step to solve the resulting system of nonlinear equations. The computational efficiency of the Newton iteration depends on the method employed to assemble the Jacobian matrix and more importantly on the conditioning of the reduced Jacobian matrix. It also depends on the method used to solve the resulting linear system at each iteration of the Newton step and generally, it takes operations for each saturation update [1, 81]. Moreover, previous studies [9, 55] pointed to the loss of stability of POD-Galerkin reduced model in several cases and it was attributed to ill-conditioning and poor spectral properties of the reduced Jacobian matrix.
In this paper, we build on the recently introduced DR-RNN  and formulate an accurate POD-DEIM reduced order models. DR-RNN is a deep RNN architecture , constructed by stacking physics aware network layers. DR-RNN could be applied to any nonlinear dynamical system of the form
where is the state variable at time , is a system parameter vector, the matrix is the linear part of the dynamical system and the vector is the nonlinear term . The state variable at different time steps is obtained by solving the nonlinear residual equation defined as
where is termed as the residual vector at time step and is the approximate solution of Eq. (22) at time step obtained by using implicit Euler time integration method. DR-RNN  approximates the solution of Eq. (22) using the following iterative update equations
where are the training parameters of DR-RNN, is the activation function, is an element-wise multiplication operator, is the residual in layer obtained by substituting into Eq. (23) and is an exponentially decaying squared norm of the residual defined by
where are fraction factors and is a smoothing term to avoid divisions by zero . In this formulation, we set
. The architecture of DR-RNN is inspired by the rmsprop algorithm which is a variant of the steepest descent method. The DR-RNN output at each time step is defined as
The formulation of DR-RNN is explicit in time and has a fixed number of iterations per time step. However, the dimension of the DR-RNN system depends on the dimension of the residual. For example, DR-RNN (Eq. (24)) can be derived from the POD-DEIM reduced model residual (). In such setting, the DR-RNN dynamics has a fixed computational budget of for each time step. In addition, DR-RNN has the prospect of employing large time step violating the numerical stability constraint . Furthermore, DR-RNN does not rely on the reduced Jacobian matrix to approximate the solution of POD-DEIM reduced model.
The DR-RNN parameters are fitted by minimizing the mean square error (mse) defined by
is commonly estimated by a technique called Backpropagation Through Time (BPTT)[89, 90, 71, 72]
, which backpropagates the gradient of the loss functionwith respect to in time over the length of the simulation.
5 Numerical Experiments
In this section, we evaluate the performance of the reduced order models based on DR-RNN against the standard implementation of POD-Galerkin reduced model. Specifically, we develop two POD-Galerkin based reduced model using DR-RNN architecture namely, DR-RNN (DR-RNN combined with POD-Galerkin) and DR-RNN (DR-RNN combined with POD-Galerkin and DEIM). The numerical evaluations are performed using two uncertainty quantification tasks involving subsurface flow models. We did not include standard POD-DEIM reduced model implementation as we expect that the standard POD reduced model results to be far superior [53, 1, 53].
The outline of this section is as follows: In subsection 5.1, we present the description of the flow problem, followed by a brief description of the finite-volume approach employed for obtaining the full-order model solution. Following that, in subsection 5.2, we outline the specific details to formulate POD reduced model. Then, we list the settings adopted to model the DR-RNN ROMs (i.e. number of layers, optimization settings, etc) in the subsection 5.3. In subsection 5.4, we provide a set of error metrics utilized to evaluate the performance of the different ROMs. In subsection 5.5, we present the numerical results for the quarter five spot model followed by results for the uniform flow model in the subsection 5.6.
5.1 Full-order model setup
We consider a two-phase (oil and water) porous media flow problem over the two-dimensional domain meters. The equations governing the two-phase flow are the pressure equation (Eq. (3)) and the saturation equation (Eq. (4)). The relative permeability is defined as a function of saturation using Corey’s model , , where , is the irreducible water saturation and is the residual oil saturation . We set and . We set the initial water saturation over the domain to the irreducible water saturation . The water and oil viscosities are 1 and 1.5 centipoise, respectively. The porosity is assumed to be a constant value of 0.2 over the entire problem domain. The uncertain permeability field is modeled as a log-normal distribution function with zero mean and exponential covariance kernel of the form
where is the variance, is the correlation length. In all test cases, we set to 1 and the correlation length to . Figure 1 shows several realizations of the log-permeability values.
For solving the full-order model, the problem domain is discretized using a uniform grid of blocks. The pressure equation is discretized using simple finite volume method (aka. Two Point Flux Approximation)  and an upwind finite-volume scheme is used to discretized the saturation equation. For the time discretization, an implicit backward Euler method combined with Newton-Raphson iterative method is used to solve the resulting system of nonlinear equations. We set the time step size to and the total number of time steps is set to . We note that, the time is measured in a non-dimensional quantity called pore volumes injected (PVI). PVI defines the net volume of water injected as a fraction of the total pore volume. As the pressure changes at much slower rate than the saturation, the pressure equation (and hence the velocity) is solved at every 8th saturation time step. For reference solutions, this system of equations is solved for random permeability realizations to estimate an ensemble based statistics using Monte-Carlo method .
5.2 POD-Galerkin based reduced model setup
The first step in formulating POD reduced model is to compute the optimal POD basis matrices and . In order to obtain these basis matrices, we initially preformed a realization clustering algorithm to enforce the diversity of the collected snapshots and clustered the random permeability realizations into clusters . Then, we randomly selected a single permeability realization from each cluster (total 45 random samples of the permeability field). The full system is then solved for each of the 45 realizations and the solution vectors are collected to build the snapshot matrices (pressure, saturation, nonlinear function). Finally, we compute the POD basis matrices from the SVD of the collected snapshot matrices.
Following that, the obtained basis vectors are used to build POD reduced model (as detailed in the section 3). We then employ the same sequential implicit technique settings adopted for obtaining the full model solutions to solve the resultant POD reduced model. For numerical evaluations, we solve the POD reduced model for the same permeability realizations to estimate an ensemble based statistics in the engineering quantities of interest.
5.3 DR-RNN setup
In all the numerical test cases, we utilize DR-RNN with six layers ( in Eq. (24)). We evaluate DR-RNN and DR-RNN for different number of POD basis, however, we fix the number of DEIM basis to . The PyTorch framework 
, a deep learning python package usingTorch library as a backend, is used to implement the DR-RNN. Further, we optimize the DR-RNN model parameters using rmsprop algorithm [82, 91] as implemented in PyTorch, where we set the weighted average parameter to and the learning rate to . The weight matrix in Eq. (24
) is initialized randomly from the uniform distribution function. The vector training parameter in Eq. (24) is initialized randomly from the uniform distribution function . The scalar training parameters in Eq. (24) are initialized randomly from the uniform distribution
. We set the hyperparametersand in Eq. (25) to and , respectively. The formulated DR-RNN and DR-RNN are trained to approximate the reduced state vector representation obtained from least–squares fits. Specifically, we collect a set of best reduced state vector representation of the saturation state vector using . The collected set of reduced state vectors is then used to train the parameters of the DR-RNN by minimizing the loss function defined in Eq. (27).
5.4 Evaluation metrics
We evaluate the performance of DR-RNN and DR-RNN using two time specific error metrics defined by
where is the realization index, and is computed from the reduced model. Additionally, we utilize two relative error metrics defined as
where all the time snapshots of saturation vectors in all realizations are used.
5.5 Numerical test case 1
In this test case, water is injected at the lower left corner () of the domain and a mixture of oil and water is produced at the top right corner of the domain (). We set the injection rate at () and at () as defined in Eq. (4). We impose a no flow boundary condition in all the four sides of the domain. We fix the number of pressure POD basis to 5 and obtain all the ROMs for a set of different number of saturation POD basis functions (). The configuration of the problem domain is shown in top left panel of Figure 2, where the blue spot in the lower left corner (0,0) corresponds to the injector well and the blue spot in the upper right corner (1,1) corresponds to the production well. Figure 2 shows the singular values of the pressure snapshot matrix in the top right panel, the saturation snapshot matrix in the bottom left panel, and the nonlinear function snapshot matrix in the bottom right panel.
The mean water saturation plots over the simulation time are shown in Figure 3, where the results in the top row corresponds to using POD basis and the results in the bottom row corresponds to using POD basis. The sub-plots in Figure 3 are arranged from left to right following the numbering of the spatial points shown in Figure 2. From these results, it is clear that DR-DR-RNN and DR-RNN results are very close to the least–square solutions (LS fit). In Figure 3, POD-Galerkin reduced model yields extremely inaccurate and unstable results. We attribute the small errors in DR-RNN and DR-RNN results to the insufficient number of POD basis vectors and we note that the error magnitude is equivalent to the optimal values obtained by least-squares projection.
show the results for the first (mean) and second (standard deviation) moments of the saturation field at timePVI obtained from the full model and from the various ROMs. In these figures (4, 5, and 6), results for 10 POD basis are shown in the top row and results for 20 POD basis are shown in the bottom row. As shown in Figure 4, the mean saturation obtained from DR-RNN ROMs are almost indistinguishable from the reference results. However, the mean saturation field obtained from POD reduced model (left panels of Figure 6) deviates significantly from the reference mean saturation.
In Figure 5, we observe small discrepancy of standard deviation results obtained the DR-RNN ROMs in comparison to the full model results especially near the location of the mean water saturation front. Figure 6 (right panels) shows the standard deviation results obtained by POD reduced model which show significant inaccuracies that could be indicative to instabilities of the obtained solutions. We note that the white spots in Figure 6 correspond to out of limits shown in colorbar.
Figure 7 compares the saturation PDF estimated from the ensemble of numerical solutions (ROMs and the full model). Figure 7 settings are similar to the one adopted in Figure 3. In Figure 7, we can see that all the plots obtained from DR-DR-RNN and DR-RNN are indistinguishable from the plots obtained from the LS fit (the best approximation). Further, we observe that the saturation PDF obtained from DR-DR-RNN and DR-RNN follow nearly the same trend of saturation PDF obtained from the full model when the reference distribution is unimodal. However, we observe some discrepancy when the distributions are multimodal. Please note that similar discrepancy is also observed in the PDF obtained from LS fit. Hence, we postulate that these discrepancies are attributed to the limited number of POD basis vectors utilized. In Figure 7, POD reduced model yields very inaccurate approximation of the saturation PDF irrespective of the number of POD basis.
Figures 8 and 9 displays samples of and errors at time PVI obtained from all the ROMs. All the ROMs use 10 POD basis to display the errors in Figure 8 and likewise 20 POD basis to display the errors in Figure 9. From these figures, we can see that the POD reduced model approximation errors are at least an order of magnitude more than the least–squares solution errors (Eq. (11)), whereas the errors obtained from DR-RNN and DR-RNN are nearly indistinguishable from the least–squares projection errors.
|Error||#Basis||Reduced Order Models|
We further list in Table 1, the and errors for the saturation field. From Table 1, we can see that the approximation errors obtained from DR-RNN and DR-RNN have the same order of magnitude as the least–squares (best approximation) errors. Further, in Table 1, the approximation errors obtained from all ROMs except POD reduced model decreases when we increase the number of POD basis. These results conform with the decay of singular values of the saturation snapshot matrix. In Table 1, the approximation errors obtained from POD reduced model are at least an order of magnitude larger than other methods. Also, we observe that POD reduced model results might be worst when we include more basis functions. These results conform with the results presented in , where it was shown that selecting large number of basis vectors based on singular values may not lead to stable POD-Galerkin reduced model. Further, it was presented in  that the relation between the stability property of POD-Galerkin reduced model and the number of basis vectors used in POD-Galerkin projection is somewhat random and that the use of more POD basis vectors do not necessarily lead to improved stability.
5.6 Numerical test case 2
In this test case, the boundary conditions are set to no flow boundary conditions on the two sides aligned in the horizontal direction (top and bottom). Water is injected from the left side of the domain boundary and fluids are produced from the right side boundary of the domain. The total inflow rate from the left side is set to 0.05 and the total outflow rate from the right side to 0.05 as the problem is incompressible. Similar to test case 1, we evaluate all the ROMs for two different number of saturation POD basis functions (). Also, we fix the number of POD basis for the pressure state vector to 5. Figure 10 shows the setup for test case 2 and the corresponding singular values of the snapshot matrices , , and .
Figure 11 shows the time plot of mean water saturation obtained from all the ROMs and from the full model. The display settings in Figure 11 are the same as defined in Figure 3. In Figure 11, we can see that all the results obtained from DR-RNN, DR-RNN, and the LS fit (the best approximation) closely approximates the full model whereas POD reduced model yields extremely inaccurate results regardless of the number of utilized POD basis.
Figures 12, 13, and 14 shows the results for the mean and standard deviation of the saturation field at PVI. From these figures, we can conclude that all the plots obtained from DR-RNN ROMs are almost indistinguishable from the LS fit (the best approximation) results, whereas the plots obtained from POD reduced model (Figure 14) exhibit significant discrepancy when compared to the plots shown in Figure 12. Again, we note that the white spots displayed in Figure 14 are the regions whose values are out of the limits marked in the respective colorbar.
Figure 15 compares the saturation PDF estimated from the ensemble of numerical solutions obtained from all the ROMs and the full model. The plotted results show that DR-RNN, DR-RNN predictions are nearly indistinguishable from the plots obtained from the full model and are very close to the best possible approximation using LS fit. Further, Figure 15 shows that all the saturation PDFs obtained from full model are uni-modal distribution. Similar to test case 1, POD reduced model yields inaccurate approximation of the saturation PDFs.
|Error||#Basis||Reduced Order Models|
We further list in Table 2, the error metrics and for the saturation fields. From Table 2, we can see that the approximation errors obtained from DR-RNN and DR-RNN are almost close to the least–squares (best approximation) approximation errors. However, the POD reduced model yields very inaccurate results due to numerical instabilities.
In this work, we extended the DR-RNN introduced in  into nonlinear multi-phase flow problem with distributed uncertain parameters. In this extended formulation, DR-RNN based on the reduced residual obtained from POD-DEIM reduced model is used to construct the reduced order model termed DR-RNN. We evaluated the proposed DR-RNN on two forward uncertainty quantification problems involving two-phase flow in subsurface porous media. The uncertainty parameter is the permeability field modeled as log-normal distribution. In the two test cases, full order model and ROMs are solved for random permeability realizations to estimate an ensemble based statistics using Monte-Carlo method. full model and POD reduced model used implicit time stepping method as the time step size violates the numerical stability condition. However, DR-RNN architecture employs explicit time stepping procedure for the same step size used in full model and POD reduced model. Hence, DR-RNN had a limited computational complexity instead of per saturation update, where is the dimension of the POD reduced model, is the number of stacked network layers in DR-RNN and is the average number of Newton iterations used in the standard POD-DEIM reduced model. The obtained numerical results shows that DR-RNN provides accurate and stable approximations of the full model in comparison to the standard POD reduced model.
Future work should consider the development of accurate and stable DR-RNN for UQ tasks involving subsurface flow simulations with the additional effects including the capillary pressure, compressibility, and the gravitational effects. In addition, it will be of interest to explore the applicability of DR-RNN for UQ tasks with the permeability fields that has randomly oriented channels or barriers. The use of DR-RNN for history matching [5, 3], where we minimize the mismatch between simulated and field observation data by adjusting the geological model parameters is also expected to show significant reduction of the computational cost.
- Nagoor Kani and Elsheikh  J. Nagoor Kani and A. H. Elsheikh. DR-RNN: A deep residual recurrent neural network for model reduction. ArXiv e-prints, September 2017.
- Petvipusit et al.  Kurt R. Petvipusit, Ahmed H. Elsheikh, Tara C. Laforce, Peter R. King, and Martin J. Blunt. Robust optimisation of CO2 sequestration strategies under geological uncertainty using adaptive sparse grid surrogates. Computational Geosciences, 18(5):763–778, Oct 2014. ISSN 1573-1499. URL https://doi.org/10.1007/s10596-014-9425-z.
- Elsheikh et al.  Ahmed H. Elsheikh, Mary F. Wheeler, and Ibrahim Hoteit. Nested sampling algorithm for subsurface flow model selection, uncertainty quantification, and nonlinear calibration. Water Resources Research, 49(12):8383–8399, 2013. ISSN 1944-7973. URL http://dx.doi.org/10.1002/2012WR013406.
- Aarnes et al.  Jørg E Aarnes, Tore Gimse, and Knut A Lie. An introduction to the numerics of flow in porous media using Matlab. In Geometric modelling, numerical simulation, and optimization, pages 265–306. Springer, 2007.
- Elsheikh et al.  Ahmed H. Elsheikh, Matthew Jackson, and Tara Laforce. Bayesian reservoir history matching considering model and parameter uncertainties. Mathematical Geosciences, 44(5):515–543, Jul 2012. ISSN 1874-8953. URL https://doi.org/10.1007/s11004-012-9397-2.
- Ibrahima  Fayadhoi Ibrahima. Probability distribution methods for nonlinear transport in heterogenous porous media. Ph.D. thesis, Stanford Univeristy, 2016.
- Frangos et al.  Michalis Frangos, Youssef Marzouk, Karen Willcox, and Bart VB Waanders. Surrogate and Reduced-Order Modeling: A Comparison of Approaches for Large-Scale Statistical Inverse Problems, pages 123–149. John Wiley & Sons, Ltd, 2010. ISBN 9780470685853. URL http://dx.doi.org/10.1002/9780470685853.ch7.
- Koziel and Leifsson  Slawomir Koziel and Leifur Leifsson. Surrogate-based modeling and optimization. Applications in Engineering, 2013.
- He  Jincong He. Reduced-order modeling for oil-water and compositional systems, with application to data assimilation and production optimization. Ph.D. thesis, Stanford Univeristy, 2013.
- Elsheikh et al.  Ahmed H. Elsheikh, Ibrahim Hoteit, and Mary F. Wheeler. Efficient bayesian inference of subsurface flow models using nested sampling and sparse polynomial chaos surrogates. Computer Methods in Applied Mechanics and Engineering, 269:515–537, 2014. URL http://dx.doi.org/10.1016/j.cma.2013.11.001.
- Josset et al.  Laureline Josset, Vasily Demyanov, Ahmed H. Elsheikh, and Ivan Lunati. Accelerating Monte Carlo Markov chains with proxy and error models. Computers & Geosciences, 85:38–48, 2015. URL https://doi.org/10.1016/j.cageo.2015.07.003.
- Bazargan et al.  Hamid Bazargan, Mike Christie, Ahmed H. Elsheikh, and Mohammad Ahmadi. Surrogate accelerated sampling of reservoir models with complex structures using sparse polynomial chaos expansion. Advances in Water Resources, 86:385–399, 2015. URL http://dx.doi.org/10.1016/j.advwatres.2015.09.009.
- Durlofsky and Chen  Louis J Durlofsky and Yuguang Chen. Uncertainty quantification for subsurface flow problems using coarse-scale models. In Numerical analysis of multiscale problems, pages 163–202. Springer, 2012.
- Li et al.  Hao Li, Zhien Zhang, and Zhijian Liu. Application of artificial neural networks for catalysis: a review. Catalysts, 7(10):306, 2017.
- Yeten et al.  Burak Yeten, Alexandre Castellini, Baris Guyaguler, and WH Chen. A comparison study on experimental design and response surface methodologies. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 2005.
- Berkooz et al.  Gal Berkooz, Philip Holmes, and John L Lumley. The proper orthogonal decomposition in the analysis of turbulent flows. Annual review of fluid mechanics, 25(1):539–575, 1993.
- Lassila et al.  Toni Lassila, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza. Model order reduction in fluid dynamics: challenges and perspectives. In Reduced Order Methods for modeling and computational reduction, pages 235–273. Springer, 2014.
- Antoulas et al.  Athanasios C Antoulas, Danny C Sorensen, and Serkan Gugercin. A survey of model reduction methods for large-scale systems. Contemporary mathematics, 280:193–220, 2001.
- Fang et al.  F. Fang, C.C. Pain, I.M. Navon, A.H. Elsheikh, J. Du, and D. Xiao. Non-linear petrov–galerkin methods for reduced order hyperbolic equations and discontinuous finite element methods. Journal of Computational Physics, 234:540––559, 2013. URL http://dx.doi.org/10.1016/j.jcp.2012.10.011.
- Babaei et al.  Masoud Babaei, Ahmed H. Elsheikh, and Peter R. King. A comparison study between an adaptive quadtree grid and uniform grid upscaling for reservoir simulation. Transport in Porous Media, 98:377–400, 2013. URL http://dx.doi.org/10.1007/s11242-013-0149-7.
- Abdi-Khanghah et al.  Mahdi Abdi-Khanghah, Amin Bemani, Zahra Naserzadeh, and Zhien Zhang. Prediction of solubility of n-alkanes in supercritical co 2 using rbf-ann and mlp-ann. Journal of CO2 Utilization, 25:108–119, 2018.
- A.Wood  David A.Wood. transparent open-box learning network provides insight to complex systems and a performance benchmark for more-opaque machine learning algorithms. Advances in Geo-Energy Research, 2(2):148–162, 2018.
- Asher et al.  Michael J Asher, Barry FW Croke, Anthony J Jakeman, and Luk JM Peeters. A review of surrogate models and their application to groundwater modeling. Water Resources Research, 51(8):5957–5973, 2015.
- Razavi et al.  Saman Razavi, Bryan A Tolson, and Donald H Burn. Review of surrogate modeling in water resources. Water Resources Research, 48(7), 2012.
- Gugercin and Antoulas  Serkan Gugercin and Athanasios C Antoulas. A survey of model reduction by balanced truncation and some new results. International Journal of Control, 77(8):748–766, 2004.
- Freund  Roland W Freund. Model reduction methods based on Krylov subspaces. Acta Numerica, 12:267–319, 2003.
- Lall et al.  Sanjay Lall, Jerrold E Marsden, and Sonja Glavaški. A subspace approach to balanced truncation for model reduction of nonlinear control systems. International journal of robust and nonlinear control, 12(6):519–535, 2002.
- Rozza et al.  Gianluigi Rozza, Dinh BP Huynh, and Anthony T Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Archives of Computational Methods in Engineering, 15(3), 2007. URL https://doi.org/10.1007/BF03024948.
- Sirovich  Lawrence Sirovich. Turbulence and the dynamics of coherent structures. i. coherent structures. Quarterly of applied mathematics, 45(3):561–571, 1987.
- Lumley  John L Lumley. The structure of inhomogeneous turbulent flows. Atmospheric turbulence and radio wave propagation, 1967.
- Zheng et al.  Daguang Zheng, Karlene A Hoo, and Michael J Piovoso. Low-order model identification of distributed parameter systems by a combination of singular value decomposition and the karhunen- loève expansion. Industrial & Engineering Chemistry Research, 41(6):1545–1556, 2002.
- Cao et al.  Yanhua Cao, Jiang Zhu, Zhendong Luo, and Ionel M Navon. Reduced-order modeling of the upper tropical pacific ocean model using proper orthogonal decomposition. Computers & Mathematics with Applications, 52(8-9):1373–1386, 2006.
- Bui-Thanh et al.  Tan Bui-Thanh, Murali Damodaran, and Karen E Willcox. Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition. AIAA journal, 42(8):1505–1516, 2004.
- Meyer and Matthies  Marcus Meyer and Hermann G Matthies. Efficient model reduction in non-linear dynamics using the karhunen-loéve expansion and dual-weighted-residual methods. Computational Mechanics, 31(1-2):179–191, 2003.
- Astrid  Patricia Astrid. Reduction of process simulation models: a proper orthogonal decomposition approach. Technische Universiteit Eindhoven Ph. D.-thesis, 2004.
- Jin and Durlofsky  Zhaoyang L Jin and Louis J Durlofsky. Reduced-order modeling of co 2 storage operations. International Journal of Greenhouse Gas Control, 68:49–67, 2018.
- Vermeulen et al.  PTM Vermeulen, Arnold W Heemink, and Chris BM Te Stroet. Reduced models for linear groundwater flow models using empirical orthogonal functions. Advances in water resources, 27(1):57–69, 2004.
- Vermeulen et al.  PTM Vermeulen, Chris BM Te Stroet, and Arnold W Heemink. Model inversion of transient nonlinear groundwater flow models using model reduction. Water resources research, 42(9):W09417, 2006.
- McPhee and Yeh  James McPhee and William WG Yeh. Groundwater management using model reduction via empirical orthogonal functions. Journal of Water Resources Planning and Management, 134(2):161–170, 2008.
- Siade et al.  Adam J Siade, Mario Putti, and William WG Yeh. Snapshot selection for groundwater model reduction using proper orthogonal decomposition. Water Resources Research, 46(8):W08539, 2010.
- Heijn et al.  T Heijn, Renato Markovinovic, and Jan D Jansen. Generation of low-order reservoir models using system-theoretical concepts. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 2003.
- Van Doren et al.  Jorn FM Van Doren, Renato Markovinović, and Jan D Jansen. Reduced-order optimal control of water flooding using proper orthogonal decomposition. Computational Geosciences, 10(1):137–158, 2006.
- Cardoso et al.  Marco A Cardoso, Louis J Durlofsky, and Pallav Sarma. Development and application of reduced-order modeling procedures for subsurface flow simulation. International journal for numerical methods in engineering, 77(9):1322–1350, 2009.
- Cardoso and Durlofsky  Marco A Cardoso and Louis J Durlofsky. Linearized reduced-order models for subsurface flow simulation. Journal of Computational Physics, 229(3):681–700, 2010.
- Pasetto et al.  Damiano Pasetto, Alberto Guadagnini, and Mario Putti. POD-based Monte-Carlo approach for the solution of regional scale groundwater flow driven by randomly distributed recharge. Advances in water resources, 34(11):1450–1463, 2011.
- Pasetto et al.  Damiano Pasetto, Mario Putti, and William WG Yeh. A reduced-order model for groundwater flow equation with random hydraulic conductivity: Application to Monte-Carlo methods. Water Resources Research, 49(6):3215–3228, 2013.
- Boyce and Yeh  Scott E Boyce and William WG Yeh. Parameter-independent model reduction of transient groundwater flow models: Application to inverse problems. Advances in water resources, 69:168–180, 2014.
- Baú  Domenico A Baú. Planning of groundwater supply systems subject to uncertainty using stochastic flow reduced models and multi-objective evolutionary optimization. Water resources management, 26(9):2513–2536, 2012.
- He et al.  Jincong He, Jonsat Sætrom, and Louis J Durlofsky. Enhanced linearized reduced-order models for subsurface flow simulation. Journal of Computational Physics, 230(23):8313–8341, 2011.
- Trehan and Durlofsky  Sumeet Trehan and Louis J Durlofsky. Trajectory piecewise quadratic reduced-order model for subsurface flow, with application to pde-constrained optimization. Journal of Computational Physics, 326:446–473, 2016.
- Rousset et al.  Matthieu AH Rousset, Chung K Huang, Hector Klie, and Louis J Durlofsky. Reduced-order modeling for thermal recovery processes. Computational Geosciences, 18(3-4):401–415, 2014.
- Jansen and Durlofsky  Jan Dirk Jansen and Louis J Durlofsky. Use of reduced-order models in well control optimization. Optimization and Engineering, 18(1):105–132, 2017.
- Chaturantabut and Sorensen  Saifon Chaturantabut and Danny C Sorensen. Nonlinear model reduction via discrete empirical interpolation. SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.
- Rewienski and White  Michal Rewienski and Jacob White. A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices. IEEE Transactions on computer-aided design of integrated circuits and systems, 22(2):155–170, 2003.
- He  Jincong He. Enhanced linearized reduced-order models for subsurface flow simulation. M.S. thesis, Stanford Univeristy, 2010.
- Bui-Thanh et al.  Tan Bui-Thanh, Karen Willcox, Omar Ghattas, and Bart VB Waanders. Goal-oriented, model-constrained optimization for reduction of large-scale systems. Journal of Computational Physics, 224(2):880–896, 2007.
- Wang et al.  Zhu Wang, Imran Akhtar, Jeff Borggaard, and Traian Iliescu. Proper orthogonal decomposition closure models for turbulent flows: a numerical comparison. Computer Methods in Applied Mechanics and Engineering, 237:10–26, 2012.
- Willcox  Karen Willcox. Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. Computers & fluids, 35(2):208–226, 2006.
- Barrault et al.  Maxime Barrault, Yvon Maday, Ngoc C Nguyen, and Anthony T. Patera. An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique, 339(9):667 – 672, 2004. ISSN 1631-073X. URL http://dx.doi.org/10.1016/j.crma.2004.08.006.
- Nguyen et al.  Ngoc C Nguyen, Anthony T. Patera, and Jaime Peraire. A ‘best points’ interpolation method for efficient approximation of parametrized functions. International journal for numerical methods in engineering, 73(4):521–543, 2008.
- Ghasemi  Mohammadreza Ghasemi. Model order reduction in porous media flow simulation and optimization. Ph.D. thesis, Texas AM Univeristy, 2015.
- Chaturantabut and Sorensen  Saifon Chaturantabut and Danny C Sorensen. A state space error estimate for pod-deim nonlinear model reduction. SIAM Journal on numerical analysis, 50(1):46–63, 2012.
- Xiao et al.  D Xiao, Fangxin Fang, Andrew G Buchan, Christopher C Pain, Ionel Michael Navon, Juan Du, and G Hu. Non-linear model reduction for the navier–stokes equations using residual deim method. Journal of Computational Physics, 263:1–18, 2014.
- Buffoni and Willcox  Marcelo Buffoni and Karen Willcox. Projection-based model reduction for reacting flows. In 40th Fluid Dynamics Conference and Exhibit, page 5008, 2010.
- Chaturantabut and Sorensen  Saifon Chaturantabut and Danny C Sorensen. Application of POD and DEIM on dimension reduction of non-linear miscible viscous fingering in porous media. Mathematical and Computer Modelling of Dynamical Systems, 17(4):337–353, 2011.
- Alghareeb and Williams  Zeid M Alghareeb and John Williams. Optimum decision-making in reservoir managment using reduced-order models. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2013.
- Yoon et al.  Seonkyoo Yoon, Zeid Alghareeb, and John Williams. Development of reduced-order oil reservoir models using localized deim. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2014.
- Carlberg et al.  Kevin Carlberg, Charbel Bou-Mosleh, and Charbel Farhat. Efficient non-linear model reduction via a least-squares Petrov–Galerkin projection and compressive tensor approximations. International Journal for Numerical Methods in Engineering, 86(2):155–181, 2011.
- Xiao et al.  Dunhui Xiao, Zhi Lin, Fangxin Fang, Christopher C Pain, Ionel M Navon, Pablo Salinas, and Ann Muggeridge. Non-intrusive reduced-order modeling for multiphase porous media flows using smolyak sparse grids. International Journal for Numerical Methods in Fluids, 83(2):205–219, 2017.
- Wang et al.  Zheng Wang, Dunhui Xiao, Fangxin Fang, Rajesh Govindan, Christopher C Pain, and Yike Guo. Model identification of reduced order fluid dynamics systems using deep learning. International Journal for Numerical Methods in Fluids, 86(4):255–268, 2017.
- Pascanu et al. [2013a] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. International Conference on Machine Learning (3), 28:1310–1318, 2013a.
- Mikolov et al.  Tomas Mikolov, Armand Joulin, Sumit Chopra, Michael Mathieu, and Marc’Aurelio Ranzato. Learning longer memory in recurrent neural networks. arXiv preprint arXiv:1412.7753, 2014.
- Pascanu et al. [2013b] Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, and Yoshua Bengio. How to construct deep recurrent neural networks. arXiv preprint arXiv:1312.6026, 2013b.
Irsoy and Cardie 
Ozan Irsoy and Claire Cardie.
Opinion mining with deep recurrent neural networks.
Empirical Methods in Natural Language Processing (EMNLP), pages 720–728, 2014.
- Hermans and Schrauwen  Michiel Hermans and Benjamin Schrauwen. Training and analysing deep recurrent neural networks. In Advances in Neural Information Processing Systems, pages 190–198, 2013.
- He et al.  Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. arXiv preprint arXiv:1512.03385, 2015.
- Hinton et al.  Geoffrey Hinton, Li Deng, Dong Yu, George E Dahl, Abdel-rahman Mohamed, Navdeep Jaitly, Andrew Senior, Vincent Vanhoucke, Patrick Nguyen, and Tara N Sainath. Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups. IEEE Signal Processing Magazine, 29(6):82–97, 2012.
- Graves  Alex Graves. Generating sequences with recurrent neural networks. arXiv preprint arXiv:1308.0850, 2013.
- Zimmermann et al.  Hans.G Zimmermann, Christoph Tietz, and Ralph Grothmann. Forecasting with recurrent neural networks: 12 tricks. In Neural Networks: Tricks of the Trade, pages 687–707. Springer, 2012.
- Bailer-Jones et al.  Coryn AL Bailer-Jones, David JC MacKay, and Philip J Withers. A recurrent neural network for modelling dynamical systems. network: computation in neural systems, 9(4):531–547, 1998.
- Bertsekas  Dimitri P Bertsekas. Nonlinear programming. Athena scientific Belmont, 1999.
- Tieleman and Hinton  Tijmen Tieleman and Geoffrey Hinton. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural Networks for Machine Learning, 4(2):26–31, 2012.
- Pletcher et al.  Richard H Pletcher, John C Tannehill, and Dale Anderson. Computational fluid mechanics and heat transfer. CRC Press, 2012.
- Chen et al.  Zhangxin Chen, Guanren Huan, and Yuanle Ma. Computational methods for multiphase flows in porous media. SIAM, 2006.
- Bastian  Peter Bastian. Numerical computation of multiphase flow in porous media. PhD thesis, habilitationsschrift Univeristät Kiel, 1999.
Matrix methods in data mining and pattern recognition, volume 4. SIAM, 2007.
- Trefethen and Bau III  Lloyd N Trefethen and David Bau III. Numerical linear algebra, volume 50. SIAM, 1997.
- Lucia et al.  David J Lucia, Philip S Beran, and Walter A Silva. Reduced-order modeling: new approaches for computational physics. Progress in Aerospace Sciences, 40(1):51–117, 2004.
- Werbos  Paul J Werbos. Backpropagation through time: what it does and how to do it. Proceedings of the IEEE, 78(Issue 10):1550–1560, 1990.
- Rumelhart et al.  David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning representations by back-propagating errors. nature, 323(6088):533, 1986.
- Paszke et al.  Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. In NIPS-W, 2017.