This paper examines the temporal evolution of a two-stage stochastic model for spherical random fields. The model uses a time-fractional stochastic hyperbolic diffusion equation, which describes the evolution of spherical random fields on $\bS^2$ in time. The diffusion operator incorporates a time-fractional derivative in the Caputo sense. In the first stage of the model, a homogeneous problem is considered, with an isotropic Gaussian random field on $\bS^2$ serving as the initial condition. In the second stage, the model transitions to an inhomogeneous problem driven by a time-delayed Brownian motion on $\bS^2$. The solution to the model is expressed through a series of real spherical harmonics. To obtain an approximation, the expansion of the solution is truncated at a certain degree $L\geq1$. The analysis of truncation errors reveals their convergence behavior, showing that convergence rates are affected by the decay of the angular power spectra of the driving noise and the initial condition. In addition, we investigate the sample properties of the stochastic solution, demonstrating that, under some conditions, there exists a local Hölder continuous modification of the solution. To illustrate the theoretical findings, numerical examples and simulations inspired by the cosmic microwave background (CMB) are presented.