Despite great advances in reservoir simulation capabilities with the introduction of high-performance computing (HPC) platforms and enhanced solvers, high fidelity grid-based simulation still remains a challenging task. This task is especially demanding for fine-resolved geological reservoirs with multiphase and multicomponents and in production optimization and uncertainty quantification frameworks where several calls of the large scale simulation model need to be performed. Model order reduction techniques have been applied in porous media flow modeling in order to obtain fast simulation surrogates to alleviate the high computational cost associated with these simulations. The difficulty of applying model reduction techniques to porous media flow arises from the nonlinearity of the system. In order to overcome the issue of nonlinearity, we introduce the bilinear form of the dynamical system which in many cases produces a satisfying approximation of the system. The bilinear approximation is a simple form of the parent system and it is linear in the input and linear in the state but it not linear in both jointly. The bilinear form is computed using truncated multidimensional Taylor's expansion of the system terms. Examples are presented to illustrate this recent approach for the case of single phase flow modeling, and comparisons are made with the case of linearized models and the full nonlinear models. In addition, we present methods to reduce the complexity of the bilinear system of equations and results are presented that compare linear model reduction (balanced truncation) and nonlinear model reductions, such as POD techniques.