A novel approach to the automatic detection of fault surface images in 3D seismic datasets is presented. Based on the premise that seismic faulting introduces discontinuities into the rock layering (that is, the horizons), a coherency measure is used to detect points of significant horizon discontinuity. A highest confidence first (HCF) merging strategy is then combined with a flexible surface model to estimate the 3D fault surfaces iteratively.