In the domain of reservoir engineering, fault stability assessment is traditionally performed by assimilating fault systems to surfaces. However, faults are complex and heterogeneous geological systems, whose compartmentalized architecture generally corresponds to an inner core surrounded by outer, often fractured, damage zones (DZ). To date, few studies have assessed the potential role of the DZ characteristics (anisotropy, complex spatial distribution, etc.) on the shear behavior of the whole system, the main challenge being the proper (numerical) integration of fracture networks in reservoir simulators. The present study was motivated by the stochastically generated fracture networks representative of DZ investigated on outcrops (Cirques de Navacelles, South of France). We propose here a methodology relying on homogenization techniques for deriving the spatial distribution of effective anisotropic hydro-poro-elastic DZ properties from these stochastic fracture networks. Using them as inputs of 2D finite-element coupled hydro-mechanical models of a reservoir-scale fault zone, the key role of the spatially-varying DZ stiffness and Biot's coefficients are highlighted: neglecting them might lead to a high over-estimation of the maximum sustainable injection pressure (by nearly 50% in the considered cases).