Reservoir Geomechanics is playing an increasingly important role in developing and producing hydrocarbon reserves. One of the main challenges in reservoir modeling is accurate and efficient simulation of arbitrary intersecting faults. In this paper, we propose a new formulation to model multiple intersecting cohesive discontinuities (faults) in reservoirs using the XFEM framework. This formulation involves construction of enrichment functions and computation of stiffness sub-matrices for bulk (rock mass), faults and their interaction DOF’s. A sub-divisional scheme has been developed, to perform volume and surface integrations. The sub-divisional scheme can efficiently handle most possible configurations of multiple intersecting/non-intersecting faults in 3D finite elements. In order to be able to model cohesive behaviour of intersecting faults, we have performed additional surface integrations which, to the knowledge of the authors, are not reported in the literature. For example, if an element is intersected by two cohesive faults, we compute 16 volumes and 8 surface integrations for stiffness matrices and 4 surface integrations for force vectors, to form elemental equations. The current distributed implementation using C++/MPI can simulate very large meshes. Three-dimensional benchmark cases are presented to validate the accuracy of the approach, and the potential benefits in applications. The present study shows that interaction DOF’s and their associated surface stiffness matrices play a significant role in accurate modeling of cohesive faults via XFEM.