A novel quantum dynamical method to simulate vibronic reaction dynamics in molecules at metal surfaces is proposed. The method is based on the hierarchical quantum master equation approach and uses a discrete variable representation of the nuclear degrees of freedom in combination with complex absorbing potentials and an auxiliary source term. It provides numerically exact results for a range of models. By taking the coupling to the continuum of electronic states of the surface properly into account, nonadiabatic processes can be described and the effect of electronic friction is included in a nonperturbative and non-Markovian way. Illustrative applications to models for desorption of a molecule at a surface and the current-induced bond rupture in single-molecule junctions demonstrate the performance and versatility of the method.