This paper investigates solutions of hyperbolic diffusion equations in $\mathbb{R}^3$ with random initial conditions. The solutions are given as spatial-temporal random fields. Their restrictions to the unit sphere $S^2$ are studied. All assumptions are formulated in terms of the angular power spectrum or the spectral measure of the random initial conditions. Approximations to the exact solutions are given. Upper bounds for the mean-square convergence rates of the approximation fields are obtained. The smoothness properties of the exact solution and its approximation are also investigated. It is demonstrated that the Hölder-type continuity of the solution depends on the decay of the angular power spectrum. Conditions on the spectral measure of initial conditions that guarantee short or long-range dependence of the solutions are given. Numerical studies are presented to verify the theoretical findings.