In order to better understand the interdependence between tectonic loading, earthquake nucleation and after-slip phenomena, geophysicists require accurate and performant simulations of sequences of earthquakes and aseismic slip (SEAS) that can cover time-scales of several centuries. The model is based on the rate-and-state friction law, which considers ageing at the tectonic fault and stems from experiments. The problem is formulated as an elliptic partial differential equation with algebraic constraints. This thesis investigates the accuracy and performance of time-adaptive explicit Runge-Kutta and implicit backward differentiation formula schemes. First, the quality of the error estimates is analyzed on a simple one-dimensional domain where an analytic solution is available. Next, a two-dimensional domain with antiplane shear motion is considered on which four different formulations of the constrained differential equation are derived. All formulations
can be solved with implicit and some additionally with explicit time integration methods. It is found that the dimension of the problem can be reduced by one and that the algebraic condition can be removed by transforming the problem into a second order ODE. The simulation can therefore be accelerated by a factor of ten. Finally, the formulations are extended to a general three-dimensional problem and similar accuracy and performance are observed.