1 Introduction
Simulation of multiphase 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)
[4]. 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].Moreover, these subsurface properties are only known at a sparse set of points (i.e. well locations) and the grid properties are populated stochastically over the entire domain [6, 5, 3]
. MonteCarlo methods are usually employed to propagate the uncertainties in the subsurface properties to the flow response. MonteCarlo 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, Markovchain MonteCarlo (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 multiquery 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 fullorder model (full model) [7]
. A number of surrogate modeling techniques have been developed and could be broadly classified into three classes: simplified physics based models
[13, 11], datafit blackbox 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 highfidelity models using approaches such as simplifying physics assumptions, using coarse grids, and/or upscaling of the model parameters [13, 7, 9, 20]. Datafit 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. [23], Frangos et al. [7], Koziel and Leifsson [8], Razavi et al. [24].In projection based reduced order models (utilized in this paper), the governing equations of the full model are projected into a lowdimensional 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 inputoutput 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 [25], and Krylov subspace methods [26] use the characteristics of the full model and have been developed mainly for linear timeinvariant problems, although much progress has been done on extensions of these methods to nonlinear problems [27]. Snapshot based methods such as reduced basis methods [28], 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 MonteCarlo 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 [30] 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. [37] introduced POD in the confined, groundwater flow problems (linear subsurface flow model). Vermeulen et al. [38] applied POD in gradientbased optimization problem involving groundwater flow model. McPhee and Yeh [39] employed POD to enhance the groundwater management optimization problem. Siade et al. [40]
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. [41] and Van Doren et al. [42] applied POD to accelerate the optimization of a waterflood process. Cardoso et al. [43] incorporated a new snapshot clustering procedure to enhance the standard POD for oil–water subsurface flow problems.In the context of MonteCarlo 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. [45] 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. [46] applied POD to accelerate the MC simulations of transient confined groundwater flow models with stochastic hydraulic conductivity. Baú [48] derived a set of POD ROMs for each MC realization of hydraulic conductivity to solve a stochastic, multiobjective, confined groundwater management problem. Boyce and Yeh [47] applied a single parameterindependent 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 [52] has done an extensive review on the use of reducedorder 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 [50], 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 PODGalerkin 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 PODGalerkin 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 nonpolynomial 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) [54], gappy POD technique [58], missing point estimation (MPE) [59], best point interpolation method [60], and discrete empirical interpolation method (DEIM) [59, 53]. Among these techniques, TPWL and DEIM are widely used for efficient treatment of nonlinearities in multiphase flow reservoir simulations [61, 55, 9].
In TPWL method [54], the nonlinear function is first approximated by a piecewise linear function obtained by linearizing the fullorder 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 [54]. Finally, the reduced model can be obtained by projecting the resultant linearized fullorder system using standard techniques like POD [54]. The TPWL method was first introduced in [54] 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
[53]. The coefficients of expansion are determined by evaluating the nonlinear term only at a small number of selected interpolation points [53]. DEIM was developed in [53]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 [65] applied DEIM for model reduction of viscous fingering problems of an incompressible fluid through a two dimensional homogeneous porous medium. Alghareeb and Williams [66] combined DEIM with POD procedures and the resultant reduced model was applied in waterflood optimization problem. Recently, Ghasemi [61] applied POD with DEIM to an optimal control problem governed by twophase flow in a porous media. Next, Ghasemi [61]used machine learning technique to construct a number of PODDEIM local reducedorder models. In that work, machine learning technique was used to construct a number of PODDEIM local reducedorder models and then a specific local reducedorder model was selected with respect to the current state of the dynamical system during the gradient based optimization task. Similarly,
Yoon et al. [67] used multiple local DEIM approximations in POD reduced model framework to reduce the computational costs of highfidelity reservoir simulations.The overall convergence and stability is another issue that limits the applicability of PODGalerkin based ROMs. PODGalerkin 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 PODGalerkin ROMs is directly influenced by the number of POD basis used in the PODGalerkin projection. However, in many applications involving nonlinear conservation laws (e.g. high Reynold number fluid flow), PODGalerkin 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. [57] showed that PODGalerkin 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. [57] 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 PODGalerkin reduced model is to compute a new set of optimal basis or to improve the POD basis vectors by solving a constrained optimization problem. BuiThanh et al. [56] 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 PODGalerkin reduced model orthogonally projects the nonlinear residual into the subspace spanned by the POD basis vectors. Unlike PODGalerkin reduced model, PetrovGalerkin projection scheme design a different set of orthonormal basis called left reduced order basis into which the nonlinear residual is projected. Carlberg et al. [68] formulated stable PetrovGalerkin 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 [55] observed that poor spectral properties of the reduced Jacobian matrix could cause numerical instabilities in PODGalerkin TPWL reduced model. Hence, He [55] 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, datafit blackbox models have been combined with POD [69] to develop nonintrusive POD based ROMs, where the datafit models are used to regress the relationship between the input parameter and the reduced representation of the full model state vector. Hence, nonintrusive ROMs do not require any knowledge of the fullorder model and is mainly developed to circumvent the shortcomings in accessing the governing equations of the full model [69]. However, it can also be used to address the stability and nonlinearity issues of POD based ROMs. Wang et al. [70] developed a nonintrusive POD reduced model using Recurrent Neural Network (RNN) as a datafit 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 datafit 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 nonintrusive ROMs is severely undermined in many realworld problems, where increasing the dimensionality of the input parameter space increases the complexity and training time of the datafit model.
In summary, among many surrogate modeling techniques, PODGalerkin reduced model is a viable option for accelerating multiquery tasks like UQ. Generally, PODGalerkin 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 DRRNN [1] to alleviate the potential limitations of PODGalerkin reduced models. More specifically, we combine DRRNN with PODGalerkin and DEIM methods to derive an accurate and computationally effective reduced model for uncertainty quantification (UQ) tasks. The architecture of DRRNN is inspired by the iterative line search methods where the parameters of the DRRNN are optimized such that the residual of the numerically discretized PDEs is minimized [81, 82, 1]. Unlike the standard RNN which is very generic, DRRNN [1] uses the residual of the discretized differential equation. In addition, the parameters of the DRRNN are fitted such that the computed DRRNN output optimally minimizes the residual of the targeted equation. In this context, DRRNN 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 DRRNN combined with PODGalerkin and DEIM algorithm has a number of salient features. First, the dynamics of DRRNN 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 DRRNN are explicit in time, there is a reduction in the computational complexity of the extracted reduced model from corresponding to implicit PODDEIM reduced order models, to , where is the size of the reduced model. Third, DRRNN requires only very few training samples (obtained by solving the full model) to optimize the parameters of the DRRNN 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 datadriven algorithms (e.g. standard RNN architectures). Moreover, DRRNN can effectively emulate the parameterized nonlinear dynamical system with a significantly lower number of parameters in comparison to standard RNN architectures [1].
In this work, we demonstrate the superior properties of DRRNN in accelerating UQ tasks for subsurface reservoir models using MonteCarlo method. As far as we are aware, the use of a single parameterindependent PODGalerkin reduced model in MonteCarlo 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 DRRNN could be used to reconstruct stable solutions emulating the full model dynamics using only a small set of POD basis functions. The proposed DRRNN technique is validated on two forward uncertainty quantification problems involving twophase 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 [4]
. In these two numerical examples, the permeability field is modeled as lognormal distribution. The obtained results demonstrate that DRRNN combined with PODDEIM 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 DRRNN to formulate an accurate and computationally effective PODDEIM reduced model for flow problems where the standard PODGalerkin reduced models are inaccurate and possibly unstable. We also note that DRRNN architecture is generic and could be used to emulate any wellposed nonlinear dynamical system
[1] 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 multiphase flow problem in a porous media. In section 3, we introduce PODGalerkin method for model reduction followed by a discussion of DEIM for handling nonlinear systems. In Section 4, we describe the architecture of DRRNN and in section 5, we evaluate the reduced model derived by combining DRRNN with PODDEIM on two uncertainty quantification test cases. Finally, in Section 6, we present the conclusions of this manuscript.
2 Problem Formulation
The equations governing twophase flow of a wetting phase (water) and nonwetting 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
(1) 
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
(2) 
where is the phase velocity. The phase relative permeabilities and the capillary pressure () are usually modeled as functions of the phase saturations [4]. 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
(3) 
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
(4) 
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
(5) 
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
(6) 
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 multiphase 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 multiquery 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 PODGalerkin reduced model (POD reduced model) and PODDEIM reduced model where PODGalerkin 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 lowdimensional 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 [54].
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(7) 
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 fullorder model by a linear combination of basis vectors as following
(8) 
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
(9) 
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
(10) 
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
(11) 
The associated error termed as least–squares errors in approximating by using only basis vectors is given by
(12) 
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
(13) 
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.
3.3 PODGalerkin
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
(14) 
where and . Similarly, POD based reduced model for the saturation equation (Eq. (6)) takes the form
(15) 
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 [53] 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 .
3.4 Deim
Discrete Empirical Interpolation Method (DEIM) was introduced in [53] 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
(16) 
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
(17) 
The matrix factor in Eq. (17) is precomputed before solving Eq. (15). The overdetermined system is approximated using the DEIM algorithm introduced in [53] by first computing a matrix that selects rows of the matrix to obtain as following
(18) 
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
(19) 
where the matrix termed as the DEIM matrix. Similarly, the Jacobian of the nonlinear term in Eq. (15) is approximated using DEIM as following
(20) 
where is the Jacobian matrix computed using the components of evaluated by the DEIM algorithm [53, 54, 1]. Finally, the PODDEIM based reduced model takes the form
(21) 
We note that PODDEIM formulation is independent of the full model dimension and that the DEIM procedure exploits the structure of the nonlinear function as componentwise operation at [53].
4 Deep Residual RNN
PODDEIM reduced order models, as introduced in the last chapter, could be used to perform parametric UQ tasks. However, the PODDEIM 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 PODGalerkin reduced model in several cases and it was attributed to illconditioning and poor spectral properties of the reduced Jacobian matrix.
In this paper, we build on the recently introduced DRRNN [1] and formulate an accurate PODDEIM reduced order models. DRRNN is a deep RNN architecture [1], constructed by stacking physics aware network layers. DRRNN could be applied to any nonlinear dynamical system of the form
(22) 
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 [1]. The state variable at different time steps is obtained by solving the nonlinear residual equation defined as
(23) 
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. DRRNN [1] approximates the solution of Eq. (22) using the following iterative update equations
(24)  
where are the training parameters of DRRNN, is the activation function, is an elementwise 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
(25) 
where are fraction factors and is a smoothing term to avoid divisions by zero [1]. In this formulation, we set
. The architecture of DRRNN is inspired by the rmsprop algorithm
[82] which is a variant of the steepest descent method. The DRRNN output at each time step is defined as(26) 
The formulation of DRRNN is explicit in time and has a fixed number of iterations per time step. However, the dimension of the DRRNN system depends on the dimension of the residual. For example, DRRNN (Eq. (24)) can be derived from the PODDEIM reduced model residual (). In such setting, the DRRNN dynamics has a fixed computational budget of for each time step. In addition, DRRNN has the prospect of employing large time step violating the numerical stability constraint [1]. Furthermore, DRRNN does not rely on the reduced Jacobian matrix to approximate the solution of PODDEIM reduced model.
The DRRNN parameters are fitted by minimizing the mean square error (mse) defined by
(27) 
where (mse) is the average distance between the reference solution and the RNN output across a number of samples with timedependent observations [1, 73]. The set of parameters
is commonly estimated by a technique called Backpropagation Through Time (BPTT)
[89, 90, 71, 72], which backpropagates the gradient of the loss function
with 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 DRRNN against the standard implementation of PODGalerkin reduced model. Specifically, we develop two PODGalerkin based reduced model using DRRNN architecture namely, DRRNN (DRRNN combined with PODGalerkin) and DRRNN (DRRNN combined with PODGalerkin and DEIM). The numerical evaluations are performed using two uncertainty quantification tasks involving subsurface flow models. We did not include standard PODDEIM 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 finitevolume approach employed for obtaining the fullorder 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 DRRNN 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 Fullorder model setup
We consider a twophase (oil and water) porous media flow problem over the twodimensional domain meters. The equations governing the twophase 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 [4]. 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 lognormal distribution function with zero mean and exponential covariance kernel of the form
(28) 
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 logpermeability values.
For solving the fullorder 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) [4] and an upwind finitevolume scheme is used to discretized the saturation equation. For the time discretization, an implicit backward Euler method combined with NewtonRaphson 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 nondimensional 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 MonteCarlo method [6].
5.2 PODGalerkin 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 [61]. 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 DRRNN setup
In all the numerical test cases, we utilize DRRNN with six layers ( in Eq. (24)). We evaluate DRRNN and DRRNN for different number of POD basis, however, we fix the number of DEIM basis to . The PyTorch framework [91]
, a deep learning python package using
Torch library as a backend, is used to implement the DRRNN. Further, we optimize the DRRNN 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 hyperparameters
and in Eq. (25) to and , respectively. The formulated DRRNN and DRRNN 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 DRRNN by minimizing the loss function defined in Eq. (27).5.4 Evaluation metrics
We evaluate the performance of DRRNN and DRRNN using two time specific error metrics defined by
(29)  
where is the realization index, and is computed from the reduced model. Additionally, we utilize two relative error metrics defined as
(30)  
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 subplots 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 DRDRRNN and DRRNN results are very close to the least–square solutions (LS fit). In Figure 3, PODGalerkin reduced model yields extremely inaccurate and unstable results. We attribute the small errors in DRRNN and DRRNN results to the insufficient number of POD basis vectors and we note that the error magnitude is equivalent to the optimal values obtained by leastsquares projection.
show the results for the first (mean) and second (standard deviation) moments of the saturation field at time
PVI 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 DRRNN 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 DRRNN 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 DRDRRNN and DRRNN are indistinguishable from the plots obtained from the LS fit (the best approximation). Further, we observe that the saturation PDF obtained from DRDRRNN and DRRNN 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 DRRNN and DRRNN are nearly indistinguishable from the least–squares projection errors.
Error  #Basis  Reduced Order Models  

LS fit  POD  DRRNN  DRRNN  
10  0.13  0.56  0.14  0.15  
20  0.10  2.7  0.11  0.13  
10  0.20  1.8  0.20  0.27  
20  0.17  5.8  0.19  0.26 
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 DRRNN and DRRNN 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 [55], where it was shown that selecting large number of basis vectors based on singular values may not lead to stable PODGalerkin reduced model. Further, it was presented in [55] that the relation between the stability property of PODGalerkin reduced model and the number of basis vectors used in PODGalerkin 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 DRRNN, DRRNN, 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 DRRNN 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 DRRNN, DRRNN 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 unimodal distribution. Similar to test case 1, POD reduced model yields inaccurate approximation of the saturation PDFs.
Error  #Basis  Reduced Order Models  

LS fit  POD  DRRNN  DRRNN  
10  0.09  1.30  0.10  0.12  
20  0.07  2.05  0.08  0.10  
10  0.19  3.5  0.21  0.22  
20  0.16  6.2  0.18  0.22 
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 DRRNN and DRRNN are almost close to the least–squares (best approximation) approximation errors. However, the POD reduced model yields very inaccurate results due to numerical instabilities.
6 Conclusion
In this work, we extended the DRRNN introduced in [1] into nonlinear multiphase flow problem with distributed uncertain parameters. In this extended formulation, DRRNN based on the reduced residual obtained from PODDEIM reduced model is used to construct the reduced order model termed DRRNN. We evaluated the proposed DRRNN on two forward uncertainty quantification problems involving twophase flow in subsurface porous media. The uncertainty parameter is the permeability field modeled as lognormal distribution. In the two test cases, full order model and ROMs are solved for random permeability realizations to estimate an ensemble based statistics using MonteCarlo method. full model and POD reduced model used implicit time stepping method as the time step size violates the numerical stability condition. However, DRRNN architecture employs explicit time stepping procedure for the same step size used in full model and POD reduced model. Hence, DRRNN 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 DRRNN and is the average number of Newton iterations used in the standard PODDEIM reduced model. The obtained numerical results shows that DRRNN 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 DRRNN 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 DRRNN for UQ tasks with the permeability fields that has randomly oriented channels or barriers. The use of DRRNN 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.
References
 Nagoor Kani and Elsheikh [2017] J. Nagoor Kani and A. H. Elsheikh. DRRNN: A deep residual recurrent neural network for model reduction. ArXiv eprints, September 2017.
 Petvipusit et al. [2014] 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 15731499. URL https://doi.org/10.1007/s105960149425z.
 Elsheikh et al. [2013] 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 19447973. URL http://dx.doi.org/10.1002/2012WR013406.
 Aarnes et al. [2007] 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. [2012] 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 18748953. URL https://doi.org/10.1007/s1100401293972.
 Ibrahima [2016] Fayadhoi Ibrahima. Probability distribution methods for nonlinear transport in heterogenous porous media. Ph.D. thesis, Stanford Univeristy, 2016.
 Frangos et al. [2010] Michalis Frangos, Youssef Marzouk, Karen Willcox, and Bart VB Waanders. Surrogate and ReducedOrder Modeling: A Comparison of Approaches for LargeScale 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 [2013] Slawomir Koziel and Leifur Leifsson. Surrogatebased modeling and optimization. Applications in Engineering, 2013.
 He [2013] Jincong He. Reducedorder modeling for oilwater and compositional systems, with application to data assimilation and production optimization. Ph.D. thesis, Stanford Univeristy, 2013.
 Elsheikh et al. [2014] 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. [2015] 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. [2015] 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 [2012] Louis J Durlofsky and Yuguang Chen. Uncertainty quantification for subsurface flow problems using coarsescale models. In Numerical analysis of multiscale problems, pages 163–202. Springer, 2012.
 Li et al. [2017] Hao Li, Zhien Zhang, and Zhijian Liu. Application of artificial neural networks for catalysis: a review. Catalysts, 7(10):306, 2017.
 Yeten et al. [2005] 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. [1993] 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. [2014] 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. [2001] Athanasios C Antoulas, Danny C Sorensen, and Serkan Gugercin. A survey of model reduction methods for largescale systems. Contemporary mathematics, 280:193–220, 2001.
 Fang et al. [2013] F. Fang, C.C. Pain, I.M. Navon, A.H. Elsheikh, J. Du, and D. Xiao. Nonlinear 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. [2013] 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/s1124201301497.
 AbdiKhanghah et al. [2018] Mahdi AbdiKhanghah, Amin Bemani, Zahra Naserzadeh, and Zhien Zhang. Prediction of solubility of nalkanes in supercritical co 2 using rbfann and mlpann. Journal of CO2 Utilization, 25:108–119, 2018.
 A.Wood [2018] David A.Wood. transparent openbox learning network provides insight to complex systems and a performance benchmark for moreopaque machine learning algorithms. Advances in GeoEnergy Research, 2(2):148–162, 2018.
 Asher et al. [2015] 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. [2012] 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 [2004] 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 [2003] Roland W Freund. Model reduction methods based on Krylov subspaces. Acta Numerica, 12:267–319, 2003.
 Lall et al. [2002] 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. [2007] 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 [1987] Lawrence Sirovich. Turbulence and the dynamics of coherent structures. i. coherent structures. Quarterly of applied mathematics, 45(3):561–571, 1987.
 Lumley [1967] John L Lumley. The structure of inhomogeneous turbulent flows. Atmospheric turbulence and radio wave propagation, 1967.
 Zheng et al. [2002] Daguang Zheng, Karlene A Hoo, and Michael J Piovoso. Loworder 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. [2006] Yanhua Cao, Jiang Zhu, Zhendong Luo, and Ionel M Navon. Reducedorder modeling of the upper tropical pacific ocean model using proper orthogonal decomposition. Computers & Mathematics with Applications, 52(89):1373–1386, 2006.
 BuiThanh et al. [2004] Tan BuiThanh, 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 [2003] Marcus Meyer and Hermann G Matthies. Efficient model reduction in nonlinear dynamics using the karhunenloéve expansion and dualweightedresidual methods. Computational Mechanics, 31(12):179–191, 2003.
 Astrid [2004] Patricia Astrid. Reduction of process simulation models: a proper orthogonal decomposition approach. Technische Universiteit Eindhoven Ph. D.thesis, 2004.
 Jin and Durlofsky [2018] Zhaoyang L Jin and Louis J Durlofsky. Reducedorder modeling of co 2 storage operations. International Journal of Greenhouse Gas Control, 68:49–67, 2018.
 Vermeulen et al. [2004] 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. [2006] 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 [2008] 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. [2010] 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. [2003] T Heijn, Renato Markovinovic, and Jan D Jansen. Generation of loworder reservoir models using systemtheoretical concepts. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 2003.
 Van Doren et al. [2006] Jorn FM Van Doren, Renato Markovinović, and Jan D Jansen. Reducedorder optimal control of water flooding using proper orthogonal decomposition. Computational Geosciences, 10(1):137–158, 2006.
 Cardoso et al. [2009] Marco A Cardoso, Louis J Durlofsky, and Pallav Sarma. Development and application of reducedorder modeling procedures for subsurface flow simulation. International journal for numerical methods in engineering, 77(9):1322–1350, 2009.
 Cardoso and Durlofsky [2010] Marco A Cardoso and Louis J Durlofsky. Linearized reducedorder models for subsurface flow simulation. Journal of Computational Physics, 229(3):681–700, 2010.
 Pasetto et al. [2011] Damiano Pasetto, Alberto Guadagnini, and Mario Putti. PODbased MonteCarlo 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. [2013] Damiano Pasetto, Mario Putti, and William WG Yeh. A reducedorder model for groundwater flow equation with random hydraulic conductivity: Application to MonteCarlo methods. Water Resources Research, 49(6):3215–3228, 2013.
 Boyce and Yeh [2014] Scott E Boyce and William WG Yeh. Parameterindependent model reduction of transient groundwater flow models: Application to inverse problems. Advances in water resources, 69:168–180, 2014.
 Baú [2012] Domenico A Baú. Planning of groundwater supply systems subject to uncertainty using stochastic flow reduced models and multiobjective evolutionary optimization. Water resources management, 26(9):2513–2536, 2012.
 He et al. [2011] Jincong He, Jonsat Sætrom, and Louis J Durlofsky. Enhanced linearized reducedorder models for subsurface flow simulation. Journal of Computational Physics, 230(23):8313–8341, 2011.
 Trehan and Durlofsky [2016] Sumeet Trehan and Louis J Durlofsky. Trajectory piecewise quadratic reducedorder model for subsurface flow, with application to pdeconstrained optimization. Journal of Computational Physics, 326:446–473, 2016.
 Rousset et al. [2014] Matthieu AH Rousset, Chung K Huang, Hector Klie, and Louis J Durlofsky. Reducedorder modeling for thermal recovery processes. Computational Geosciences, 18(34):401–415, 2014.
 Jansen and Durlofsky [2017] Jan Dirk Jansen and Louis J Durlofsky. Use of reducedorder models in well control optimization. Optimization and Engineering, 18(1):105–132, 2017.
 Chaturantabut and Sorensen [2010] 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 [2003] Michal Rewienski and Jacob White. A trajectory piecewiselinear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices. IEEE Transactions on computeraided design of integrated circuits and systems, 22(2):155–170, 2003.
 He [2010] Jincong He. Enhanced linearized reducedorder models for subsurface flow simulation. M.S. thesis, Stanford Univeristy, 2010.
 BuiThanh et al. [2007] Tan BuiThanh, Karen Willcox, Omar Ghattas, and Bart VB Waanders. Goaloriented, modelconstrained optimization for reduction of largescale systems. Journal of Computational Physics, 224(2):880–896, 2007.
 Wang et al. [2012] 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 [2006] Karen Willcox. Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. Computers & fluids, 35(2):208–226, 2006.
 Barrault et al. [2004] Maxime Barrault, Yvon Maday, Ngoc C Nguyen, and Anthony T. Patera. An empirical interpolation method: application to efficient reducedbasis discretization of partial differential equations. Comptes Rendus Mathematique, 339(9):667 – 672, 2004. ISSN 1631073X. URL http://dx.doi.org/10.1016/j.crma.2004.08.006.
 Nguyen et al. [2008] 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 [2015] Mohammadreza Ghasemi. Model order reduction in porous media flow simulation and optimization. Ph.D. thesis, Texas AM Univeristy, 2015.
 Chaturantabut and Sorensen [2012] Saifon Chaturantabut and Danny C Sorensen. A state space error estimate for poddeim nonlinear model reduction. SIAM Journal on numerical analysis, 50(1):46–63, 2012.
 Xiao et al. [2014] D Xiao, Fangxin Fang, Andrew G Buchan, Christopher C Pain, Ionel Michael Navon, Juan Du, and G Hu. Nonlinear model reduction for the navier–stokes equations using residual deim method. Journal of Computational Physics, 263:1–18, 2014.
 Buffoni and Willcox [2010] Marcelo Buffoni and Karen Willcox. Projectionbased model reduction for reacting flows. In 40th Fluid Dynamics Conference and Exhibit, page 5008, 2010.
 Chaturantabut and Sorensen [2011] Saifon Chaturantabut and Danny C Sorensen. Application of POD and DEIM on dimension reduction of nonlinear miscible viscous fingering in porous media. Mathematical and Computer Modelling of Dynamical Systems, 17(4):337–353, 2011.
 Alghareeb and Williams [2013] Zeid M Alghareeb and John Williams. Optimum decisionmaking in reservoir managment using reducedorder models. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2013.
 Yoon et al. [2014] Seonkyoo Yoon, Zeid Alghareeb, and John Williams. Development of reducedorder oil reservoir models using localized deim. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2014.
 Carlberg et al. [2011] Kevin Carlberg, Charbel BouMosleh, and Charbel Farhat. Efficient nonlinear model reduction via a leastsquares Petrov–Galerkin projection and compressive tensor approximations. International Journal for Numerical Methods in Engineering, 86(2):155–181, 2011.
 Xiao et al. [2017] Dunhui Xiao, Zhi Lin, Fangxin Fang, Christopher C Pain, Ionel M Navon, Pablo Salinas, and Ann Muggeridge. Nonintrusive reducedorder 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. [2017] 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. [2014] 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 [2014]
Ozan Irsoy and Claire Cardie.
Opinion mining with deep recurrent neural networks.
In
Empirical Methods in Natural Language Processing (EMNLP)
, pages 720–728, 2014.  Hermans and Schrauwen [2013] 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. [2015] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. arXiv preprint arXiv:1512.03385, 2015.
 Hinton et al. [2012] Geoffrey Hinton, Li Deng, Dong Yu, George E Dahl, Abdelrahman 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 [2013] Alex Graves. Generating sequences with recurrent neural networks. arXiv preprint arXiv:1308.0850, 2013.
 Zimmermann et al. [2012] 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.
 BailerJones et al. [1998] Coryn AL BailerJones, 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 [1999] Dimitri P Bertsekas. Nonlinear programming. Athena scientific Belmont, 1999.
 Tieleman and Hinton [2012] Tijmen Tieleman and Geoffrey Hinton. Lecture 6.5rmsprop: 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. [2012] Richard H Pletcher, John C Tannehill, and Dale Anderson. Computational fluid mechanics and heat transfer. CRC Press, 2012.
 Chen et al. [2006] Zhangxin Chen, Guanren Huan, and Yuanle Ma. Computational methods for multiphase flows in porous media. SIAM, 2006.
 Bastian [1999] Peter Bastian. Numerical computation of multiphase flow in porous media. PhD thesis, habilitationsschrift Univeristät Kiel, 1999.

Eldén [2007]
Lars Eldén.
Matrix methods in data mining and pattern recognition
, volume 4. SIAM, 2007.  Trefethen and Bau III [1997] Lloyd N Trefethen and David Bau III. Numerical linear algebra, volume 50. SIAM, 1997.
 Lucia et al. [2004] David J Lucia, Philip S Beran, and Walter A Silva. Reducedorder modeling: new approaches for computational physics. Progress in Aerospace Sciences, 40(1):51–117, 2004.
 Werbos [1990] 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. [1986] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning representations by backpropagating errors. nature, 323(6088):533, 1986.
 Paszke et al. [2017] 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 NIPSW, 2017.
Comments
There are no comments yet.