We apply Massively Parallel Interface for MPFA-O scheme with state-of-the-art Operator-Based Linearization (OBL) approach for multiphase flow in porous media. The implementation of MPFA-O scheme enhances the modelling capabilities for non-K-orthogonal mesh. A fully implicit scheme is applied to guarantee the stability of solutions when a mass-based formulation is involved to keep the flexibility of the framework for general-purpose reservoir simulation. As the MPFA-O introduces more non-zeros elements in the Jacobian matrix than the traditional TPFA, massively parallel computations via Message Passing Interface (MPI) in this work help to guarantee competitive computational efficiency for high-fidelity geological models. Concerning the Jacobian assembly hassle, we apply the OBL approach which introduces operators combining the fluid and rock properties in the conservation equations and discretizes the operators in the physical parameter space. By computing values and derivatives of the operators via a multilinear interpolation, the assembly of Jacobian matrix and residual vector could be drastically simplified. Another benefit of the OBL is that by only evaluating operator values on the predefined nodes in the physical parameter space, the overhead related to complex phase behavior and property evaluation is significantly reduced. In the end, we present several benchmark cases to rigorously demonstrate the accuracy, convergence, and robustness of the framework and two challenging field-scale cases to further prove its computing performance and parallel scalability.