Introducing flexible covariance functions is critical for interpolating spatial data since the properties of interpolated surfaces depend on the covariance function used for Kriging. An extensive literature is devoted to covariance functions on Euclidean spaces, where the Matérn covariance family is a valid and flexible parametric family capable of controlling the smoothness of corresponding stochastic processes. Many applications in environmental statistics involve data located on spheres, where less is known about properties of covariance functions, and where the Matérn is not generally a valid model with great circle distance metric. In this paper, we advance the understanding of covariance functions on spheres by defining the notion of and proving a characterization theorem for m times mean square differentiable processes on d-dimensional spheres. Stochastic processes on spheres are commonly constructed by restricting processes on Euclidean spaces to spheres of lower dimension. We prove that the resulting sphere-restricted process retains its differentiability properties, which has the important implication that the Matérn family retains its full range of smoothness when applied to spheres so long as Euclidean distance is used. The restriction operation has been questioned for using Euclidean instead of great circle distance. To address this question, we construct several new covariance functions and compare them to the Matérn with Euclidean distance on the task of interpolating smooth and non-smooth datasets. The Matérn with Euclidean distance is not outperformed by the new covariance functions or the existing covariance functions, so we recommend using the Matérn with Euclidean distance due to the ease with which it can be computed.