The muonic component of extensive air showers (EAS) contains information about the mass
composition of cosmic rays and hadronic interactions occurring early in the EAS development.
While propagating through the atmosphere, muons lose a small fraction of their energy and
deviate little from a straight-line trajectory, retaining information about their production points.
It is, therefore, possible to reconstruct the Muon Production Depth (MPD) from the arrival time
and position of muons measured by ground arrays of cosmic-ray observatories. The estimation
of the kinematic delay, a quantity directly related to the energy of muons, is key for the MPD
reconstruction. The existing kinematic delay parametrization, tailored to showers with zenith
angles θ ~ 60° and muons arriving far from the shower core, represents the dominant contribution
to the method’s systematic uncertainties. In this contribution, we extend the MPD reconstruction
to showers with θ < 60° while considerably reducing the present radial cut and applying it to
energies between the second knee and ankle of the cosmic-ray spectrum, where overlap with the
nominal energy at the LHC exists. To filter out the electromagnetic component of EAS, we select
muons energetic enough to reach underground detectors. We use Gradient-Boosted Decision
Trees, a machine learning algorithm suited for structured, heterogeneous datasets, to reconstruct
the MPD. We report an unbiased MPD reconstruction (< 10 g cm$^{-2}$), with a muon-by-muon
resolution of ~ 80 (~ 170) g cm$^{-2}$ for θ = 0° (60°) showers. The method’s applicability to higher
energies and different primary species is also investigated.