SUMMARYOur understanding of the seismicity of continental interiors, far from plate margins, relies on the ability to account for behaviours across a broad range of time and spatial scales. Deformation rates around seismic faults range from the slip-on-fault during earthquakes to the long-term viscous deformation of surrounding lithosphere, thereby presenting a challenge to modelling techniques. The aim of this study was to test a new method to simulate seismic faults using a continuum approach, reconciling the deformation of viscoelastoplastic lithospheres over geological timescales. A von Mises yield criterion is adopted as a proxy for the frictional shear strength of a fault. In the elastoplastic fault models a rapid change in strength occurs after plastic yielding, to achieve stress–strain equilibrium, when the coseismic slip and slip velocity from the strain-rate response and size of the fault are calculated. The cumulative step-function shape of the slip and temporally partitioned slip velocity of the fault demonstrated self-consistent discrete fault motion. The implementation of elastoplastic faults successfully reproduced the conceptual models of seismic recurrence, that is strictly periodic and time- and slip-predictable. Elastoplastic faults that include a slip velocity strengthening and weakening with reduction of the time-step size during the slip stage generated yield patterns of coseismic stress changes in surrounding areas, which were similar to those calculated from actual earthquakes. A test of fault interaction captured the migration of stress between two faults under different spatial arrangements, reproducing realistic behaviours across time and spatial scales of faults in continental interiors.