  Underground experiments searching for rare events, such as interactions from
dark matter, need to exhibit background as low as possible. One source of
background is from cosmic ray muons and muon-induced neutron production.
Presently these background are not fully understood. In this study Geant4 is
used to model cosmic ray muon induced neutron multiplicity production and
compare the modeling with data collected using an $^3$He instrumented Pb-target
detector system. The neutron event multiplicity production is taken from the
2002 NMDS-II data sets, consisting of 6504 hrs collected at 583 m.w.e. and 1440
hrs, with the identical detector system, collected at 1166 m.w.e.. The detector
consists of a 30 cm cube Pb-target surrounded by 60 $^3$He tubes. The single
particle detection efficiency is 23.2\%$\pm$1.2\% calibrated using a $^{252}$Cf
neutron source. The highest neutron multiplicity event, observed at 583 m.w.e.
was 54 observed neutrons corresponding to $\sim$ 233 produced neutrons. The
neutron multiplicity, n, distributions fit well a 2-parameter power law fit,
$k\times n^{-p}$. For the Monte Carlo simulations at both depths and the data
collected at both depths, all are consistent with a single slope parameter p.
For the simulation at 583 m.w.e., p=2.37$\pm0.01$ and for the data collected at
583 m.w.e, p=2.36$\pm0.10$. At 1166 m.w.e., p=2.31$\pm0.01$ for the simulation,
and for the data with only 6 detected events above multiplicity 5, p=$2.50 \pm
0.35$ predicted using a Maximum Likelihood Estimation method. At both depths,
the power law amplitudes of the Geant4 simulations differ by a factor of 2
larger than the data sets. However, the disagreement is within the estimated
systematic error of the simulations.
