Multiphase flow in fractured porous media is relevant in many geoscience applications, including hydrocarbon reservoirs, underground water resources and subsurface storage of fluids (including hydrogen and carbon dioxide). Depositional environments are geometrically complex anisotropic media, because the sedimentary layers can be deposited in different ways, giving different preferential directions to the fluid flow. Therefore, accurate and scalable modeling of fluid flow through these formations still remains as a great challenge. The presence of explicitly-modeled fractures with flowing conductivities ranging from sealing barriers to highly conductive channels significantly increases the simulation complexity. If these fractures are modeled through a conforming mesh, the imposed constraints on the mesh geometry makes it too difficult to be applicable to complex 3D highly-fractured media. On the other hand, if a non-conforming method is used, the fractures can cross the rock matrix grid cells, making it feasible to model field-relevant cases. In this context, the embedded discrete fracture model (EDFM) employs independent meshes for fractures and matrix domains. With the extension of the projection-based embedded discrete fracture model (pEDFM), additional non-neighboring connections are constructed along with modified matrix-matrix transmissibilities, allowing to capture all types of fractures from low to highly permeable. Just recently pEDFM has been extended to discrete models with Corner-Point-Grid geometries, which are still structured elements. In the present work, we develop a fully implicit strategy to numerically simulate two-phase immiscible flow in 3-D naturally fractured reservoirs by extending the pEDFM method to general unstructured tetrahedral meshes. To discretize the diffusive (pressure) terms, we devise a 3-D extension of the multipoint flux approximation approach that uses the "diamond stencil" (MPFA-D), which is very robust and can deal with full tensors on arbitrary tetrahedral meshes. While the advective terms are discretized through the classical first order upwind method (FOUM). The performance of the proposed method is investigated for several test cases. It is found that our approach is robust and capable to accurately capture the effects of both high and low permeability fractures for different unstructured meshes and heterogeneous and highly anisotropic permeability tensors.