Surface waves are modulated and scattered by topographic variations, rough interfaces, and medium heterogeneities with the propagation around the Earth. Scattered surface waves appear as wave trains between multiple direct arrivals in long-period seismograms. The single isotropic-scattering model for Rayleigh waves of the fundamental mode successfully explained observed mean square (MS) envelopes of vertical-component seismograms. However, it is necessary to consider multiple scattering process for more accurate simulation of envelopes for a wide range of lapse time. Extending the radiative transfer equation for a flat surface, we mathematically formulate the multiple isotropic-scattering process on a spherical surface. We can solve this integral equation using the Laplace transform in time and the spherical harmonics expansion in angular space. The time trace of energy density corresponding to the MS envelope can be obtained by using the inverse Laplace transform. Analyzing 12 IRIS station data of the 1999 Kocaeli, Turkey earthquake for 80--180 s periods, we estimated the total scattering coefficient g0 ≈ 2 ¿ 10-6 1/km and the total attenuation Q-1 ≈ 8.475 ¿ 10-3 for the fundamental-mode Rayleigh waves. The multiple scattering model explains observed envelopes better than the single scattering model; however, synthesized envelopes are still smaller than observed envelopes at lapse times larger than about 20,000 s. It suggests that we need to study scattering and dispersion of higher modes in addition to the multiple scattering of the fundamental mode especially at large lapse times, because higher-mode Rayleigh waves attenuate more slowly and spread out more rapidly compared with those of the fundamental mode.