Reduced-Order Modeling of Subsurface Multi-phase Flow Models Using Deep Residual Recurrent Neural Networks

We present a reduced-order modeling technique for subsurface multi-phase ﬂow problems building on the recently introduced deep residual recurrent neural network (DR-RNN) (Nagoor Kani et al. in DR-RNN: a deep residual recurrent neural network for model reduction. ArXiv e-prints, 2017). DR-RNN is a physics-aware recurrent neural network for modeling the evolution of dynamical systems. The DR-RNN architecture is inspired by iterative update techniques of line search methods where a ﬁxed number of layers are stacked together to minimize the residual (or reduced residual) of the physical model under consideration. In this manuscript, we combine DR-RNN with proper orthogonal decomposition (POD) and discrete empirical interpolation method (DEIM) to reduce the computational complexity associated with high-ﬁdelity numerical simulations. In the presented formulation, POD is used to construct an optimal set of reduced basis functions and DEIM is employed to evaluate the nonlinear terms independent of the full-order model size. We demonstrate the proposed reduced model on two uncertainty quantiﬁcation test cases using Monte Carlo simulation of subsurface ﬂow with random permeability ﬁeld. The obtained results demonstrate that DR-RNN combined with POD–DEIM provides an accurate and stable reduced model with a ﬁxed computational budget that is much less than the computational cost of standard POD–Galerkin reduced model combined with DEIM for nonlinear dynamical systems.


Introduction
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 (Petvipusit et al. 2014;Elsheikh et al. 2013 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) (Aarnes et al. 2007). 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 exhibits multi-scale features (Elsheikh et al. 2012;Petvipusit et al. 2014).
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 (Ibrahima 2016;Elsheikh et al. 2012Elsheikh et al. , 2013. 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 (Petvipusit et al. 2014;Elsheikh et al. 2013;Ibrahima 2016). 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 (Elsheikh et al. 2012. 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 (Frangos et al. 2010;Koziel and Leifsson 2013;He 2013;Elsheikh et al. 2014;Josset et al. 2015;Bazargan et al. 2015). Surrogate models are computationally efficient mathematical models that can effectively approximate the main characteristics of the full-order model (full model) (Frangos et al. 2010). A number of surrogate modeling techniques have been developed and could be broadly classified into three classes: simplified physics-based models (Durlofsky and Chen 2012;Josset et al. 2015), data-fit black-box models (Frangos et al. 2010;Li et al. 2017;Yeten et al. 2005), and projection-based reduced-order models commonly referred to as reduced model (Berkooz et al. 1993;Lassila et al. 2014;Antoulas et al. 2001;Fang et al. 2013). 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 (Durlofsky and Chen 2012;Frangos et al. 2010;He 2013;Babaei et al. 2013). Data-fit models are generated using the detailed simulation data to regress the relation between the input and the corresponding output of interest (Frangos et al. 2010;Yeten et al. 2005;Abdi-Khanghah et al. 2018;Wood 2018). For a complete review of various surrogate modeling techniques, we refer the readers to the following papers by Asher et al. (2015), Frangos et al. (2010), Koziel and Leifsson (2013) and Razavi et al. (2012).
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 (Lassila et al. 2014;Antoulas et al. 2001). Projectionbased ROMs rely 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 (Frangos et al. 2010;Lassila et al. 2014;Antoulas et al. 2001). Projection-based reduced-order models are broadly categorized into system-based methods and snapshot-based methods. System-based methods like balanced truncation realization methods (Gugercin and Antoulas 2004) and Krylov subspace methods (Freund 2003) 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 (Lall et al. 2002). Snapshot-based methods such as reduced basis methods (Rozza et al. 2007) and proper orthogonal decomposition (POD) (Sirovich 1987;Berkooz et al. 1993) 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 are 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 (Lassila et al. 2014;Sirovich 1987) (see Sect. 3 for further details). Lumley (1967) was the first to apply POD techniques in fluid flow simulations. Since then, POD procedures have successfully been applied in a number of application areas (e.g., Sirovich 1987;Zheng et al. 2002;Cao et al. 2006;Bui-Thanh et al. 2004;Meyer and Matthies 2003;Astrid 2004;Jin and Durlofsky 2018).
In the context of fluid flow in porous media, Vermeulen et al. (2004) introduced POD in the confined, groundwater flow problems (linear subsurface flow model). Vermeulen et al. (2006) applied POD in gradient-based optimization problem involving groundwater flow model. McPhee and Yeh (2008) employed POD to enhance the groundwater management optimization problem. Siade et al. (2010) 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. (2003) and Van Doren et al. (2006) applied POD to accelerate the optimization of a waterflood process. Cardoso et al. (2009) 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) (Cardoso and Durlofsky 2010;Pasetto et al. 2011Pasetto et al. , 2013Boyce and Yeh 2014). Pasetto et al. (2011) employed POD-based reduced model to construct MC realizations of twodimensional steady-state confined groundwater flow subject to a spatially distributed random recharge. Pasetto et al. (2013) applied POD to accelerate the MC simulations of transient confined groundwater flow models with stochastic hydraulic conductivity. Baú (2012) 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 (2014) 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) (Cardoso and Durlofsky 2010;He et al. 2011;Trehan and Durlofsky 2016;Rousset et al. 2014;Jansen and Durlofsky 2017). Recently, Jansen and Durlofsky (2017) 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 Trehan and Durlofsky (2016), which impose some additional computational overhead. Although POD combined with Galerkin projection has been applied more frequently to nonlinear flow problems (Bui-Thanh et al. 2004;Berkooz et al. 1993;Rousset et al. 2014), 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 (Chaturantabut and Sorensen 2010;Rewienski and White 2003;Cardoso and Durlofsky 2010), and the second factor is related to maintaining the overall stability of the resulting reduced model (Cardoso and Durlofsky 2010;He 2010He , 2013Bui-Thanh et al. 2007;Wang et al. 2012).
In relation to computing reduced non-polynomial nonlinear functions, POD-based ROMs are 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) (Rewienski and White 2003), gappy POD technique (Willcox 2006), missing point estimation (MPE) (Barrault et al. 2004), best point interpolation method (Nguyen et al. 2008), and discrete empirical interpolation method (DEIM) (Barrault et al. 2004;Chaturantabut and Sorensen 2010). Among these techniques, TPWL and DEIM are widely used for efficient treatment of nonlinearities in multi-phase flow reservoir simulations (Ghasemi 2015;He 2010He , 2013. In TPWL method (Rewienski and White 2003), 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 (Rewienski and White 2003). Finally, the reduced model can be obtained by projecting the resultant linearized fullorder system using standard techniques like POD (Rewienski and White 2003). The TPWL method was first introduced in Rewienski and White (2003) for modeling nonlinear circuits and micromachined devices. In the context of subsurface flow problems, TPWL procedures were applied in Cardoso and Durlofsky (2010), He et al. (2011), Trehan and Durlofsky (2016) and Rousset et al. (2014) 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 (Chaturantabut and Sorensen 2010). The coefficients of expansion are determined by evaluating the nonlinear term only at a small number of selected interpolation points (Chaturantabut and Sorensen 2010). DEIM was developed in Chaturantabut and Sorensen (2010) for model reduction of general nonlinear system of ordinary differential equations (ODEs) and has been used in several areas (Chaturantabut and Sorensen 2012;Xiao et al. 2014;Buffoni and Willcox 2010). Within the context of subsurface flow problems, Chaturantabut and Sorensen (2011) applied DEIM for model reduction of viscous fingering problems of an incompressible fluid through a two-dimensional homogeneous porous medium. Alghareeb and Williams (2013) combined DEIM with POD procedures, and the resultant reduced model was applied in waterflood optimization problem. Recently, Ghasemi (2015) applied POD with DEIM to an optimal control problem governed by two-phase flow in a porous media. Next, Ghasemi (2015) 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. (2014) 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 (Wang et al. 2012;He 2010;Bui-Thanh et al. 2007) 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 Reynolds number fluid flow), POD-Galerkin reduced-order models have shown poor performance even after retaining a sufficient number of POD basis (Wang et al. 2012;Sirovich 1987;Berkooz et al. 1993).
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 (Berkooz et al. 1993;Wang et al. 2012). 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 (Berkooz et al. 1993;Wang et al. 2012). Wang et al. (2012) 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 Reynolds number of 1000. Wang et al. (2012) 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. (2007) 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 designs a different set of orthonormal basis called left reduced-order basis into which the nonlinear residual is projected. Carlberg et al. (2011) formulated stable Petrov-Galerkin reduced model in which the left reducedorder basis vectors were computed from an optimization problem at every iteration of the Gauss Newton method. He (2010) observed that poor spectral properties of the reduced Jacobian matrix could cause numerical instabilities in POD-Galerkin TPWL reduced model. Hence, He (2010) 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 (He 2010(He , 2013. 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 are 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. (2017) 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 (Pascanu et al. 2013a;Mikolov et al. 2014) which has at least one feedback connection in addition to the feedforward connections (Pascanu et al. 2013a, b;Irsoy and Cardie 2014). 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 (Hermans and Schrauwen 2013;He et al. 2015;Hinton et al. 2012;Graves 2013). Additionally, RNN has been applied to emulate the evolution of nonlinear dynamical systems in a number of applications (Zimmermann et al. 2012;Bailer-Jones et al. 1998) 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 (Cardoso and Durlofsky 2010;He et al. 2011;Trehan and Durlofsky 2016;Ghasemi 2015). 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 (Wang et al. 2012;He 2010He , 2013 are not cost-effective and ultimately do not guarantee stability of the extracted reduced-order models. In this paper, we use DR-RNN (Nagoor Kani and Elsheikh 2017) 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 (Bertsekas 1999;Tieleman and Hinton 2012;Nagoor Kani and Elsheikh 2017). Unlike the standard RNN which is very generic, DR-RNN (Nagoor Kani and Elsheikh 2017) 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 (Nagoor Kani and Elsheikh 2017; Pletcher et al. 2012). 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 O(r 3 ) corresponding to implicit POD-DEIM reduced-order models, to O(r 2 ), where r 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 with standard RNN architectures (Nagoor Kani and Elsheikh 2017).
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 (Cardoso and Durlofsky 2010;He et al. 2011;Boyce and Yeh 2014;Ghasemi 2015). 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 Sect. 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 (Aarnes et al. 2007). 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 (Nagoor Kani and Elsheikh 2017) 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 Sect. 2, we present the formulation of multi-phase flow problem in a porous media. In Sect. 3, we introduce POD-Galerkin method for model reduction followed by a discussion of DEIM for handling nonlinear systems. In Sect. 4, we describe the architecture of DR-RNN, and in Sect. 5, we evaluate the reduced model derived by combining DR-RNN with POD-DEIM on two uncertainty quantification test cases. Finally, in Sect. 6, we present the conclusions of this manuscript.

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 (Aarnes et al. 2007;He 2013;Chen et al. 2006;Bastian 1999). The continuity equation for each phase α takes the form where the subscript α = w denotes the water phase, the subscript α = o denotes the oil phase, K is the absolute permeability tensor, λ α = k r α /μ α is the phase mobility, with k r α the relative permeability to phase α and μ α the viscosity of phase α, p α is the phase pressure, ρ α is the density of phase α, g is the gravitational acceleration, h is the depth, φ is the porosity, s α is the saturation of the phase α and q α is the phase source and sink terms (Aarnes et al. 2007;Chen et al. 2006). Further, the phase saturations are constrained by s w + s o = 1, since the oil and the water jointly fill the void space (Aarnes et al. 2007;He 2013). The phase velocities are modeled by the multi-phase Darcy's law to relate the phase velocities to the phase pressures and take the form where v α is the phase velocity. The phase relative permeabilities k r α and the capillary pressure ( p cow = p o − p w ) are usually modeled as functions of the phase saturations (Aarnes et al. 2007). 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 (Aarnes et al. 2007;He 2013;Bastian 1999). The simplified global pressure equation takes the form where p = p o = p w is the global pressure, λ = λ w + λ o is the total mobility, q = q w + q o is the source and sink term. The saturation equation for the water phase takes the following form is a function of saturation termed as the fractional flow function for the water phase, v = −Kλ ∇ p is the total velocity vector, and s = s w is the water saturation (Aarnes et al. 2007;Chen et al. 2006). In the rest of the paper, we write the water phase saturation as s = s w for simplicity. Coupled equations Eqs. (3) and (4) could then be solved for the evolution of the saturation by providing the appropriate initial and boundary conditions. Equations (3) and (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 n grid blocks and then applying the finite-volume method to discretize the spatial derivatives of Eqs. (3) and (4). The discretized pressure equation takes the form where A ∈ R n×n , b ∈ R n , and y p ∈ R n is the pressure vector in which each component y p i of y p represent the pressure value at the ith grid block. Similarly, the spatially discretized saturation equation takes the form where B ∈ R n×n , d ∈ R n , v is the total velocity vector, and y s ∈ R n is the saturation vector in which each component y s i of y s is the saturation value at the ith grid block.
Equations (5) and (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 A on the mobilities λ(y s (t)) in the pressure full model [Eq. (5)] and from the dependence of the matrix B on the velocity vector v(y p ) in the saturation full model [Eq. (6)]. In this paper, we adopt an implicit sequential splitting method to solve the full model [Eqs. (5) and (6)]. In this method, the saturation vector y s (t) from the present time step is used to assemble the matrix A in Eq. (5) and then the pressure full model [Eq. (5)] is solved for the pressure vector y p . Following that, the velocity vector v (computed from the pressure vector y p ) is used to assemble the matrix B 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 Eqs. (5) and (6), when n (the number of grid block) is large (Chaturantabut and Sorensen 2010;Ghasemi 2015). nonlinear terms. Both methods are introduced to reduce the computational effort associated with solving the full model [Eqs. (5) and (6)].

POD Basis
As stated in Sect. 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 U spanned by a small set of r 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 (Rewienski and White 2003).
The optimal basis vectors {u i } r i=1 in POD are computed by singular value decomposition (SVD) of the solution snapshot matrix X. The solution snapshot matrix X is obtained from a set of solution vectors of size n s obtained by solving the full model at selected points in the input parameter space. The SVD of X is expressed as where X ∈ R n×n s , U = [u 1 u 2 u 3 · · · u n ] ∈ R n×n is the left singular matrix and Σ = diag(σ 1 > σ 2 > σ 3 > · · · σ ns ≥ 0) is the diagonal matrix containing the singular values σ i of the snapshot matrix X in descending order. The dominant left singular vectors corresponding to the first r largest singular values represent the basis vectors to span the optimal subspace U of POD-based reduced model. Thus, the first step in deriving the PODbased reduced model is to express the state vector y of the full-order model by a linear combination of r basis vectors as follows: whereỹ ∈ R r is the reduced state vector representation of full-dimensional state vector y, and U r = [u 1 · · · u r ] ∈ R n×r is the matrix that contains r orthonormal basis vectors in its columns. By following this step, for example, the optimal basis vectors for the saturation state vector y s are obtained from the SVD of the saturation snapshot matrix X s = (y s 1 . . . y s T ) 1 . . . (y s 1 . . . y s T ) L , where T is the number of time steps and L is the number of samples of input parameter used to build the snapshot matrix. The SVD of X s is expressed as where U s ∈ R n×n is the left singular matrix, and Σ s is the diagonal matrix containing the singular values of the snapshot matrix X s in descending order. The saturation state vector y s is optimally expressed as y s ≈ U r sỹ s whereỹ s ∈ R r is the reduced state vector representation of y s , and U r s ∈ R n×r is the matrix that contains r orthonormal basis vectors in its columns. Similarly, we can represent the pressure state vector y p from its reduced state vector representationỹ p using optimal basis matrix U p obtained from the SVD of the pressure snapshot matrix X p .

Least-Squares Approximation
The capacity of a set of basis functions to represent a new solution vector could be tested using least-squares fitting (Eldén 2007;Trefethen and Bau III 1997). For example, the least-squares solution for approximating a saturation state vector y * s ∈ R n is defined as The associated error termed as least-squares errors in approximating y s by y * s using only r basis vectors is given by The least-squares error ε s [Eq. (12)] is equivalent to the omitted energy Ω s = n i=r +1 σ s i (Lucia et al. 2004;Berkooz et al. 1993). In practice, r is commonly chosen as the smallest integer such that the relative omitted energy ν is less than a preset value (e.g., 0.01), where the omitted energy is defined by the following equation Similar expressions mentioned in Eqs. (11), (12), and (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 (Frangos et al. 2010;Lucia et al. 2004). The least-squares solution is the best approximation of the state variables in the sense that, for the chosen lowdimensional subspace U, no other low-dimensional approximation can represent the given snapshot set with a lower least-squares error (Lassila et al. 2014;Sirovich 1987;Berkooz et al. 1993). 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 Sect. 5.

POD-Galerkin
Once the POD basis vectors are obtained, the reduced representation of the pressure vector y p is substituted into the pressure full model [Eq. (5)], followed by Galerkin projection of the pressure equation into the subspace spanned by U r p . The resulting POD-based reduced model for the pressure equation then takes the following form whereÃ = U r p A U r p ∈ R r ×r andb = U r p b ∈ R r . Similarly, POD-based reduced model for the saturation equation [Eq. (6)] takes the form whered = U r s d andd ∈ R r . The POD-based reduced model formulated by Eqs. (14) and (15) is of the reduced dimension r . However, the nonlinear function f w in Eq. (15) is still of the order of full dimension n. Moreover, the reduced Jacobian matrixJ =Ĩ − U r s B J f (f w (U r sỹ s ))U r s ∈ R r ×r needed for Newton-like iterations to solve this nonlinear equation is also of order n (Chaturantabut and Sorensen 2010) as it relies on evaluating the full-order nonlinear function f w . 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 n.

DEIM
Discrete empirical interpolation method (DEIM) was introduced in Chaturantabut and Sorensen (2010) to approximate the nonlinear terms in POD-based reduced model using a limited number of points that are independent of the full system dimension n. Similar to POD, the first step of DEIM is to approximate the nonlinear function f w in Eq. (15) using a separate set of basis vectors wheref is the coefficient of expansion of the nonlinear function f w in the reduced subspace spanned by {v i } m i=1 , V m ∈ R n×m is the matrix containing the first m columns of the left singular matrix V ∈ R n×n obtained from the SVD of the snapshot matrix X f of the nonlinear function f w . We note that no additional computational costs are associated with collecting the snapshot matrix of the nonlinear terms X f 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 (U r s B V m ) ∈ R r ×m in Eq. (17) is precomputed before solving Eq. (15). The overdetermined systemf = V m f w is approximated using the DEIM algorithm introduced in Chaturantabut and Sorensen (2010) by first computing a matrix P ∈ R n×m that selects m rows of the matrix V m to obtainf as follows: Using this expression off to approximate the nonlinear function in Eq. (17), we obtain a nonlinear term that is independent of n that takes the form where the matrix D = U r s B V m (P V m ) −1 ∈ R r ×m termed as the DEIM matrix. Similarly, the Jacobian of the nonlinear term in Eq. (15) is approximated using DEIM as follows: whereĴ f (f w (P U r sỹ s )) ∈ R m×m is the Jacobian matrix computed using the m components of f w evaluated by the DEIM algorithm (Chaturantabut and Sorensen 2010;Rewienski and White 2003;Nagoor Kani and Elsheikh 2017). Finally, the POD-DEIM-based reduced model takes the form dỹ s dt We note that POD-DEIM formulation is independent of the full model dimension n and that the DEIM procedure exploits the structure of the nonlinear function f w as component-wise operation at U r sỹ s (Chaturantabut and Sorensen 2010).
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 O(r 3 ) operations for each saturation update (Nagoor Kani and Elsheikh 2017; Bertsekas 1999). Moreover, previous studies (He 2010(He , 2013 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 (Nagoor Kani and Elsheikh 2017) and formulate an accurate POD-DEIM reduced-order models. DR-RNN is a deep RNN architecture (Nagoor Kani and Elsheikh 2017), constructed by stacking K physicsaware network layers. DR-RNN could be applied to any nonlinear dynamical system of the form dy dt = A y + F(y) where y(a, t) ∈ R n is the state variable at time t, a ∈ R d is a system parameter vector, the matrix A ∈ R n×n is the linear part of the dynamical system, and the vector F(y) ∈ R n is the nonlinear term (Nagoor Kani and Elsheikh 2017). The state variable y(t) at different time steps is obtained by solving the nonlinear residual equation defined as where r(t) is termed as the residual vector at time step t and y(t + 1) is the approximate solution of Eq. (22) at time step t + 1 obtained by using implicit Euler time integration method. DR-RNN (Nagoor Kani and Elsheikh 2017) approximates the solution of Eq. (22) using the following iterative update equations where U, w, η k are the training parameters of DR-RNN, φ h is the tanh activation function, • is an element-wise multiplication operator, r (k) t+1 is the residual in layer k obtained by substituting y t+1 = y (k−1) t+1 into Eq. (23), and G k 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 (Nagoor Kani and Elsheikh 2017). In this formulation, we set y (k=0) t+1 = y t . The architecture of DR-RNN is inspired by the rmsprop algorithm (Tieleman and Hinton 2012) 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 K 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 (r t+1 = −ỹ s t+1 +ỹ s t + D f w (P U r sỹ s t+1 ) +d). In such setting, the DR-RNN dynamics has a fixed computational budget of O(r 2 ) for each time step. In addition, DR-RNN has the prospect of employing large time step violating the numerical stability constraint (Nagoor Kani and Elsheikh 2017). Furthermore, DR-RNN does not rely on the reduced Jacobian matrix to approximate the solution of POD-DEIM reduced model.
The DR-RNN parameters θ = {U, w, η k } are fitted by minimizing the mean square error (mse) defined by where J MSE (mse) is the average distance between the reference solution y t and the RNN output y RNN t across a number of samples L with time-dependent observations (t = 1 . . . T and = 1 . . . L) (Nagoor Kani and Elsheikh 2017;Pascanu et al. 2013b). The set of parameters θ is commonly estimated by a technique called backpropagation through time (BPTT) (Werbos 1990;Rumelhart et al. 1986;Pascanu et al. 2013a;Mikolov et al. 2014), which backpropagates the gradient of the loss function J MSE with respect to θ in time over the length of the simulation.

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 p (DR-RNN combined with POD-Galerkin) and DR-RNN pd (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 (Chaturantabut and Sorensen 2010;Nagoor Kani and Elsheikh 2017;Chaturantabut and Sorensen 2010).
The outline of this section is as follows: In Sect. 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 Sect. 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 Sect. 5.3. In Sect. 5.4, we provide a set of error metrics utilized to evaluate the performance of the different ROMs. In Sect. 5.5, we present the numerical results for the quarter five spot model followed by results for the uniform flow model in the Sect. 5.6.

Full-Order Model Setup
We consider a two-phase (oil and water) porous media flow problem over the two-dimensional domain (4)]. The relative permeability is defined as a function of saturation using Corey's model k r w (s) = s * 2 , k ro = (1 − s * ) 2 , where s * = (s − s wc )/(1 − s or − s wc ), s wc is the irreducible water saturation, and s or is the residual oil saturation (Aarnes et al. 2007). We set s or = 0.2 and s wc = 0.2. We set the initial water saturation over the domain to the irreducible water saturation s wc = 0.2. 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 σ k is the variance, L k is the correlation length. In all test cases, we set σ k to 1 and the correlation length L k to 0.1 m. 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 64×64 blocks. The pressure equation is discretized using simple finite-volume method (aka. two-point flux approximation) (Aarnes et al. 2007), 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 0.015, and the total number of time steps is set to 160. 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 eighth saturation time step. For reference solutions, this system of equations is solved for 2000 random permeability realizations to estimate an ensemble-based statistics using Monte Carlo method (Ibrahima 2016).

POD-Galerkin-Based Reduced Model Setup
The first step in formulating POD reduced model is to compute the optimal POD basis matrices U r p and U r s . 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 2000 random permeability realizations into 45 clusters (Ghasemi 2015). 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, nonlin-ear 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 Sect. 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 2000 permeability realizations to estimate an ensemble-based statistics in the engineering quantities of interest.

DR-RNN Setup
In all the numerical test cases, we utilize DR-RNN with six layers [K = 6 in Eq. (24)]. We evaluate DR-RNN p and DR-RNN pd for different number of POD basis; however, we fix the number of DEIM basis to 35. The PyTorch framework (Paszke et al. 2017), a deep learning python package using Torch library as a backend, is used to implement the DR-RNN. Further, we optimize the DR-RNN model parameters using rmsprop algorithm (Tieleman and Hinton 2012;Paszke et al. 2017) as implemented in PyTorch, where we set the weighted average parameter to 0.9 and the learning rate to 0.001. The weight matrix U in Eq. (24) is initialized randomly from the uniform distribution function U[0.01, 0.02]. The vector training parameter w in Eq. (24) is initialized randomly from the uniform distribution function U[0.1, 0.5]. The scalar training parameters η k in Eq. (24) are initialized randomly from the uniform distribution U[0.1, 0.4]. We set the hyperparameters ζ and γ in Eq. (25) to 0.9 and 0.1, respectively. The formulated DR-RNN p and DR-RNN pd 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ỹ * s of the saturation state vector using y * s = U r s y s . 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).

Evaluation Metrics
We evaluate the performance of DR-RNN p and DR-RNN pd using two time specific error metrics defined by where l is the realization index, and y (RM) t is computed from the reduced model. Additionally, we utilize two relative error metrics defined as

Numerical Test Case 1
In this test case, water is injected at the lower left corner (0, 0) of the domain and a mixture of oil and water is produced at the top right corner of the domain (1, 1). We set the injection rate q = 0.05 at (0, 0) and q = −0.05 at (1, 1) 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 (r = 10, 20). The configuration of the problem domain is shown in top left panel of Fig. 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 X p in the top right panel, the saturation snapshot matrix X s in the bottom left panel, and the nonlinear function snapshot matrix X f in the bottom right panel.
The mean water saturation plots over the simulation time are shown in Fig. 3, where the results in the top row correspond to using 10 POD basis and the results in the bottom row correspond to using 20 POD basis. The subplots in Fig. 3 are arranged from left to right following the numbering of the spatial points shown in Fig. 2. From these results, it is clear that DR-DR-RNN p and DR-RNN pd results are very close to the least-square solutions (LS fit). In Fig. 3, POD-Galerkin reduced model yields extremely inaccurate and unstable results. We attribute the small errors in DR-RNN p and DR-RNN pd 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.
Figures 4, 5, and 6 show the results for the first (mean) and second (standard deviation) moments of the saturation field at time = 0.3 PVI obtained from the full model and from the various ROMs. In these Figs. 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 Fig. 4, the mean saturation obtained from DR-RNN ROMs is almost indistinguishable from the reference results. However, the mean saturation field obtained from POD reduced model (left panels of Fig. 6) deviates significantly from the reference mean saturation.
In Fig. 5, we observe small discrepancy of standard deviation results obtained in the DR-RNN ROMs in comparison with 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 Fig. 6 correspond to out of limits shown in colorbar.  Figure 7 settings are similar to the one adopted in Fig. 3. In Fig. 7, we can see that all the plots obtained from DR-DR-RNN p and DR-RNN pd 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 p and DR-RNN pd follows 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 Fig. 7, POD reduced model yields very inaccurate approximation of the saturation PDF irrespective of the number of POD basis. Figures 8 and 9 display samples of log(L 2 l,t ) and log(L ∞ l,t ) errors at time 0.3 PVI obtained from all the ROMs. All the ROMs use 10 POD basis to display the errors in Fig. 8 and likewise 20 POD basis to display the errors in Fig. 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 p and DR-RNN pd are nearly indistinguishable from the least-squares projection errors.
We further list in Table 1, the L rel 2 and L rel 2,max errors for the saturation field. From Table 1, we can see that the approximation errors obtained from DR-RNN p and DR-RNN pd 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 decrease 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 He (2010), 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 He (2010) 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.

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 (r = 10, 20). 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 X p , X s , and X f . Figure 11 shows the time plot of mean water saturation obtained from all the ROMs and from the full model. The display settings in Fig. 11 are the same as defined in Fig. 3. In Fig. 11, we can see that all the results obtained from DR-RNN p , DR-RNN pd , and the LS fit (the best approximation) closely approximate the full model whereas POD reduced model yields extremely inaccurate results regardless of the number of utilized POD basis. Figures 12, 13, and 14 show the results for the mean and standard deviation of the saturation field at 0.4 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 (Fig. 14) exhibit significant discrepancy when compared to the plots shown in Fig. 12. Again, we note that the white spots displayed in Fig. 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 p , DR-RNN pd 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, Fig. 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.
We further list in Table 2, the error metrics L rel 2 and L rel 2,max for the saturation fields. From Table 2, we can see that the approximation errors obtained from DR-RNN p and DR-RNN pd are almost close to the least-squares (best approximation) approximation errors. However, the POD reduced model yields very inaccurate results due to numerical instabilities.

Conclusion
In this work, we extended the DR-RNN introduced in Nagoor Kani and Elsheikh (2017) 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 pd . We evaluated the proposed DR-RNN pd on two forward uncertainty quantification problems involving twophase 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 2000 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 pd architecture employs explicit time stepping procedure for the same step size used in full model and POD reduced model. Hence, DR-RNN pd had a limited computational complexity O(K × r 2 ) instead of O( p × r 3 ) per saturation update, where r is the dimension of the POD reduced model, K p is the number of stacked network layers in DR-RNN and p is the average number of Newton iterations used in the standard POD-DEIM reduced model. The obtained numerical results show that DR-RNN pd provides accurate and stable approximations of the full model in comparison with the standard POD reduced model. Future work should consider the development of accurate and stable DR-RNN pd 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 pd for UQ tasks with the permeability fields that has randomly oriented channels or barriers. The use of DR-RNN pd for history matching (Elsheikh et al. 2012, 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. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.