Fault zones are composed of two structural domains: the fault core, which accommodates most of the deformation, and a damage zone, with less intense deformation. The damage zone may act as a preferential flow path due to the presence of fractures, or as a barrier due to deformation bands. Therefore, the characterization of geological fault zones is essential for the adoption of adequate production strategies in oil fields. Generally, geophysical methods are used to characterize geological faults in the field. However, they hardly allow the identification of damage zones due to low seismic resolution. Alternatively, damage zones are analyzed through surface outcrops. Nonetheless, there is a wide dispersion of data in this type of study, which may be related to various factors, such as the properties of the host rock and the acting deformation mechanisms. Therefore, it is interesting to carry out this type of study in conjunction with numerical modeling to understand better the damage zone formation process. In this study, we present two methodologies based on the finite element method (FEM) to analyze the formation and evolution of damage zones at a reservoir scale. In the first methodology, the entire fault zone is represented through a continuum medium, while in the second methodology, the fault core is represented as a plane through a discontinuity. In both approaches, the damage zone is defined through the regions where plastic deformations were triggered. The numerical results obtained were close to field observations. They enabled the identification of the advantages and limitations of the two approaches based on the MEF. Finally, the results also allowed to identify the main parameters that influence the development of the damage zones and the different deformation mechanisms that occur along the damage zone.