A model-operator approach to fully relativistic calculations of the nuclear recoil effect on energy levels in many-electron atomic systems is worked out. The one-electron part of the model operator for treating the normal mass shift beyond the Breit approximation is represented by a sum of semilocal and nonlocal potentials. The latter ones are constructed by employing the diagonal and off-diagonal matrix elements rigorously evaluated for hydrogenlike ions to first order in the electron-to-nucleus mass ratio. The specific mass shift beyond the lowest-order relativistic approximation has a form which can be directly employed in calculations. The capabilities of the method are probed by comparison of its predictions with the results of ab initio QED calculations. The proposed operator can be easily incorporated into any relativistic calculation based on the Dirac-Coulomb-Breit Hamiltonian.