We present a split-node finite difference (FD) method for modelling shear ruptures that is consistently fourth-order accurate in its spatial discretization, both in the interior of the model and at the fault. The method, called mimetic operator split node (MOSN), uses a staggered grid, and the fault plane is discretized using split nodes for both discontinuous displacement and discontinuous stress components. The method differs in several ways from previous FD methods for rupture modelling. (1) It uses one-sided differentiation (in the fault-normal direction) to retain high-order accuracy at the discontinuity. (2) The one-sided differentiation is implemented with mimetic FD operators that conserve the integration-by-parts formula (including boundary terms) in discrete form. (3) The mimetic operators lead naturally to the evaluation of conjugate tractions at the split-node displacement sites, as required to apply the frictional-sliding boundary conditions. In the interior of the grid, the mimetic operators reduce to the conventional fourth-order staggered-grid FD operators. We verify the accuracy of the MOSN method using two test problems: a fixed-speed rupture case (Kostrov's problem), and the case of a spontaneous rupture with slip-dependent friction. In the former test, comparisons of slip and shear stress waveforms confirm the convergence of numerical results to the analytical solution as the grid is refined. In the latter test, assessment is based on both qualitative waveform characteristics and a quantitative convergence analysis that uses error metrics in rupture time, final slip and peak slip rate. Compared with a previously verified second-order scheme, the MOSN method shows significant reduction in artefacts attributable to numerical dispersion and also yields a considerable improvement in rupture arrival times. While from a practical standpoint the performance gains of MOSN over the second-order scheme are considerable, MOSN has comparatively poorer convergence rates, and loses its advantage in accuracy in the limit of very small spatial steps. We attribute the relatively low convergence rates to our use of a friction law that induces gradient discontinuities in the fault traction, violating the smoothness assumptions underlying the accuracy estimates, a limitation which may not apply to more complex, experimentally derived friction laws with smooth dependence on slip-rate and state variables. Our current implementation is 2-D and assumes a straight fault, but the method is extensible to more general 3-D geometries. © 2007 The Authors Journal compilation © 2007 RAS.