Finite-fault source inversions reveal the spatial complexity of earthquake slip over the fault plane. We develop a stochastic characterization of earthquake slip complexity, based on published finite-source rupture models, in which we model the distribution of slip as a spatial random field. The model most consistent with the data follows a von Karman autocorrelation function (ACF) for which the correlation lengths a increase with source dimension. For earthquakes with large fault aspect ratios, we observe substantial differences of the correlation length in the along-strike (ax) and downdip (az) directions. Increasing correlation length with increasing magnitude can be understood using concepts of dynamic rupture propagation. The power spectrum of the slip distribution can also be well described with a power law decay (i.e., a fractal distribution) in which the fractal dimension D remains scale invariant, with a median value D = 2.29 ¿ 0.23, while the corner wave number kc, which is inversely proportional to source size, decreases with earthquake magnitude, accounting for larger slip patches for large-magnitude events. Our stochastic slip model can be used to generate realizations of scenario earthquakes for near-source ground motion simulations.