We present a software framework tailored to seismic applications embedding high performance features. Especially, we discuss the spatial cache-blocking algorithm for the Finite-Difference Time-Domain (FDTD) method and the protocol to find optimal parameters. We apply our modeling engine to perform 3D elastic seismic modeling on a large-scale 877 km2 land survey. The velocity and density models are built with geological stratigraphic forward modeling. In total, 162300 shot gathers are computed with a maximum frequency of 20Hz and a maximum offset of 5 km. The main trait of the computed data is the complexity related to the surface-waves, especially their long and dispersive wave-train obscuring other wave arrivals. Computations were performed on the Kaust supercomputer Shaheen II. The modeling campaign took one month to complete using in average 2434 compute nodes in parallel. This achievement represents a workload of 59.6 ExaFLOP on a current PetaFLOP/s machine. This indicates that the next generation of supercomputers targeting the ExaFLOP/s sustained performance, would allow reducing the run-time of our application to one hour or less. With such performance, it is reasonable to predict that 3D elastic imaging will be a routinely used algorithm by seismic exploration in the years to come, while nowadays it still requires leading-edge hardware.