OpenMP/MPI hybrid parallel programming models were implemented to 3D finite-volume based simulation code for groundwater flow problems through heterogeneous porous media using parallel iterative solvers with multigrid preconditioning. Performance and robustness of the developed code has been evaluated on the "T2K Open Supercomputer (T2K/Tokyo)" using up to 8,192 cores through both of weak and strong scaling computations. Furthermore new strategy for solving equations at the coarsest level (coarse grid solver) was proposed and evaluated. New coarse grid solver improved the scalability of the multigrid solver dramatically. OpenMP/MPI hybrid parallel programming model, where one MPI process is applied to a single socket (HB 4x4) demonstrated the best performance and robustness for large-scale ill-conditioned problems with appropriate optimization and coarse grid solvers.