  Motivated by the recent detection of a large number of embedded young stellar
objects (YSOs) with mass accretion rates that are inconsistent with the
predictions of the standard model of inside-out collapse (Shu 1977), we perform
a series on numerical hydrodynamic simulations of the gravitational collapse of
molecular cloud cores with various initial masses, rotation rates, and sizes.
We focus on the early Class I stage of stellar evolution when circumstellar
disks are exposed to high rates of mass deposition from infalling envelopes.
Our numerical modeling reproduces the large observed spread in accretion rates
inferred for embedded YSOs in Perseus, Serpens, and Ophiuchus star forming
regions by Enoch et al. (2009), yielding 37%--75% of objects with "sub-Shu"
accretion rates \dot{M} \la 10^{-6} Msun/yr and 1%--2% of objects with
"super-Shu" accretion rates \dot{M}>10^{-5} Msun/yr. Mass accretion rates in
the Class I stage have a log-normal distribution, with its shape controlled by
disk viscosity and disk temperature. The spread in $\dot{M}$ is greater in
models with lower viscosity and smaller in models with higher viscosity and
higher disk temperature, suggesting that gravitational instability may be a
dominant cause of the observed diversity in $\dot{M}$ in embedded YSOs. Our
modeling predicts a weak dependence between the mean mass accretion rates and
stellar masses in the Class I stage, in sharp contrast to the corresponding
steep dependence for evolved T Tauri stars and brown dwarfs.
