We address the problem of statistical simulation of a scalar real Gaussian random field inside the unit 3D ball. Two different methods are studied: (i) the method based on the known homogeneous isotropic power spectrum developed by Meschede and Romanowicz [M. Meschede and B. Romanowicz, Non-stationary spherical random media and their effect on long-period mantle waves, Geophys. J. Int. 203 2015, 1605-1625] and (ii) the method based on known radial and angular covariance functions suggested in this work. The first approach allows the extension of the simulation technique to the inhomogeneous or anisotropic case. However, the disadvantage of this approach is the lack of accurate statistical characterization of the results. The accuracy of considered methods is illustrated by numerical tests, including a comparison of the estimated and analytical covariance functions. These methods can be used in many applications in geophysics, geodynamics, or planetary science where the objective is to construct spatial realizations of 3D random fields based on a statistical analysis of observations collected on the sphere or within a spherical region.