A hybridized discontinuous Galerkin method with high order spectral element discretization is presented for modeling multiphase immiscible 2D fluid flow in petroleum reservoir simulations. The time dependent nonlinear governing equations for reservoir simulation are presented, and a fully coupled implicit time stepping procedure is implemented using a high order W-Rosenbrock algorithm. For incompressible fluids considered in the present study, the governing PDE's include an elliptic PDE for the pressure field and a parabolic convective diffusion equation for the saturation (water volume fraction) in a two phase oil-water system. The PDE systems are coupled through the highly nonlinear relative permeability and capillary pressure functions relating the fluid velocities of individual species to the individual species pressures and saturations. With HDG/spectral element (HDGSEM) spatial discretization of the coupled governing equations, the system constitutes an index 1 differential algebraic equation. It is shown that the W-Rosenbrock algorithm in combination with the HDGSEM method yields an efficient and effective solution procedure which demonstrates robust, high order temporal accuracy.Numerous physical tests cases are presented demonstrating that the versatility of this approach is sufficient to address the challenges in modeling heterogeneous permeability fields and geologic features of arbitrary shape with discontinuous permeability functions. It is also capable of resolving the sharp time dependent displacement fronts in oil-water systems encountered in secondary oil recovery. It is demonstrated that the high order polynomial basis functions (up to 12th order) used in the HDGSEM algorithm have the ability to resolve complex geometry and achieve accurate results utilizing a simple Cartesian mesh for non-conforming geological features. Eliminating the need to re-grid and conform the mesh to the boundaries of different geological features greatly simplifies the workflow for petroleum engineers.