Modeling the dynamics of a quantum system coupled to a dissipative environment becomes particularly challenging when the system's dimensionality is too high to permit the computation of its eigenstates. This problem is addressed by introducing an eigenstate-free formalism, where the open quantum system is represented as a mixture of high-dimensional, time-dependent wave packets governed by coupled Schrödinger equations, while the environment is described by a multi-component quantum master equation. An efficient computational implementation of this formalism is presented, employing a variational mixed Gaussian/multiconfigurational time-dependent Hartree (G-MCTDH) ansatz for the wave packets and propagating the environment dynamics via hierarchical equations, truncated at the first or second level of the hierarchy. The effectiveness of the proposed methodology is demonstrated on a 61-dimensional model of phonon-driven vibrational relaxation of an adsorbate. G-MCTDH calculations on 4- and 10-dimensional reduced models, combined with truncated hierarchical equations for the mean fields, nearly quantitatively replicate the full-dimensional quantum dynamical results on vibrational relaxation while significantly reducing the computational time. This approach thus offers a promising quantum dynamical method for modeling complex system-bath interactions, where a large number of degrees of freedom must be explicitly considered.