Large footprint waveform-recording laser altimeters (lidars) have demonstrated a potential for accurate remote sensing of forest biomass and structure, important for regional and global climate studies. Currently, radiative transfer analyses of lidar data are based on the simplifying assumption that only single scattering contributes to the return signal, which may lead to errors in the modeling of the lower portions of recorded waveforms in the near-infrared spectrum. In this study we apply time-dependent stochastic radiative transfer (RT) theory to model the propagation of lidar pulses through forest canopies. A time-dependent stochastic RT equation is formulated and solved numerically. Such an approach describes multiple scattering events, allows for realistic representation of forest structure including foliage clumping and gaps, simulates off-nadir and multiangular observations, and has the potential to provide better approximations of return waveforms. The model was tested with field data from two conifer forest stands (southern old jack pine and southern old black spruce) in central Canada and two closed canopy deciduous forest stands (with overstory dominated by tulip poplar) in eastern Maryland. Model-simulated signals were compared with waveforms recorded by the Scanning Lidar Imager of Canopies by Echo Recovery (SLICER) over these regions. Model simulations show good agreement with SLICER signals having a slow decay of the waveform. The analysis of the effects of multiple scattering shows that multiply scattered photons magnify the amplitude of the reflected signal, especially that originating from the lower portions of the canopy.