Subsurface flow modeling is an indispensable task for reservoir management, but the associated computational cost is burdensome owing to model complexity and the fact that many simulation runs are required for its applications such as production optimization, uncertainty quantification, and history matching. To relieve the computational burden in reservoir flow modeling, a reduced-order modeling procedure based on hyper-reduction is presented. The procedure consists of three components: state reduction, constraint reduction, and nonlinearity treatment. State reduction based on proper orthogonal decomposition (POD) is considered, and the impact of state reduction, with different strategies for collecting snapshots, on accuracy and predictability is investigated. Petrov- Galerkin projection is used for constraint reduction, and a hyper-reduction that couples the Petrov-Galerkin projection and a 'gappy' reconstruction is applied for the nonlinearity treatment. The hyper-reduction method is a Gauss-Newton framework with approximated tensors (GNAT), and the main contribution of this study is the presentation of a procedure for applying the method to subsurface flow simulation. A fully implicit oil-water two-phase subsurface flow model in three-dimensional space is considered, and the application of the proposed hyper-reduced-order modeling procedure achieves a runtime speedup of more than 300 relative to the full-order method, which cannot be achieved when only constraint reduction is adopted. In addition, two types of sequential Bayesian filtering for history matching are considered to investigate the performance of the developed hyper-reduced-order model to relive the associated computational cost. First, an ensemble Kalman filter (EnKF) is considered for Gaussian system and a procedure embedding the hyper-reduced model (HRM) into the EnKF is presented. The use of the HRM for the EnKF significantly reduces the computational cost without much loss of accuracy, but the combination requires a few remedies such as clustering to find an optimum reduced-order model according to spatial similarity of geological condition, which causes an additional computation. For non-Gaussian system, an advanced particle filter, known as regularized particle filter (RPF), is considered because it does not take any distributional assumptions. Particle filtering has rarely been applied for reservoir history matching due to the fact that it is hard to locate the initial particles on highly probable regions of state spaces especially when large scale system is considered, which makes the required number of particles scale exponentially with the model dimension. To resolve the issues, reparameterization is adopted to reduce the order of the geological parameters. For the reparameterization, principal component analysis (PCA) is used to compute the reduced space of the model parameters, and by constraining the filtering analysis with the computed subspace the required number of initial particles can be reduced down to a manageable level. Consequently, a huge computational saving is achieved by embedding the HRM into the RPF. Furthermore, the additional cost of clustering required to identify the geospatially optimum reduced-order model is saved because the advanced particle filter allows to easily identify the groups of geospatially similar particles.