In this paper, shear wave propagation in soils is examined in a stochastic context considering spatial variability of the shear modulus soil parameter. To this purpose, the recently established concept of dynamic mean and variability response functions (DMRF, DVRF) is reformulated in the framework of stochastic finite element analyses of shear wave propagation problems in order to efficiently calculate the response time history statistics of the soil surface. Similarly to the approximation formulas of classical VRFs, a fast Monte Carlo simulation procedure is implemented to numerically evaluate the above functions in the time domain. The main advantage of the proposed methodology lies on the independence of the DMRF and DVRF on the marginal probability density function and correlation structure of the stochastic system parameter, which in our case is assumed to be the inverse of the soil shear modulus 1/G. By integrating the product of the spectral density of 1/G with the DMRF and DVRF, the mean and variance of the ground response are obtained at each time step of the dynamic analysis. The method also allows for the estimation of time dependent but spectral and probability distribution free upper bounds of the response mean and variance. To illustrate the efficiency and applicability of the proposed approach, stochastic finite element analyses of wave propagation of a Ricker synthetic wavelet as well as a recorded earthquake motion in 1D and 2D soil domains are performed and a sensitivity analysis is carried out with respect to various correlation structures of the underlying random fields representing 1/G. The accuracy of the proposed methodology is validated with comparison to direct Monte Carlo simulation. Useful conclusions regarding the sensitivity of the system response to the spectral characteristics of the underlying random fields representing 1/G are drawn.