A computer-implemented method of modeling a subsurface formation, comprising: computing an iso-surface for an iso-value from a three-dimensional stratigraphic function for a volume of interest in the subsurface formation; computing first and second strike traces (454, 456), being fault parallel, following a topography of the computed iso-surface on respective first and second sides of a fault in the volume of interest; extracting seismic data (458, 460) from a seismic volume in the volume of interest along the first and second strike traces (454, 456); correlating the extracted seismic data (458, 460) along the first and second strike traces (454, 456); and computing (402, 446) a three-dimensional fault displacement vector for the fault from the correlated extracted seismic data (458, 460) along the first and second strike traces (454, 456), wherein the fault displacement vector characterizes the net-slip of a horizon along the fault and includes one or more of a dip-slip component, a strike-slip component, vertical and horizontal throw components and a heave component and wherein said computing includes determining a magnitude of the fault displacement vector based upon an offset between matching signatures (462, 464) in the extracted seismic data (458. 460) along the first and second strike traces (454, 456); and computing a fault damage zone for the fault by performing the fault displacement vector calculation operation for a plurality of pairs of first and second strike traces, wherein each pair of strike traces is separated by a different distance to identify a distance proximate an edge of the fault damage zone. computing an iso-surface for an iso-value from a three-dimensional stratigraphic function for a volume of interest in the subsurface formation; computing first and second strike traces (454, 456), being fault parallel, following a topography of the computed iso-surface on respective first and second sides of a fault in the volume of interest; extracting seismic data (458, 460) from a seismic volume in the volume of interest along the first and second strike traces (454, 456); correlating the extracted seismic data (458, 460) along the first and second strike traces (454, 456); and computing (402, 446) a three-dimensional fault displacement vector for the fault from the correlated extracted seismic data (458, 460) along the first and second strike traces (454, 456), wherein the fault displacement vector characterizes the net-slip of a horizon along the fault and includes one or more of a dip-slip component, a strike-slip component, vertical and horizontal throw components and a heave component and wherein said computing includes determining a magnitude of the fault displacement vector based upon an offset between matching signatures (462, 464) in the extracted seismic data (458. 460) along the first and second strike traces (454, 456); and computing a fault damage zone for the fault by performing the fault displacement vector calculation operation for a plurality of pairs of first and second strike traces, wherein each pair of strike traces is separated by a different distance to identify a distance proximate an edge of the fault damage zone.