This paper develops a fractional stochastic partial differential equation (SPDE) to model the evolution of a random tangent vector field on the unit sphere. The SPDE is governed by a fractional diffusion operator to model the Lévy-type behaviour of the spatial solution, a fractional derivative in time to depict the intermittency of its temporal solution, and is driven by vector-valued fractional Brownian motion on the unit sphere to characterize its temporal long-range dependence. The solution to the SPDE is presented in the form of the Karhunen-Loève expansion in terms of vector spherical harmonics. Its covariance matrix function is established as a tensor field on the unit sphere that is an expansion of Legendre tensor kernels. Approximations to the solutions are studied and convergence rates of the approximation errors are given. It is demonstrated how these convergence rates depend on the decay of the power spectrum and variances of the fractional Brownian motion.