A computer-implemented method for subsurface characterization from seismic data, the method comprising: a. receiving, at a computer processor, effective parameters and properties from an anisotropic seismic velocity model, at least one velocity function, and key descriptors from seismic data, wherein the effective parameters and properties used comprise an effective moveout velocity V nmo and an effective η, η eff, with η representing a difference between the horizontal velocity of the seismic energy and the near offset moveout velocity v nmo of the seismic energy, wherein the at least one velocity function comprises the vertical velocity v o as a function of depth at a particular spatial location, and wherein the key descriptors used comprise a plurality of maximum usable frequencies, a minimum offset, a maximum offset, and an angle mute; b. determining, via a computer processor, a detectability threshold for moveout in a seismic data gather based on the seismic data, wherein the detectability threshold is determined based on the key descriptors of the seismic data and the effective parameters and properties of the velocity model; and c. computing, via a computer processor, a depth uncertainty function Δz(z) for the particular spatial location based on the detectability threshold, the effective parameters and properties, and the at least one velocity function, wherein the depth uncertainty function Δz(z) expresses depth uncertainty Δz as a function of depth z and represents an error estimate that is used to analyze an interpretation of the seismic data, wherein the computing of the depth uncertainty function Δz(z) for the particular spatial location comprises generating bounding effective moveout velocity and η curves, (V nmoLO, η effLO) and (V nmoHI, η effHI), for a plurality of zero-offset intercept times t o using the detectability threshold to obtain vertical velocities v oLO and v oHI as functions of t o for the low and high bounds of detectability, computing the depth functions z LO (t o) and z HI (t o) from v oLO and v oHI by integrating over to and obtaining therefrom the interval depth models z LO (z) and z HI (z) as functions of depth, and computing the depth uncertainty function Δz(z) using Δz(z) = z HI (z) - z LO (z). a. receiving, at a computer processor, effective parameters and properties from an anisotropic seismic velocity model, at least one velocity function, and key descriptors from seismic data, wherein the effective parameters and properties used comprise an effective moveout velocity V nmo and an effective η, η eff, with η representing a difference between the horizontal velocity of the seismic energy and the near offset moveout velocity v nmo of the seismic energy, wherein the at least one velocity function comprises the vertical velocity v o as a function of depth at a particular spatial location, and wherein the key descriptors used comprise a plurality of maximum usable frequencies, a minimum offset, a maximum offset, and an angle mute; b. determining, via a computer processor, a detectability threshold for moveout in a seismic data gather based on the seismic data, wherein the detectability threshold is determined based on the key descriptors of the seismic data and the effective parameters and properties of the velocity model; and c. computing, via a computer processor, a depth uncertainty function Δz(z) for the particular spatial location based on the detectability threshold, the effective parameters and properties, and the at least one velocity function, wherein the depth uncertainty function Δz(z) expresses depth uncertainty Δz as a function of depth z and represents an error estimate that is used to analyze an interpretation of the seismic data, wherein the computing of the depth uncertainty function Δz(z) for the particular spatial location comprises generating bounding effective moveout velocity and η curves, (V nmoLO, η effLO) and (V nmoHI, η effHI), for a plurality of zero-offset intercept times t o using the detectability threshold to obtain vertical velocities v oLO and v oHI as functions of t o for the low and high bounds of detectability, computing the depth functions z LO (t o) and z HI (t o) from v oLO and v oHI by integrating over to and obtaining therefrom the interval depth models z LO (z) and z HI (z) as functions of depth, and computing the depth uncertainty function Δz(z) using Δz(z) = z HI (z) - z LO (z).