Numerical simulations of slip evolution along a cellular vertical strike-slip fault in an elastic half-space are performed for several models representing discrete fault systems embedded in three-dimensional elastic continua. The geometry and imposed boundary conditions correspond approximately to the central San Andreas fault. The simulations incorporate brittle and creep deformations in series; the net fault zone deformation rate is the sum of creep rate and frictional slip rate. Brittle fault properties are given by various distributions of earthquake stress drops on failing segments (numerical cells). The assumed distributions represent two idealized situations corresponding to different extreme states along an evolutionary path of a fault: (1) a strongly disordered state characterized by a wide range of size scales, representing immature fault zones and extended spatial domains, and (2) a relatively regular state having a narrow range of size scales, representing mature highly-slipped faults. The assumed creep properties are identical in all cases. These are prescribed in terms of coefficients characterizing a power law dependency of creep-slip rate on stress. The combined brittle-creep process and employed parameters lead to an overall ''pine-tree'' stress-depth profile with a ''brittle-ductile'' transition depth of about 12.5 km, and variable stress-along-strike profiles with ''brittle-creep'' transition around 65 km NW of the 1857 rupture. The spatial patterns of simulated hypocenters are statistically similar to observed data. The results indicate that the range of size scales characterizing strong fault zone heterogeneities has important manifestations on the seismic response of a fault system. A narrow range of size scales leads to frequency-size statistics of earthquakes resembling the characteristic earthquake distribution, and quasi-periodic temporal distribution of large events as in the seismic gap hypothesis. On the other hand, a wide range of size scales leads to Gutenberg-Richter power law frequency-size statistics, and random or clustered temporal distribution of large events. The simulations demonstrate that treatment of the various observed forms of frequency-size and temporal statistics of earthquakes can be unified through the concept of range of size scales characterizing fault zone heterogeneities. This has a clear physical interpretation in terms of structural properties of a given fault zone or broad lithospheric domain, and is supported by observed earthquake and fault data. In some simulated cases the frequency-size statistics of small earthquakes fall sharply below the self-similar Gutenberg-Richter line. The results indicate that small earthquakes prepare the fault for the occurrence of a large event by smoothing, during gradual tectonic loading, the long-wavelength components of stress on the fault. This is donethrough short-wavelength stress roughening associated with the numerous ruptures of the small events. The above pattern of smoothing/roughening of long/short wavelengths of stress on the fault is reversed during large-scale ruptures of the big events. ¿ American Geophysical Union 1996