Earthquake ruptures are shear fractures that propagate along pre-existing zones of weakness in the earth's crust (faults). Models couple seismic wave propagation (and stress transfer) through rocks surrounding the fault with highly nonlinear interface jump conditions relating fault strength to slip velocity (the discontinuity in the particle velocity field across the fault), temperature, and pressure of pore fluids within the core of the fault zone. The elastic wave equation is discretized using summation-by-parts (SBP) finite differences, and the nonlinear jump conditions are enforced weakly using the simultaneous approximation term (SAT) method. This approach is provably stable and accurate, but the semi-discretized equations can become prohibitively stiff unless the jump conditions are formulated in a particular way. Coordinate transformation techniques allow us to model ruptures on nonplanar fault surfaces. A combination of numerical models and asymptotic analyses using boundary perturbation techniques are used to investigate the role of fault roughness (observations show natural faults are self-similar fractal surfaces) in determining the resistance to slip and the production of incoherent high frequency seismic waves that are ubiquitous is seismograms.