We demonstrate an approach to simulate rock glacier dynamics based on a three-dimensional stick-slip model coupled to a simple model of melt water dispersion. Simulation results account for observed kinematic discontinuities and exhibit a distribution of slip event sizes that is characterized by large fluctuations following a power law. The model presents a feedback loop of dynamics and internal structure: (1) the evolution of shear horizons, which are correlated to melt water occurrence, can be seen; and (2) melt water distribution narrows to regions of maximum slip. Although this approach is promising, improvements in some details must be made. The solution for stress transfer accompanying slip has been taken from a model of a sheared body that is vertically constrained. In future models, the boundary conditions must be modified to that of a body confined only at its base. Additionally, a heuristic treatment of hydrology currently used must be further quantified.