  In this paper we study the depth at which a cosmic ray shower reaches its
maximum (\Xmax) as predicted by Monte Carlo simulation. The use of \Xmax in the
determination of the primary particle mass can only be done by comparing the
measured values with simulation predictions. For this reason it is important to
study the differences between the available simulation models. We have done a
study of the first and second moments of the \Xmax distribution using the
\Corsika and \Conex programs. The study was done with high statistics in the
energy range from $10^{17}$ to \energy{20.4}. We focus our analysis in the
different implementations of the hadronic interaction models \Sibyll and
\Qgsjet in \Corsika and \Conex. We show that the predictions of the \meanXmax
and \sigmaXmax depend slightly on the combination of simulation program and
hadronic interaction model. Although these differences are small, they are not
negligible in some cases (up to 5 g/cm$^2$ for the worse case) and they should
be considered as a systematic uncertainty of the model predictions for
\meanXmax and \sigmaXmax. We have included a table with the suggested
systematic uncertainties for the model predictions. Finally, we present a
parametrization of the \Xmax distribution as a function of mass and energy
according to the models \Sibyll and \Qgsjet, and showed an example of its
application to obtain the predicted \Xmax distributions from cosmic ray
propagation models.
