Simulation of subsurface flows in porous media is difficult due to the nonlinearity of the operators and the high heterogeneity of material coefficients. In this paper, we present a scalable fully implicit solver for incompressible two-phase flows based on overlapping domain decomposition methods. Specifically, an inexact Newton-Krylov algorithm with analytic Jacobian is used to solve the nonlinear systems arising from the discontinuous Galerkin discretization of the governing equations on 3D unstructured grids. The linear Jacobian system is preconditioned by additive Schwarz algorithms, which are naturally suitable for parallel computing. We propose a hybrid two-level version of the additive Schwarz preconditioner consisting of a nested coarse space to improve the robustness and scalability of the classical one-level version. On the coarse level, a smaller linear system arising from the same discretization of the problem on a coarse grid is solved by using GMRES with a one-level preconditioner until a relative tolerance is reached. Numerical experiments are presented to demonstrate the effectiveness and efficiency of the proposed solver for 3D heterogeneous medium problems. We also report the parallel scalability of the proposed algorithms on a supercomputer with up to 8,192 processor cores.