SUMMARY We derive two sets of new wave equations involving the fractional time derivatives or fractional Laplacians for simulating seismic wave propagation in viscoelastic anisotropic (VA) media based on the Kjartansson's constant-Q model. The approximate fractional Laplacian wave equation is developed under the assumptions of the small velocity anisotropy parameter $\delta$ and attenuation anisotropy parameter ${\delta _Q}$ (or weak velocity and attenuation anisotropy). Both the formulas have the advantages of simple form and can accurately describe the constant-Q (i.e. frequency-independent quality factor) attenuation and arbitrary attenuation anisotropy behaviours compared to the widely used VA theory based on the generalized standard linear solids (GSLS) model and memory variables. Under the assumption of homogeneous plane wave and the constant-Q attenuation mechanism, we further derive exact analytical expressions for the phase velocities, attenuation coefficients, quality factors of the P- and SV-waves, and analyse their direction dependence in different attenuation anisotropy cases. For numerical modelling, we implement the fractional finite-difference (FD) method with the Grünwald–Letnikov (GL) approximation to solve the fractional time wave equation, and the generalized Fourier pseudospectral (PS) method to solve the fractional Laplacian wave equation, respectively. The PS method is highly efficient compared with the fractional FD method since it avoids additional memory to store the past wavefields. Numerical results of the homogeneous VTI (transversely isotropic with a vertical symmetry axis) model validate the accuracy of the two proposed schemes and illustrate the influence of attenuation and attenuation anisotropy on seismic wavefields, which are consistent with theoretical analysis. Finally, the modelling of the Hess VTI model shows the applicability of our formulations and algorithms in heterogeneous media.