The variations of mechanical properties in soft tissues are biomarkers used for clinical diagnosis and disease monitoring. Optical coherence elastography (OCE) has been extensively developed to investigate mechanical properties of various biological tissues. These methods are generally based on time-domain data and measure the time-of-flight of the localized shear wave propagations to estimate the group velocity. However, there is considerable information that can be obtained from examining the mechanical properties such as wave propagation velocities at different frequencies. Here we propose a method to evaluate phase velocity, wave velocity at various frequencies, in four-dimensional space (x, y, z, f), called 4D-OCE phase velocity. The method enables local estimates of the phase velocity of propagating mechanical waves in a medium. We acquired and analyzed data with this method from a homogeneous reference phantom, a heterogeneous phantom material with four different excitation cases, and ex vivo porcine kidney tissue. The 3D-OCE group velocity was also estimated to compare with 4D-OCE phase velocity. Moreover, we performed numerical simulation of wave propagations to illustrate the boundary behavior of the propagating waves. The proposed 4D-OCE phase velocity is capable of providing further information in OCE to better understand the spatial variation of mechanical properties of various biological tissues with respect to frequency.