  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\'{e}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\`{e}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.
