The accurate modeling of elastic wave propagation in TTI media is significantly important for the study of the anisotropy, amplitude versus offset (AVO) analysis, and other inverse problems. A new numerical algorithm named parallel weighted nearly-analytical discrete (WNAD) method is presented in this paper. As an extension of 2D WNAD, the parallel WNAD method is a fast and efficient algorithm for simulating the elastic wave propagation in TTI media. We first show how to formulate the 3D WNAD method, and then compare the wavefield generated by this method with the wavefield computed by a conventional staggered-grid method. The numerical results show that the WNAD algorithm is more suitable for the simulations of the large-scale seismic wavefield by using the high-performance computers. Using the parallel WNAD algorithm, we study the elastic wave propagation in the TTI media and observe the important feature of TI media: shear-wave splitting, the quasi body wave coupling and velocity anisotropy. The reflected, refracted, and converted waves are generated at the interfaces. This makes the wavefield complex. To better understand the anisotropic media induced by fracture, it’s useful to study and identify those waves.