In the last years, numerous seismological evidences have shown a strict correlation between fluid injection and seismicity. An important topic that is currently under discussion in the scientific community concerns the prediction of the earthquake magnitude that may be triggered by fluid injection activities. Coupled fluid flow and geomechanical deformation models can aim at understanding the evolution of pore pressure and rock deformation due to fluid injection in the subsurface. To perform an accurate numerical study of the correlation among fluid injection, seismicity rates and maximum earthquake magnitude, it is necessary to characterize the model with two fundamental features: first, the presence of a system of faults possibly intersecting among each other; second, the variability of the hydro-mechanical properties across the region surrounding each fault plane (fault zone). The novelty of this work is to account for these two aspects combining together two different numerical techniques that have been proposed in literature for the fault's modelling: for the first feature, interface elements are used to describe the frictional contacts occurring on the fault surfaces; for the second feature, solid elements are adopted to describe the heterogeneous hydro-mechanical behavior across the fault zone. Moreover, we account for a spatial variation of the permeability in the fault zone both along the dip and the normal direction with respect to the fault plane. We compute the numerical solution for six models among which we vary the permeability contrasts between protolith rocks and damage zone and between damage zone and fault core. We demonstrate that the anisotropy of permeability of the fault zone has a strong impact both on the timing and on the magnitude of triggered earthquakes. We suggest that a similar approach, which includes the entire architecture of the fault zone, shall be included in fluid-flow-geomechanical simulations to improve fault stability analysis during fluid injection.