This study presents a methodology to integrate two simulation codes, i.e., cohesive element-based numerical manifold method (Co-NMM) and unified pipe network method (UPM), for 2D coupled two-phase seepage-stress analysis of waterflooding processes in fractured rock masses. The Co-NMM is an extension of the numerical manifold method (NMM) and designed for efficiently addressing continuity-discontinuity problems of rock especially complex multiple-fracture propagation, while the UPM is a relatively novel approach for underground seepage analysis with single−/multiphase fluid flow in both fractured and fractured porous reservoirs. To link the Co-NMM and UPM properly, a two-phase seepage-stress coupling model is proposed for considering interactions between the fractured rock mass and two-phase fluids. The mechanical behaviour of the fractured rock mass responds to the variation in skeletal effective stress induced by the average two-phase fluid pressure in both the rock matrix and fractures, including deformation of the rock matrix, change in fracture aperture and the initiation and propagation of fractures. Accordingly, the induced changes in hydraulic behaviours such as the porosity, permeability and capillary pressure are captured by nonlinear stress-dependent coupling functions. Furthermore, an implicit sequential solution is adopted to execute the two solvers, consisting of two iteration loops, i.e., the two-phase flow solving loop and seepage-stress coupling solving loop. The relative performance of the developed 2D hybrid NMM-UPM method is illustrated for four cases related to the waterflooding process in porous and fractured porous mediums.