Seismic and geodetic observations show that slow slip events (SSEs) in subduction zones can happen at all temporal and spatial scales and propagate at various velocities. Observation of rapid tremor reversals (RTRs) indicates back-propagating fronts traveling much faster than the main rupture front. Heterogeneity of fault properties, such as fault roughness, is a ubiquitous feature often invoked to explain this complex behavior, but how roughness affects SSEs is poorly understood. Here we use quasi-dynamic seismic cycle simulations to model SSEs on a rough fault, using normal stress perturbations as a proxy for roughness and assuming rate-and-state friction, with strengthening behavior at high slip rate. SSEs exhibit temporal clustering, large variations in rupture length and propagation speed, and back-propagating fronts at different scales. We identify a mechanism for back propagation: as ruptures propagate through low-normal stress regions, a rapid increase in slip velocity combined with rate-strengthening friction induces stress oscillations at the rupture tip, and the subsequent ”delayed stress drop’ induces secondary back-propagating fronts. Moreover, on rough faults with fractal elevation profiles, the transition from pulse to crack can also lead to the re-rupture of SSEs due to local variations in the level of heterogeneity. Our study provides a possible mechanism for the complex evolution of SSEs inferred from geophysical observations and its link to fault roughness.