We present a workflow for geostatistical modelling of the faults and damage zones. The approach is based on the combination of the numerical simulation of geological faults formation using meshless Discrete Element Method with further estimation of the statistics of strains distribution in the damage zones. After that, this information is used to incorporate the damage zones to a structural surfaces representing faults (as results of conventional seismic interpretation). Further on the mechanical properties of the rocks in the damage zones are updated to construct the grid-based geological model for either seismic of hydrodynamic simulations.