Vibrational dynamics of adsorbates near surfaces plays both an important role for applied surface science and as a model lab for studying fundamental problems of open quantum systems. We employ a previously developed model for the relaxation of a D–Si–Si bending mode at a D:Si(100)-(2 × 1) surface, induced by a “bath” of more than 2000 phonon modes [Lorenz and P. Saalfrank, Chem. Phys. 482, 69 (2017)], to extend previous work along various directions. First, we use a Hierarchical Effective Mode (HEM) model [Fischer et al., J. Chem. Phys. 153, 064704 (2020)] to study relaxation of higher excited vibrational states than hitherto done by solving a high-dimensional system–bath time-dependent Schrödinger equation (TDSE). In the HEM approach, (many) real bath modes are replaced by (much less) effective bath modes. Accordingly, we are able to examine scaling laws for vibrational relaxation lifetimes for a realistic surface science problem. Second, we compare the performance of the multilayer multiconfigurational time-dependent Hartree (ML-MCTDH) approach with that of the recently developed coherent-state-based multi-Davydov-D2 Ansatz [Zhou et al., J. Chem. Phys. 143, 014113 (2015)]. Both approaches work well, with some computational advantages for the latter in the presented context. Third, we apply open-system density matrix theory in comparison with basically “exact” solutions of the multi-mode TDSEs. Specifically, we use an open-system Liouville–von Neumann (LvN) equation treating vibration-phonon coupling as Markovian dissipation in Lindblad form to quantify effects beyond the Born–Markov approximation.