Computing the hazard posed by the next large earthquake on the San Andreas fault is best carried out with a simulation-based approach. Here we discuss a numerical simulation, Virtual California, that includes many of the physical processes known to be important in earthquake dynamics. These include elastic interactions among the faults in the model, driving at the correct plate tectonic rates, and frictional physics on the faults using the physics obtained from laboratory models with parameters consistent with the occurrence of historic earthquakes. We report progress on a variety of problems relating to the construction and use of increasingly realistic models for earthquakes on the San Andreas fault system, which will be required as the Working Group on California Earthquake Probabilities moves into its next phase of hazard and risk analysis. One of the important issues is to construct a fault system model based upon current WGCEP data. Here we use Deformation Model 2.2 ( http://www.relm.org/models/WGCEP/ ) to produce the most realistic model to date. We also have more fine-scale versions of previous models, including one model having 3 x 2**12 fault elements. While previous versions of Virtual California used only vertical strike slip faults, we are now incorporating dipping rectangular faults having arbitrary rake anble into the model as well. Versions of the basic code are available in Fortran, C, and object-oriented C++. One of the issues that we have encountered is the existence of a dynamical instability that arises as a direct result of the basic interactions between the fault elements, combined with the requirement that the long term slip rate on all fault elements match the observed field- derived average. Using these new models and simulations, we are engaged in novel types of data assimilation, using a method of "scoring" the simulation in comparison to observed paleoseismic data. In this paper, we summarize these results and discuss the implications for numerical forecasting methods.