Abstract ALICE is a large experiment at the CERN Large Hadron Collider. Located 52 meters underground, its detectors are suitable to measure muons produced by cosmic-ray interactions in the atmosphere. In this paper, the studies of the cosmic muons registered by ALICE during Run 2 (2015–2018) are described. The analysis is limited to multimuon events defined as events with more than four detected muons (Nμ > 4) and in the zenith angle range 0° < θ < 50°. The results are compared with Monte Carlo simulations using three of the main hadronic interaction models describing the air shower development in the atmosphere: QGSJET-II-04, EPOS-LHC, and SIBYLL 2.3d. The interval of the primary cosmic-ray energy involved in the measured muon multiplicity distribution is about 4 × 1015 < E prim < 6 × 1016 eV. In this interval none of the three models is able to describe precisely the trend of the composition of cosmic rays as the energy increases. However, QGSJET-II-04 is found to be the only model capable of reproducing reasonably well the muon multiplicity distribution, assuming a heavy composition of the primary cosmic rays over the whole energy range, while SIBYLL 2.3d and EPOS-LHC underpredict the number of muons in a large interval of multiplicity by more than 20% and 30%, respectively. The rate of high muon multiplicity events (Nμ > 100) obtained with QGSJET-II-04 and SIBYLL 2.3d is compatible with the data, while EPOS-LHC produces a significantly lower rate (55% of the measured rate). For both QGSJET-II-04 and SIBYLL 2.3d, the rate is close to the data when the composition is assumed to be dominated by heavy elements, an outcome compatible with the average energy E prim ∼ 1017 eV of these events. This result places significant constraints on more exotic production mechanisms.