Accurate wave runup estimations are of great interest for coastal risk assessment and engineering design. Phase-resolving depth-integrated numerical models offer a promising alternative to commonly used empirical formulae at relatively low computational cost. Several operational models are currently freely available and have been extensively used in recent years for the computation of nearshore wave transformations and runup. However, recommendations for best practices on how to correctly utilize these models in computations of runup processes are still sparse. In this work, the Boussinesq-type model BOSZ is applied to calculate runup from irregular waves on intermediate and reflective beaches. The results are compared to an extensive laboratory data set of LiDAR measurements from wave transformation and shoreline elevation oscillations. The physical processes within the surf and swash zones such as the transfer from gravity to infragravity energy and dissipation are accurately accounted for. In addition, time series of the shoreline oscillations are well captured by the model. Comparisons of statistical values such as R2% show relative errors of less than 6%. The sensitivity of the results to various model parameters is investigated to allow for recommendations of best practices for modeling runup with phase-resolving depth-integrated models. While the breaking index is not found to be a key parameter for the examined cases, the grid size and the threshold depth, at which the runup is computed, are found to have significant influence on the results. The use of a time series, which includes both amplitude and phase information, is required for an accurate modeling of swash processes, as shown by computations with different sets of random waves, displaying a high variability and decreasing the agreement between the experiment and the model results substantially. The infragravity swash SIG is found to be sensitive to the initial phase distribution, likely because it is related to the short wave envelope.