The absorption of compressional and shear waves in many viscoelastic solids has been experimentally shown to follow a frequency power law. It is now well established that this type of loss behavior can be modeled using fractional derivatives. However, previous fractional constitutive equations for viscoelastic media are based on temporal fractional derivatives. These operators are non-local in time, which makes them difficult to compute in a memory efficient manner. Here, a fractional Kelvin-Voigt model is derived based on the fractional Laplacian. This is obtained by splitting the particle velocity into compressional and shear components using a dyadic wavenumber tensor. This allows the temporal fractional derivatives in the Kelvin-Voigt model to be replaced with spatial fractional derivatives using a lossless dispersion relation with the appropriate compressional or shear wave speed. The model is discretized using the Fourier collocation spectral method, which allows the fractional operators to be efficiently computed. The field splitting also allows the use of a k-space corrected finite difference scheme for time integration to minimize numerical dispersion. The absorption and dispersion behavior of the fractional Laplacian model is analyzed for both high and low loss materials. The accuracy and utility of the model is then demonstrated through several numerical experiments, including the transmission of focused ultrasound waves through the skull.