We develop a Bayesian model to invert surface seismic refraction data with depth
constraints from boreholes for characterization of aquifer geometry and apply it to seismic
and borehole data sets collected at the contaminated Oak Ridge National Laboratory site
in Tennessee. Rather than the traditional approach of first inverting the seismic arrival times
for seismic velocity and then using that information to aid in the spatial interpolation
of wellbore data, we jointly invert seismic first arrival time data and wellbore‐based
information, such as depths of key lithological boundaries. We use a staggered‐grid finite
difference algorithm with second‐order accuracy in time and fourth‐order accuracy in
space to model seismic full waveforms and use an automated method to pick the first arrival
times. We use Markov Chain Monte Carlo methods to draw many samples from the joint
posterior probability distribution, on which we can estimate the key interfaces and their
associated uncertainty as a function of horizontal location and depth. We test the developed
method on both synthetic and field case studies. The synthetic studies show that the
developed method is effective at rigorous incorporation of multiscale data and the Bayesian
inversion reduces uncertainty in estimates of aquifer zonation. Applications of the approach
to field data, including two surface seismic profiles located 620 m apart from each other,
reveal the presence of a low‐velocity subsurface zone that is laterally persistent. This
geophysically defined feature is aligned with the plume axis, suggesting it may serve as
an important regional preferential flow pathway.