The light scattering properties of particles play important roles in radiative transfer in many dispersed systems, such as turbid atmosphere, ocean water, nanofluids, composite coatings and so on. As one of the scattering property parameters, the scattering phase functions of particles are strongly dependent on the particle size, size distribution, and morphology, as well as on the complex refractive indices of the particles and surrounding media. For the sake of simplicity, the empirical phase function models are widely used in many practical applications. In this work, we focus on the radiative transfer problem in dispersed systems composed of spherical particles, and give quantitative analyses of the impact of scattering phase functions on the radiative transfer process. We fit the scattering phase functions of four different types of practical dispersed systems using four previously proposed empirical phase function models, including the Henyey–Greenstein (HG) model, Cornette Shanks (CS) model, Reynold and McCormick (RM) model and two-term Reynolds–McCormick (TTRM) model. By comparing the radiative transfer characteristics (i.e., hemispherical reflectance, hemispherical transmittance and total absorptance) of dispersed layers calculated using the Monte Carlo method, the relative errors caused by using the empirical phase functions are systematically investigated. The results demonstrate that the HG, CS and RM models cause obvious errors in the calculation of hemispherical reflectance in many cases. Meanwhile, the induced errors show no obvious regularity, but are related to the particle size and layer optical thickness. Due to the good fitting effect in both forward and backward directions, the TTRM model provides significantly higher performances in fitting the phase functions of all considered cases than the widely used single-term parametrizations. Moreover, for different particle sizes and layer optical thicknesses, the induced errors of the TTRM model in radiative transfer characteristics are very small, especially for the case of polydisperse particles. Our results can be used to guide the design, analysis and optimization of dispersed systems in practical optics and photonics applications.