We present a reduced-order modeling technique for subsurface multi-phase flow 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 fixed 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-fidelity 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 quantification test cases using Monte Carlo simulation of subsurface flow with random permeability field. The obtained results demonstrate that DR-RNN combined with POD-DEIM provides an accurate and stable reduced model with a fixed computational budget that is much less than the computational cost of standard POD-Galerkin reduced model combined with DEIM for nonlinear dynamical systems.