In this work, the model of rate-and-state friction, which can be viewed as central to the numerical simulation of earthquakes, is considered from a mathematical point of view. First, a framework is presented through which a general class of such friction laws can be understood and analysed. A prototypical viscoelastic problem of earthquake rupture is then formulated, both in strong and in variational form. Analysis of this problem is difficult, since the incorporation of rate-and-state friction leads to a coupling of variables. In a time-discrete setting, nonetheless, results on existence, uniqueness, and continuous parameter dependence of solutions can be obtained. The principal idea is to reformulate the variable interdependence as a fixed point problem and to prove convergence for a corresponding iteration. With that in mind, next, a numerical algorithm is proposed that resolves the coupling through a fixed point iteration. Since it puts a state-of-the-art solver and adaptive time stepping to use, it is not only stable but also fast. Its applicability to problems of interest is demonstrated in the penultimate chapter, which focuses on simulations of megathrust earthquakes that form at the base of a subduction zone. The main assumptions made throughout this work are summarised and discussed in the last chapter.