The deep roots of subduction megathrusts exhibit aseismic slow slip events, commonly accompanied by tectonic tremor. Observations from exhumed rocks suggest this region of the subduction interface is a shear zone with frictional lenses embedded in a viscous matrix. Here we use numerical models to explore the transient slip characteristics of finite-width frictional-viscous megathrust shear zones. Our model utilizes an invariant, continuum-based, regularized form of rate- and state-dependent friction (RSF) and simulates earthquakes along spontaneously evolving faults embedded in a 2D heterogeneous continuum. The setup includes two elastic plates bounding a viscoelastoplastic shear zone (subduction interface melange) with inclusions (clasts) of varying distributions and viscosity contrasts with respect to the surrounding weaker matrix. The entire shear zone exhibits the same velocity-weakening RSF parameters, but the lower viscosity matrix has the capacity to switch between RSF and viscous creep as a function of local stress state. Results show that for a range of matrix viscosities near the frictional-viscous transition, viscous damping and stress heterogeneity in these shear zones both 1) sets the ‘speed limit’ for earthquake ruptures that nucleate in clasts such that they propagate at slow velocities; and 2) permits the transmission of slow slip from clast to clast, allowing slow ruptures to propagate substantial distances over the model domain. For reasonable input parameters, modeled events have moment-duration statistics, stress drops, and rupture propagation rates that match natural slow slip events. These results provide new insights into how geologic observations from ancient analogs of the slow slip source may scale up to match geophysical constraints on modern slow slip phenomena.