John B. Rundle, University of Colorado, Boulder, CO,Lead Investigator

Geoffrey Fox, NPAC, Syracuse University, Syracuse, NY, Co-Lead Investigator

William Klein, Boston University, Boston, MA

John Reynders, ACL, Los Alamos National Laboratory, Los Alamos, NM

Earth Scientists:

Thomas H. Jordan, Massachusetts Institute of Technology, Cambridge, MA

Hiroo Kanamori, California Institute of Technology, Pasadena, CA

Karen Carter Krogh, EES, Los Alamos National Laboratory, NM

Christopher J. Marone, Massachusetts Institute of Technology, Cambridge, MA

J. Bernard Minster, University of California, San Diego, CA

John B. Rundle, University of Colorado, Boulder, CO

Ross Stein, United States Geological Survey, Menlo Park, CA

Donald Turcotte, Cornell University, Ithaca, NY

Condensed Matter Scientists:

William Klein, Boston University, Boston, MA

James Crutchfield, Santa Fe Institute, Santa Fe, NM

Computer/Computational Scientists:

Geoffrey Fox, NPAC, Syracuse University, Syracuse, NY

Julian Cummings, ACL, Los Alamos National Laboratory, Los Alamos, NM

Roscoe Giles, Boston University, Boston, MA

Steven Karmesin, ACL, Los Alamos National Laboratory, Los Alamos, NM

Paul Messina, CACR, California Institute of Technology, Pasadena, CA

John Salmon, CACR, California Institute of Technology, Pasadena, CA

Bryant York, Northeastern University, Boston, MA

John Reynders, ACL, Los Alamos National Laboratory, Los Alamos, NM

Michael Warren, ACL, Los Alamos National Laboratory, Los Alamos, NM

Tim Williams, ACL, Los Alamos National Laboratory, Los Alamos, NM

Unfunded Collaborators:

Claude Allegre, Institute de Physique du Globe, Paris, France

Shawn Larsen, Lawrence Livermore National Laboratory, Livermore, CA

Lawrence Hutchings, Lawrence Livermore National Laboratory, Livermore, CA

William Foxall, Lawrence Livermore National Laboratory, Livermore, CA

Charles Carrigan, Lawrence Livermore National Laboratory, Livermore, CA

CIRES CB 216

University of Colorado

Boulder, CO 80309

Department of Physics

Syracuse Univresity

Syracuse, NY 13244

William Klein

Department of Physics

Boston University

Boston, MA 02215

John Reynders

Tim Williams

Julian Cummings

Steven Karmesin

Karen Carter Krogh

Advanced Computing Lab

Los Alamos National Lab

Los Alamos, NM 87545

Santa Fe Institute

1399 Hyde Park Road

Santa Fe, NM 87501

Thomas Jordan

Department of EAPS, MIT

Cambridge, MA 02139

Seismological Lab

California Inst. of Technology

Pasadena, CA 91125

IGPP - SIO

University of California

La Jolla, CA 92093

Department of Physics, CACR

California Inst. of Technology

Pasadena, CA 91125

Earthquake Hazards Team

United States Geological Survey

Menlo Park, CA 94025

Paul Messina

CACR

California Inst. of Technology

Pasadena, CA 91125

Dept. of Geology

Cornell Univ., Ithaca, NY 14853

College of Computer Science

Northeastern University

Boston, MA 02115

IPG Paris, 4 Place Jussieu

Tour 14-24, 3e etage

75252 Paris Cedex 05

France

Shawn Larsen

Lawrence Hutchings

William Foxall

Lawrence Livermore Nat. Lab.

Livermore, CA 94550

Our primary objective is to develop the computational capability to carry out large scale numerical simulations of the physics of earthquakes in southern California and elsewhere. These simulations will produce temporal and spatial patterns of earthquakes, surface deformation and gravity change, seismicity, stress, as well as, in principle, other variables, including pore fluid and thermal changes, for comparison to field and laboratory data. To construct of the simulations, we will develop a state-of-the-art problem solving environment that will facilitate: 1) Construction of numerical and computational algorithms and specific environment(s) needed to carry out large scale simulations of these nonlinear physical processes over a geographically distributed, heterogeneous computing network; and 2) Development of a testbed for earthquake "forecasting" & "prediction" methodologies which uses modern Object Oriented techniques and scalable systems, software and algorithms which are efficient for both people and computational execution time. The GEM approach will allow the physics of large networks of earthquake faults to be analyzed within a general theoretical framework for the first time. It will stimulate changes in earthquake science in much the same way that General Circulation Models have changed the way atmospheric science is carried out. The computational techniques developed by the project will have important applications in many other large, computationally hard problems, such as 1) Large particle-in-cell systems in ASCI; 2) Statistical physics approaches to random field spin systems; and 3) Simulating large neural networks with learning and cognition.

A. Proposal Abstract 3

B. Table of Contents 4

1 One Page Executive Summary for Reviewers 5

2 The GEM Grand Challenge 6

3 Introduction to the GEM Approach 11

4 Recent Observations of Earthquakes with Significance for GEM 13

5 Scientific Framework for GEM 14

6 Data Issues for Validation of Simulations 19

7 Overview of Computational Science Issues 26

8 Current Status of Earthquake Simulation Research 27

9 Project Plan for the GEM Grand Challenge 32

10 Individual Work Statements 37

the Southern California Earthquake Center 45

E. Bibliography 46

F. Computational Facilities 57

G. Project Personnel 59

H. Current and Pending Support 61

I. Budget 68

Objectives: The primary objective is to develop the computational capability to carry out large scale numerical simulations of the physics of earthquakes in southern California and elsewhere. To meet this objective, we will develop a state of the art problem solving environment that will facilitate: 1) The construction of numerical and computational algorithms and specific environment(s) needed to carry out large scale simulations of these nonlinear physical processes over a geographically widely distributed, heterogeneous computing network; and 2) The development of a testbed for earthquake "forecasting" & "prediction" methodologies which uses modern Object Oriented techniques and scalable systems, software and algorithms which are efficient for both people and computational execution time.

Method: We will base our work on numerical simulation techniques initially developed by Stuart (1986), Rundle (1988) and Ward and Goes (1993), using these efforts as starting points to model the physics of earthquake fault systems in southern California. The problem solving environment will be built using the best parallel algorithms and software systems. It will leverage DOE and other state of the art national activities in simulation of both cellular automata and large scale particle systems.

Scientific and Computational Foci: We will focus on developing the capability to carry out large scale simulations of complex, multiple, interacting fault systems using a software environment optimized for rapid prototyping of new phenomenological models. The software environment will require: 1) Developing algorithms for solving computationally difficult nonlinear problems involving ("discontinuous") thresholds and nucleation events in a networked parallel (super) computing environment, adapting new "fast multipole" methods previously developed for general N-body problems; 2) Leveraging the Los Alamos ACL Infrastructure, including the POOMA object oriented framework, which is already being applied to computationally related adaptive particle mesh problems; and 3) Developing a modern problem solving environment to allow researchers to rapidly integrate simulation data with field and laboratory data (visually and quantitatively). This task is becoming increasingly important with the exponential increase in seismic, as well as space-based deformation data.

Significance of Anticipated Results: The GEM approach will allow the physics of large networks of earthquake faults to be analyzed within a general theoretical framework for the first time. It will stimulate changes in earthquake science in much the same way that General Circulation Models have changed the way atmospheric science is carried out. The computational techniques developed by the project will have important applications in many other large, computationally hard problems, such as 1) Large particle-in-cell systems in ASCI; 2) Statistical physics approaches to random field spin systems; as well as 3) Simulating large neural networks with learning and cognition.

Investigator Team and Project Organization: We have assembled a group of internationally recognized experts in the three areas of 1) observational and theoretical earth science 2) statistical mechanics and complex systems and 3) computer and computational science. The latter include world experts in the critical algorithms, software and systems required. We are represented by key people both inside DOE and in universities in all areas, so technology transfer to related projects, as well as educational benefits, will follow easily. Rundle will serve as Principal Investigator. The Investigators will participate in periodic workshops at which 1) results will discussed; and 2) specific research avenues will be formulated on a regular and timely basis.

Status: Many basic facts about earthquakes have been known since ancient times (see e.g., Richter, 1958; Scholz, 1990; and the review by Rundle and Klein, 1995), but only during the last hundred years have earthquakes been known to represent the release of stored elastic strain energy along a fault (Reid, 1908). Although a great deal of data has accumulated about the phenomenology of earthquakes in recent years, these events remain one of the most destructive and poorly understood of the forces of nature. The importance of studying earthquakes, as well as the development of techniques to eventually predict or forecast them, has been underscored by the fact that an increasing proportion of global population lives along active fault zones (Bilham, 1996). Recent research indicates that earthquakes exhibit a wealth of complex phenomena, including space-time clustering of events, scaling, as represented by the Gutenberg-Richter relation (e.g., Turcotte, 1992), and space-time migration of activity along fault systems (e.g., Scholz, 1990). Analyzing and understanding these events is made particularly difficult due to the long time scales, usually of the order of tens to hundreds of years between earthquakes near a given location on a fault, and by the space-time complexity of the historical record.

During the last two decades, a series of national policy decisions and programs culminated, five years ago, in the establishment of the Southern California Earthquake Center (html://www.usc.edu/dept/earth/quake/). Togther with efforts undertaken several decades ago by the United States Geological Survey (html://www.usgs.gov/themes/earthqk.html), the accuracy, reliability, and availability of observational data for earthquakes, particularly in southern California, have increased enormously. Moreover, our knowledge about earthquake related processes and the earthquake source itself is at a level far beyond what we knew two decades ago, prior to the National Earrthquake Hazard Reduction Program. Finally, it is likely that the data available will increase dramatically when an operational Synthetic Aperature Radar Interferometry mission begins routine operation (see, e.g., Massonet et al., 1993).

Despite this, we still find ourselves, as a scientific community, unable to even approximately predict the time, date, and magnitude of earthquakes. At the moment, the best that can be done is embodied in the Phase II report of earthquake probabilities published by the Southern California Earthquake Center (SCEC, 1995; please refer to SCEC web page). These probabilities are based on "scenario earthquakes" and crude assumptions about whether, for example, contiguous segments of faults ("characteristic earthquakes") do or do not tend to rupture together to produce much larger events.

Earthquakes, even those significantly smaller than the largest magnitude events of about magnitude 9.5, are capable of producing enormous damage. For example, the magnitude ~ 7 January 17, 1994 Kobe, Japan earthquake was responsible for an estimated $200 Billion in damages, accounting for lost economic income as well as direct damage. This event was a complete surprise, inasmuch as the immediate region had been relatively inactive in historic time (see, e.g., Trans. AGU, v 76 Supplement, 1995).

Earthquake science has by its nature been observational in character. The inability to carry out repetitive experiments on fault systems means that it is extremely difficult to test hypotheses in a systematic way, via controlled experiments. In addition, time intervals between major earthquakes at a given location are typically much longer than human life spans, forcing reliance on unreliable historical records. This state of affiars is in many respects similar to that of the atmospheric science community in the late 1940's and early 1950's. At that time, John von Neumann and Jule Charney decided to initiate a program to develop what were to become the first General Circulation Models (GCM's) of the atmosphere (see, e.g., Weatherwise, June/July 1995 for a retrospective). These modeling and simulation efforts revolutionized research within the atmospheric science community, because for the first time a feedback loop was established between observations and theory. Hypotheses could be formulated on the basis of computer simulations that could be tested by field observations. Moreover questions arising from new observations could be answered with computer "experiments".

There is thus a growing consensus in the scientific community that the time has come to establish a similar feedback loop between observations, theory and computer simulations within the field of earthquake science. The goals of the General Earthquake Model (GEM) project are similar to those of the GCM community: 1) to develop sophisticated computer simulations based upon the best available physics, with the goal of understanding the physical processes that determine the evolution of active fault systems, and 2) to develop a Numerical Earthquake Forecast capability that will allow current data to be projected forward in time so that future events can at best be predicted, and at worst perhaps anticipated to some degree.

C.2.2 Computational Science Status

This grand challenge not only tackles a problem of great interest to society and the earth science section of DoE but also will make use of the talent and novel infracture of the computer science, algorithm and physics strengths of Los Alamos. The team must create a computational environment to understand the properties and predict the behavior of critical seismic activity - earthquakes. The complexity of the seismic grand challenge will entail the development of novel algorithms based upon innovative computational architectures with new software paradigms. The development of a predictive capability will require an enormous computational resource - one which the Los Alamos Advanced Computing Laboratory (ACL) will be able to provide through a clustered Shared Memory Processor (SMP) approach to a TeraFlop. The computational resource alone, however, is useless without the system substrates necessary to integrate the clustered SMP system into a coherent platform and a layered framework of class libraries which separate the essential abstractions of the problem domain. There are several levels of complexities in the seismic grand challenge including parallelism, seismic physics, object-orientation, adaptive-mesh-refinement, and multipole methods. The many levels of complexites present in the current and future generations of Grand Challenge simulations will demand an interactive team of physicists and computer scientists working together to attack the problem.

Over the next three years, the team will create an object-oriented Seismic FrameWork of layered class libraries for high-performance predictive seismic simulation. This effort will build upon the Parallel Object-Oriented Methods and Applications (POOMA) FrameWork at the Los Alamos Advanced Computing Laboratory. We have already suggested why POOMA is particularly attractive. The multipole approach we describe in a later section, which has been widely applied to particle simulation, is not well suited to languages like Fortran (and the parallel HPF) because it uses complex hierarchical data structures provided by POOMA's C++ base. Further, the need for rapid prototyping of a wide variety of physics modules in GEM calls out for POOMA's object oriented capabilities. We believe POOMA is recognized as the most sophisticated system of its type in the HPCC community and Los Alamos has unique capabilities that are needed for the support of GEM.

FrameWork description: The Seismic FrameWork will provide an integrated layered system of objects in which each object higher in the FrameWork is composed of or utilizes objects which are lower in the FrameWork. In the case of the POOMA FrameWork (upon which the Seismic FrameWork will be built), the higher layers provide objects which capture the main abstractions of scientific problem domains (particles, fields, matrices) in a representation of mathematical expressions preserved in the application source code. The objects which are lower in the FrameWork focus more on capturing the abstractions relevant to parallelism and efficient node-level simulation such as communications, domain decomposition, threading, and load balancing.

This layered approach provides a natural breakdown of responsibility in application development. Computer science specialists can focus on the lower realms of the FrameWork optimizing kernels and message-passing techniques without having to know the details of the application being constructed. The physicists, on the other hand, can construct algorithms, components, and simulations with objects in the higher levels of the FrameWork without having to know the specific implementation details of the objects which compose the application. This design allows for efficient team code development with no single member of the team having to know all the details of the computer and physics problem domains. Indeed, with the Seismic FrameWork it will be possible for a scientist to build an application which runs on a variety of distributed memory architectures without any understanding of the hierarchical memory and communication present in the ACL clustered SMP system.

Leverage of other DoE Programs: The Seismic FrameWork will leverage with current missions driving the POOMA FrameWork. Most notably, the parallel dipole simulation on AMR geometries will benefit directly from current efforts in the ASCI program to build a particle-based transport simulation with AMR meshes upon the POOMA FrameWork. This effort in turn is built upon particle based abstractions within the POOMA FrameWork originally from the ACL Object-Oriented Particles Simulation (OOPS) class library - a FrameWork for gyrokinetic plasma simulation built for the Numerical Tokamak Grand Challenge. Indeed, most of the key abstractions have already been developed granting both the Seismic Grand Challenge effort and the ASCI transport effort an enormous benefit in code reuse and program leveraging.

Systems Integration: The software infrastructure goals of the Seismic Grand Challenge align directly with the infrastructure goals of the ACL. In particular, the object-based framework approach to building application, analysis, and visualization components are a superb fit for the Seismic FrameWork. To this end, our approach to the Seismic Grand Challenge will be object-oriented and modular, enabling many of our unique capabilities such as multipole-AMR to be encapsulated within reusable components. This will enable rapid prototyping of new predictive capabilities into simulations and the leveraging of techniques across application domains. We seek to develop and share simulation, visualization, and analysis components for the Seismic FrameWork which will both build and be built upon capabilities at the ACL and other grand challenges.

Leverage of Center for Nonlinear Systems, LANL: It should be noted that Hans Frauenfelder, Director of the Center for Nonlinear Science at LANL is an unfunded collaborator. We intend to use personnel at CNLS as a resource in our investigations, although the specific modes of interaction and details of topics will depend on the nature of the problems that arise during the course of the research.

The investigators in the GEM project have been selected because their background

and expertise enables them all to fill necessary functions in the initial development of the

GEM simulations. These functions and roles are organized in the following table.

Funded Investigators

==============

Name Institution Role*

===== ========== =====

Rundle Colorado Lead Earth Science -- Develop earthquake models,

stat. mech approaches, validation of

simulations (AL, PSE, AN, VA, SCEC)

Fox Syracuse Lead Computer Science -- Develop multipole

algorithms and integrate projects internally

and with external projects including HPCC and WWW communities (AL, PSE, AN))

Klein Boston Statistical mechanics analogies and methods:

Langevin equations for fault systems

dynamics, especially meanfield models (AL,

AN)

Crutchfield Santa Fe Inst Pattern recognition algorithms for earthquake

forecasting, phase space reconstruction of

attractors, predictability (AN)

Giles Boston Object oriented friction model algorithms, Cellular

Automata computations (AL, AN, PSE)

York Northeastern Cellular Automata, implementation of computational

approaches (AL, PSE)

Jordan MIT Validating models against "slow earthquake" data

(VA, SCEC)

Marone MIT Validating friction models against laboratory data

(VA)

Kanamori CalTech Validating models against broadband earthquake

source mechanism data (VA, SCEC)

Krogh LANL Validating fault geometries and interactions using

field data, including stress patterns (VA)

Minster UCSD Implementation/analysis of rate & state friction

laws (AN, SCEC))

Messina, Salmon

Caltech Parallel multipole algorithms, linkage of model validation with simulation ( AL, VA, PSE) Stein USGS Visualization requirements for simulations, validating

models against geodetic data (PSE, SCEC)

Turcotte Cornell Nature of driving stresses from mantle processes,

scaling symmetries, validating models with

data on seismicity patterns and SAR

interferometry data (AN, VA)

Reynders, Williams, White, Warren

ACL, LANL Multipole algorithms, elliptic solvers, adaptive

multigrid methods, optimal use of SMP

systems, load balancing (AL, PSE)

Unfunded Collaborators. They will provide advice on:

=====================================

Frauenfelder & Colleagues

CNLS, LANL Analysis of simulations for patterns, phase space reconstruction, analogies to other systems,

scaling symmetries (AN)

Allegre IPG Paris Statistical models for friction based upon

scale invariance (AN)

Larsen, Hutchings,Foxall, Carrigan

LLNL Green's functions, finite element models for

quasistatic and dynamic problems, faulting in

porous-elastic media, applications to San

Francisco Bay area, visualization

requirements (AL, VA)

*Roles: AL) Algorithms; PSE) Problem Solving Environment; AN) Analysis by statistical mechanics/statistical mechanics; VA) Validation; SCEC) Interaction with SCEC/CalTech and other earthquake data bases

C.2.4 GEM Team: Management Plan

We base our management plan on experience with 1) current Grand Challenges and the well known NSF center CRPC -- Center for Research in Parallel Computation -- where

Syracuse and Los Alamos are members; and 2) current NSF STC centers, most especially the Southern California Earthquake Center, for which Rundle chairs the Advisory Council. Rundle, the PI, will have full authority and responsibility for making decisions as to appropriate directions for the GEM Grand Challenge. In particular he will approve budgets and work plans by each contractor and subcontractor. These must be aligned with the general and specific team goals. The PI will be advised by an executive committee made up of a subset of the PI's representing the key subareas and institutions. This committee will meet approximately every 6 months in person and use the best available collaboration technologies for other discussions. Fox will have particular responsibility for overall integration of the computer science activities and Reynders for interaction with DOE and LANL. Other subareas include Statistical Physics Models (Klein), Cellular Automata (Giles), Data Analysis & Model Validation (Jordan), Parallel Algorithms (Salmon). The expectation is that the executive committee will operate on a consensus basis. Note that the goals of the Grand Challenge are both Scientific (simulation of Earth Science phenomena) and Computational (development of an object oriented Problem Solving Environment). The needs of both goals will be respected in all planning processes and contributions in both areas will be respected and viewed as key parts for the mission of the project.

The executive committee will be expanded to a full technical committee comprising at least all the funded and unfunded PI's. The technical committee will be responsible for developing the GEM plan which will be discussed in detail at least every 12 months at the major annual meeting that we intend to hold for scientists inside and outside this project. As well as this internal organization, we expect DOE will naturally set up an external review mechanism. However we suggest that a GEM external advisory committee consisting of leading Earth and Computer Scientists should be set up and that it will attend GEM briefings and advise the PI as to changes of direction and emphasis.

C.2.5 Summary of Project and Proposal

In the following sections, the GEM project is summarized and the proposed work is described in detail. We begin with an introduction to the GEM (section C.3) that describes the rationale for the simulations we propose, then evaluate some of the decisions that must be made in executing the computations. Because all such simulations must be continually validated with real world data, we then enumerate (section C.4) some of the recent observations that are suggestive of the kind of physics that must be included in the model. We next describe the scientific framework for the GEM approach, focusing in particular on the fundamental equations and the Green's functions that will be needed. We finally arrive at one of the most important computational results expected to flow from our work, development and application of a set of fast multipole algorithms to computationally represent the network of faults in the simulations. We then discuss the nature of friction on fault surfaces, and present a summary of the six most important ideas in this area, followed by mention of problems relating to scaling symmetries and renormalization. In section C.6, we dicuss issues relating to use of data to validate our simulations, including both existing and anticipated data types. We then present an overview of the important computational science issues (C.7) followed by a review of the current status of earthquake simulation research (C.8). Finally, we present the project plan for the GEM Grand Challenge (C.9), and the work plans of the various investigators (C.10).

Rationale for Simulations. Understanding the physics of earthquakes is complicated by the fact that the large events of greatest interest recurr at a given location along an active fault only on time scales of the order of hundreds of years (e.g., Richter, 1958; Kanamori, 1983; Pacheco et al., 1992). To acquire an adequate data base of similar large earthquakes therefore requires the use of historical records, which are known to possess considerable uncertainty. Moreover, instrumental coverage of even relatively recent events is often inconsistent, and network coverage and detection levels can change with time (Haberman, 1982). Understanding the details of the rupture process is further complicated by the spatial heterogenity of elastic properties, the vagaries of near field instrumental coverage, and other factors (see for example Heaton, 1990; Kanamori, 1993). We are therefore motivated to use methods that provide insight into details of the rupture process which are complementary to the usual observational techniques (e.g., Kanamori, 1993).

Given these constraints, numerical simulation is a practical means of obtaining experimental information about the behavior of such a nonlinear system over many cycles of activity (Binder, 1979; Binder, 1984; Mouritsen, 1984; Ma, 1985; Yeomans, 1992). In fact, other than paleoseismology, it is the only means of obtaining information about the dynamics of a given fault system over many earthquake cycles. And even paleoseismology has clear limitations, because information on slip is only available at the surface, and will not reflect any of the complexity known to characterize the depth dependence of slip.

If adequate theoretical understanding cannot be developed for the simulations, it is unlikely that the dynamics of real faults can be similarly understood. Numerical simulation has been used extensively to study earthquakes in the recent past (Rundle, 1988, 1989; Rundle & Klein, 1989; Carlson and Langer, 1989; Bak and Tang, 1989; Rundle & Turcotte, 1993; Sahimi et al., 1993). In the case of earthquake faults, complications such as 1) the long time scales involved (comparable to the human lifetime) and 2) the (at present) unpredictable nature of the events, make earthquakes difficult to study systematically in the field. For this reason, it is advantageous to develop computer methods so that the physics of earthquakes can be studied more easily by simulation in the computer. In that environment, predictions that result from alternative hypotheses can be readily examined. Once predictions are made, tests of the hypotheses can be devised that demand specific observational programs or laboratory experiments. In turn, if interesting and unexplained new observations are made, simulations can be carried out with a view towards explaining them.

Simulations will lead to a much more detailed understanding of the physics of earthquakes, well beyond the present state of affairs in which information exists only on time and space averaged processes. An improved understanding of the physics of earthquake processes would, among other objectives, allow predictions of phenomena such as:

(1) The probability that a given event is a foreshock of a coming mainshock (e.g., Scholz, 1990)

(2) The validity of the "gap" and "antigap" models for earthquake forecasting (e.g., Kagan and Jackson, 1991; Nishenko et al., 1993))

(3) Spatial correlation lengths and temporal correlation times for faults, which would allow the estimation of the effect of activity on one fault upon the activity on the same or other faults (e.g., Rundle et al., 1996)

(4) One or more earthquake forecast algorithms, based upon log-periodic (Sornette et al., 1996), bond probability (Coniglio and Klein, 1980), or other possible methods to be developed.

(5) The physical origin of the "action at a distance", and "time delayed triggering", i.e., seismicity that is correlated over much larger distances and time intervals than previously thought (see for example, Hill et al., 1993)

This large scale project should be viewed as a "testbed" for HPCC technologies. Over the past decade, an increasing number of research studies, described below, have investigated the physics of failure on individual faults. The vast majority of this work was done on single computers, either workstations or supercomputers. This work has now reached a level of maturity at which the next logical step in the research progression is to adopt a "systems approach", to build simulations of systems of faults, and to look at the influences of one fault upon another. In order to make this step, simulations on a scale not previously attempted will be necessary, of the order of solving 106 to 109 simultaneous nonlinear equations. These simulations will demand, for example: 1) significant advances in numerical algorithm efficiency; 2) advances in messaging and high speed communication hardware and software for networking supercomputers, both with each other and with other distributed processors; and 3) the development of ideas about how to deal with physical processes on smaller spatial and temporal scales than is practical to include in the simulations.

Finally, while the methods to be developed in the proposed work will be generally applicable, the initial modeling focus will be on one of the most well-studied earthquake fault systems in the world, southern California. For that reason, close cooperation and coordination is envisaged with the staff and management of the Southern California Earthquake Center. Coordination will be greatly facilitated since Rundle (PI) is Chair of the Advisory Council of SCEC, Jordan (Co-I) is a member of the Advisory Council, Minster (Co-I) is Chair of the Board of Directors, and Stein (Co-I) is an active member of SCEC. Moreover, Rundle's Ph.D. dissertation advisor was Professor David Jackson, who is now the Science Director of SCEC.

Construction of Models. In developing the numerical simulation techniques proposed here, there are many decisions to be made. Foremost among these is the physics to be included in the simulations, and the associated approximations to be used. As a general principle, we plan to examine the simplest physical models that can be shown to capture a given observable effect. We are motivated to take this approach because even the existence of some observational data are often still subjects for research, an example being Heaton-type pulses rather than Kostrov crack-like slip processes (e.g., Heaton, 1990; Kostrov, 1964). In the Heaton model, healing takes place immediately after a site on the fault slips, rather than at a delayed time after slip everwhere is complete, as happens in the Kostrov model.

In complex systems such as earthquake faults, it is often found that a given physical observable may have its origin in a variety of physical processes, some quite unexpected. An example is the question of whether the Gutenberg-Richter magnitude frequency relation arises from self-organization processes on a single fault with no inherent scale (Burridge and Knopoff, 1967; Rundle and Jackson, 1977), or from a scale invariant distribution of faults (Rice, 1993). For single faults, the two types of models are very different. The first (BK model) includes inertia, whereas the second, a massless Cellular Automaton (CA) model, does not. For the massive BK models, the dynamics must be obtained by solving a set of coupled nonlinear differential equations. It was shown by Nakanishi (1990) that massless CA's can be constructed that give the same quantitative results as the massive slider block models examined by Carlson and Langer (1989). These results include not only the rupture patterns, but also the scaling statistics. The decision about whether or not to include seismic radiation is another example. Kanamori and Anderson (1975) estimated that the seismic efficiency h, which measures the fraction of energy in the earthquake lost to seismic radiation, is less than 5%-10%. The absence of inertia does not preclude a time evolution process in the CA, with a separation of loader plate and source time scales, just as in models with inertia and waves. The physical interpretation of the CA time scales is discussed by Gabrielov et al. (1994).

Other decisions include the following. 1) Interactions between fault segments might be short range, 1/r3 "continuum interactions", or long range of some other variety. 2) Friction might be a) simple stick slip (Mohr-Coulomb), b) have a simple velocity dependence as in the original Burridge-Knopoff model, c) be based on idealized laboratory friction experiments on clean, dry, dust-free sliding surfaces such as in the rate and state models, or d) be based on new physical ideas assuming the evolution of systems toward states of lower energy and greater stability (Rundle et al., 1996). 3) How to include "asperities" (Lay et al., 1982) and other regions of heterogeneous friction in the various friction models is an open question. 4) Similarly, how to include unmodeled faults, particularly at small length scales is another open question, but presumably these might be treated as some kind of "noise" in the model. 5) The effects of pore water might be included or ignored, along with the effects associated with variable normal stress on the fault. Which observable data are to be examined will determine the physics to be included in a given simulation.

We refer the interested reader to a recent review of the literature compiled by Rundle and Klein (1995). A brief summary of some relevant material is provided here. While much research in the area of earthquake physics is fundamental in character, the results have many potential applications in the areas of earthquake risk and hazard analysis, and seismic zonation (see e.g., Cornell et al., 1993; Jackson et al., 1995). In an important earthquake, the June 28, 1992 Landers event (Sieh et al., 1993), observed effects (Hill et al., 1993) indicated that microearthquake activity in geothermal areas as remote as 1000 km epicentral distance (Cascades, Yellowstone) were triggered by the sudden onset of stress increments associated with the Landers event. These data may imply very long correlation lengths for fault systems, of the order of a thousand km or more.

A nonevent with important implications was the missing Parkfield, California, M ~ 6 earthquake, which was originally forecast to occur by January, 1993 (Bakun and Lindh, 1984). Retrospective studies have concluded that the pattern of occurrence of the the historic Parkfield events is consistent with considerable variation in event interoccurrence time (Savage, 1991; 1993). This theme was echoed in work on the validity of the seismic gap hypothesis by Kagan and Jackson (1992), who concluded that the dominating effect is the tendency of earthquakes to cluster in time. Space-time clustering of events is an important signature of a nonlinear complex system (Stauffer and Aharony, 1994; Ma, 1985). The tendency for earthquakes to cluster in time is supported also by paleoseismic observations of Grant and Sieh (1994).

Of major significance was the widespread realization that sliding on faults at roughly constant stress drop does not necessarily imply a Kostrov-type slip distribution as mentioned in the preceeding section. The Kostrov distribution arises when no healing takes place until progression of the rupture front stops. By contrast, the "Heaton pulse" model (Heaton, 1990) allows healing just behind the rupture front, and in fact, there is no a priori reason to exclude other more general modes of rupture as well (Rundle and Klein, 1994). Evidence for rupture pulses and non uniform healing behind rupture fronts is now accumulating (Wald et al., 1991; Beroza, 1991; Mendoza et al., 1994).

An important set of observations relates to the existence of "slow" earthquakes precursory to mainshock events that generate seismic radiation, the first of which was observed by Kanamori and Cipar (1974). A number of authors (Jordan, 1991; Ihmle et al., 1993; Kanamori and Hauksson, 1993; Kanamori and Kikuchi, 1993; Kikuchi et al., 1993; Wallace et al., 1991) have found evidence for slow events preceeding mainshocks. The degree to which fault plane structure (asperities and barriers) control slip has been investigated by Beck and Christensen (1991) and Ruff (1992), who provide maps of the spatial distribution of enhanced slip on the fault plane. The physical interpretation of asperities in terms of strong regions on the fault depends critically on assumptions about the constancy of rupture velocity and directivity (Rundle and Klein, 1994).

Recent observational work continues to confirm the importance of fault interactions (Bilham and Bodin, 1992; Hill et al., 1993; Sanders, 1993), which of course also exist between patches on the same fault. These interactions cause slip on the fault to self-organize (Rundle, 1988), so that slip occurs in earthquakes obeying the Gutenberg-Richter relation, and not as random events (the existence of power law behavior is incompatible with truly random dynamics). Other work points out that the existence of scaling implies that small earthquakes are as important as large events in redistributing stress (Hanks, 1992), a fact that is clearly seen in simulations.

Several new models and ideas have had a significant impact on thinking about the physical origins of the observations discussed above. A new Traveling Density Wave model (Rundle et al., 1996; Gross et al., 1996), provides a unifying new framework for friction, earthquakes, critical phenomena and nucleation, and identifies characteristic earthquakes on fault systems as spinodal nucleation events. Other work attempts to put earthquakes into the same category as the self-organized criticality model for sandpiles that has no tuning parameter. Still other researchers have proposed that earthquakes are chaotic sequences of events (Narkounskaia et al., 1992; Huang and Turcotte, 1992).

Fundamental Equations: The basic problem to be solved in GEM can be simply stated (e.g., Rundle 1988a). Given a network of faults embedded in an earth with a given rheology, subject to loading by distant stresses, the state of slip s(x,t) on a fault at (x,t) is determined from the equilibrium of stresses according to Newton's Laws. This can be written schematically as:

sint[x,t; s(x',t'); p] = sf(x,t; s(x,t)) (1)

where 1) sint is the interaction stress provided by transmission of stress through the earth's crust due to background tractions p applied to the earth's crust, as well as stresses arising from slip on other faults at other sites x' at times t'; and 2) sf is the cohesive fault frictional (f) stress at the site (x,t) due to processes arising from the state of slip s(x,t). The transmission of stress through the earth's crust involves both dynamic effects arising from the transient propagation of seismic waves, as well as static effects that remain after wave motion has ceased. For the GEM project, it will suffice to assume for this proposal that, exterior to the fault zones, the earth's crust has a rheology (deformation constitutive relation) that is linear elastic, linear viscoelastic, or linear poroelastic. The first assumes only small strain in regions external to the faults; the second allows rock to deform by bulk flow in response to applied stresses; and the third allows crustal pore water to flow in response to applied stress. The linear nature of all these constitutive laws allows stress from various sources to be superposed, and the assumption of linearity is expected to be appropriate in all crustal volumes except very near a fault. The nature of the frictional stress is the subject of current debate, but several candidate models have emerged that we will discuss below.

There are a few other model techniques that are being used to simulate faults, the primary example being the molecular dynamics simulations of Mora and Place (1994). At the moment, these seem to be primarily useful for simulating dynamics of individual faults, so these will not be considered further.

Green's Functions: Focusing on GEM models that have a linear interaction rheology between the faults implies that the interaction stress can be expressed as a spatial and temporal convolution of a stress Green's function Tijkl (x-x',t - t') with the slip deficit variable f(x,t) = s(x,t) - Vt, where V is the far field motion of the crust driving slip on the fault. Once the slip deficit is known, the displacement Green's function Gijk(x-x',t - t') can be used to compute, again by convolution, the deformation anywhere in the surrounding medium exterior to the fault surfaces (see Rundle 1988a for details).

The three kinds of rheologic models that describe to a good approximation the physical nature of the earth's crust between the faults and different sites on the same fault are all linear in nature (e.g., Rundle and Turcotte, 1993). The three rheologies describe 1) a purely elastic material on both long and short time scales; 2) a material whose instantaneous response is elastic but whose long term deformation involves bulk flow (viscoelastic); and 3) a material that is again elastic over short times, but whose long term response involves stress changes due to the flow of pore fluids through the rock matrix (poroelastic). We plan to use all of these models eventually for calculation of the interactions stresses between faults, sint.

In the first implementation of GEM models, we will further specialize to the case of quasistatic interactions. For this case, the time dependence of the Green's function typically enters only implicitly through time dependence of the elastic moduli (see, e.g., Lee, 1955). Because of the linearity property, the fundamental problem is reduced to that of calculating the stress and deformation Green's function for the rheology of interest. For materials that are homogenous in composition within horizontal layers, numerical methods to compute these Green's functions for layered media have been developed (e.g., Okada, 1985, 1992; Rundle, 1982a,b, 1988; Rice and Cleary, 1976; Cleary, 1977; Burridge and Varga, 1979; Maruyama, 1994). The problem of heterogeneous media, especially media with a distribution of cracks too small and numerous to model individually is often solved by using effective medium assumptions or self-consistency assumptions (Hill, 1965; Berryman and Milton, 1985; Ivins, 1995a,b.). Suffice to say that a considerable amount of effort has gone into constructing quasistatic Green's functions for these types of media, and while the computational problems present certain challenges, the methods are straightforward because the problems are linear. These will not be discussed further here, and the interested reviewer is referred to the literature.

In later GEM models, the inclusion of elastic waves at least approximately, will be important. Stresses are initially transmitted in the earth by means of elastic waves (e.g., Aki and Richards, 1980; Zheng et al., 1995; Kanamori, 1993; Beroza, 1995; Jordan, 1991; Imhle et al., 1993). During an earthquake, these will affect the details of the temporal development of sliding during the event. This is an important effect, and should eventually be included in maximally realistic simulations of earthquakes. Since the time scales over which earthquakes occur are short relative to relaxation times for viscous and porous effects, the effects of waves need only be calculated in linear elastic, attenuating, dispersive media. However, as inclusion of these effects will be severely limited by available computational capability, it is anticipated that it may be practical to include only the longest wavelength, or largest spatial scale, effects. This computational plan is consistent with our philosophical approach to the problem of spatial scales that is discussed below.

Multipolar Representation of Fault Systems: A primary key to a successful implementation of GEM models of systems of faults will be to utilize a computationally efficient set of algorithms for updating the interactions between faults and fault segments. These would normally be of order N*N for times between earthquakes, and possibly worse for segments of faults experiencing an earthquake. As discussed above, the nature of the quasistatic interactions is such that the Green's functions Tijkl and Gijk for linear elasticity have a simple time dependence. Moreover, the Green's functions for linear viscoelasticity and for linear poroelasticity can be obtained easily from the elastic Green's functions using the Principle of Correspondence (e.g., Lee, 1955; Rundle 1982a,b). These simplifications strongly suggest that an approach based on multipole expansions (Goil, 1994; Goil and Ranka, 1995) will be the most computationally efficient algorithm.

In fact, the stress and displacement Green's functions Tijkl and Gijk actually represent the tensor stress and vector displacement at x due to a point double couple located at x' (see for example, Steketee, 1958). The orientation of the equivalent fault surface normal vector, and of the vector displacement on that fault surface, are described by the indices i and j. Integration of Tijkl and Gijk over the fault surface then correspond to a distribution of double couples over the fault surface. For that reason, representation of the stress over segments of fault in terms of a multipole expansion is the natural basis to use for the GEM computational problem. Details of the computational implementation of this approach will be described in the following sections.

C.5.2 Friction Models

Friction Models: The nature of the nonlinear cohesive stress sf on faults is the subject of considerable debate, as will be evident in the discussion below. At the present time, there are six basic friction laws that have been proposed for faults, these are discussed briefly below. Our goal in this proposal is to develop computational methods to allow systematic comparison of the predictions of these friction models in the context of models with multiple, interacting fault systems similar to those found in southern California.

Phenomenological Friction Models:

There are two basic classes of friction models that arise from laboratory experiments. These are 1) Slip weakening models, and 2) Rate and state dependent friction models. The former has a long history that originates from experiments carried out over a number of decades, the latter is a more recent development:

1) Slip Weakening model. This friction law (Rabinowicz, 1965; Bowdon and Tabor, 1950; Stuart, 1988; Li, 1987; Rice, 1993; Stuart and Tullis, 1995) assumes that the frictional stress at a site on the fault sf = sf[s(x,t)] is a functional of the state of slip. In general, sf[s(x,t)] is peaked at regular intervals. The state of slip at a site on a fault in elastic medium, for example, is then obtained by solving the quasistatic equation sint[f(x,t)] = sf[s(x,t)] at every instant of time. Implementations of this model usually use a series of sweeps through the lattice similar to the CA models, successively solving the force balance equation on each sweep. The advantage of this model is that it has a strong basis in an enormous number of laboratory experiments under many different kinds of conditions, on many different sliding surfaces (Rabinowicz, 1965; Bowdon and Tabor, 1950; Beeler et al., 1996). The disadvantage is that the quasistatic equation does not really capture details of the time dependence of the sliding process.

2) Rate and State Dependent Friction Laws. These friction laws are based on certain kinds of laboratory sliding experiments in which two frictional surfaces are slid over each other at varying velocities, without ever experiencing arrest (Dieterich, 1972; 1978; 1981; Ruina, 1983; Rice and Ruina, 1983). The rate dependence of these friction laws refers to a dependence on sliding velocity somewhat the same as in the velocity weakening laws described above, although in this case the dependence is on the logarithm of the velocity. The state dependence refers to a dependence on one or more state variables qi(t), each of which follows its own relaxation equation. The effect of both the velocity and the state variable dependence is to introduce a "memory effect", so that changes in sliding velocity are reflected in friction variations only after a certain time delay. This friction law has been implemented in simulations of single faults (Dieterich and Linker, 1992; Rice, 1993). The advantage of using this friction model is that it's formulation is based on laboratory experiments, however idealized. Among the disadvantages are the lack of understanding of the physical interpretation of the state variables, particularly in the multi-state variable formulation. In addition, the laboratory experiments are not particularly representative of the physics of sliding on faults, and certainly not of inertially dominated sliding. These problems lead to predictions from the simulations that are never observed in nature. Among these predictions are 1) infinite friction at arrest of the sliding surfaces; 2) high shear stress and heat generation on the fault plane; 3) only periodic large events ever occur; and 4) significant nonzero precursory and interseismic quasistatic slip between events. More recent work by Beeler et al. (1994) and Beeler et al. (1996) demonstrates that sliding surfaces do experience arrest, and this substantially modifies the friction law and the associated results.

Models Based on Computational Simplifications:

These models are essentially simplifications of the two kinds of models described above, the major benefit being that they are much simpler to implement computationally. The two major models in this category are 1) Cellular automaton model, and 2) Velocity weakening model. The former is the easiest to implement inasmuch as only simple update rules are needed, the latter was first proposed by Burridge and Knopoff (1967):

1) Cellular Automaton (CA) models. These are the simplest ("freshman physics") models to implement on the computer, and for that reason have been widely used (e.g., Rundle and Jackson, 1977; Nakanishi, 1991; Brown et al., 1991; Rundle and Brown, 1991; Rundle and Klein, 1992). A static failure threshold, or equivalently a coefficient of static friction mS is prescribed, along with a residual strength, or equivalently a dynamic coefficient of friction mD. When the stress on a site increases, either gradually or suddenly, to equal or exceed the static value, a sudden jump in slip (change of state) occurs that takes the stress at the site down to the residual value. The exact details of the jump process, and whether it occurs synchronously or asynchronously with the jumps at other sites, is an arbitrary decision. These models assume a separation of time scales, the time scale for long term motion being far larger than the time scale on which adjustment of the slip state occurs. To implement the CA dynamics, the time is calculated at which the least stable site becomes unstable due to far field motion. Time is advanced to this point, after which it is "frozen" until slip is complete. The slip state of the unstable site jumps to a new value according to the prescribed rule, after which a sweep is conducted through the lattice to find other sites that have become unstable (stress larger than failure threshold) as a result of direct elastic interaction with the first site. Any unstable sites are also then allowed to jump synchronously, after which further sweeps are carried out, followed by more jumps, etc. The result is an avalanche of failed sites that are analagous to an earthquake. While the advantage of this friction model is ease of computational implementation, the disadvantage is it captures the physics of stick-slip sliding only crudely, and allows considerable lattitude in specification of the jump mechanism.

2) Velocity Weakening model. This model (Burridge and Knopoff, 1967; Carlson and Langer, 1989) is based on the observation that frictional strength diminshes as sliding proceeds. A constant static strength sf = sF is used as above, after which the assumption is made that during sliding, frictional resistance must be inversely proportional to sliding velocity. In the implementation of this model, a differential equation is written expressing the balance of forces to be solved. The advantages of this model are that the details of the sliding process can be calculated to arbitrary numerical precision, and that inertia of the fault, although not waves, are included. The disadvantages are the ad hoc nature of the friction law, and that practical numerical solutions of these equations are limited to one-dimensional models with nearest neighbor interactions having only a few hundred sites.

Models Based on Statistical Mechanics

These models for friction are based on the fundamental principles that the dynamics of faults and fault systems can be described by the construction of probability distributions involving the physical variables that characterize stress accumulation and failure. The basic goal is to construct a series of nonlinear stochastic equations whose solutions can be approached by numerical means, and which have embedded within them the kinds of statistical distributions that are often observed in nature. An interesting result is that broad classes of these models have been found that, at least at the meanfield level, predict essentially the same kinds of Ito-Langevin equations for the dynamics. Recent work (Ferguson et al., manuscript in preparation, 1996) indicates that most mean field models can be reduced in a generic way to Ginzburg-Landau type functionals, giving rise to metastability and nucleation processes similar to well-understood equilibrium systems. These nonequilibrium models seem to imply that time can be regarded as a kind of scaling field, giving a strong theoretical underpinning to the observation of broad classes of temporal scaling laws. The two models in this category are 1) Traveling density wave models, and 2) Heirarchical statistical models:

1) Traveling Density Wave Model. This new model (e.g., Rundle et al., 1996; Gross et al., 1996) is based on the same experiments that led to the formulation of the slip weakening model. The interaction stress sint is assumed to derived from elasticity, sint = se. The equality between elastic stress based on slip deficit and and frictional stress based on slip, se[f] = sf[s] can be trivially rewritten as se[f] = sf[f + Vt], indicating that the solution lies at the intersection of the curve se[f] and a leftward moving wave sf[f + Vt]. This equation is then regarded as the Euler-Lagrange equation derived from a Lyapunov functional U[f]. The principle of evolution towards maximum stability then implies a kinetic equation in which the rate of change of f is proportional to the functional derivative U[f]. Temporal changes in the friction law itself, related to wearing of the frictional contact surfaces is then modeled by random changes in the friction parameters, such as the energy of cohesion. The advantages of this model is that it is based on a large body of data in the tribology literature, and that it represents the kind of statistical field theory about which a great deal is known, specifically the scaling properties and nucleation rates. The primary disadvantage at the moment is that the equations are of the nonlinear integral-differential kind that are known to display chaos and intermittancy. Thus solution schemes for these equations will be difficult to implement computationally. It should be pointed out that recent work on CA models described above (Rundle et al., 1994; Klein et al., 1995) has shown that the internal energy of CA models with long range interactions ("mean field" models) are characterized by Boltzmann fluctuations. Coarse graining of mean field CA models in space and time yields a Langevin equation for the evolution of slip deficit f(x,t) that is essentially identical to that obtained for the TDW model.

2) Hierarchical Statistical Models. Classes of hierachical, statistical models have been developed recently. Examples of these are the models by Allegre et al. (1982); Smalley et al. (1985); Blanter et al. (1996); Allegre and Le Mouel (1994); Allegre et al. (1996); Heimpel (1996); Newman et al. (1996); and Gross (1996). These are probabilistic models in which hierarchies of blocks or asperity sites are assigned probabilities for failure. As the level of external stress rises, probabilities of failure increase, and as a site fails, it influences the probability of failure on nearby sites. The advantages of these kinds of models is that they are frequently easy to program on a computer, have a simple, intuitive appeal, and often can be shown to lead to a scaling behavior, both for the events themselves, and for aftershocks and foreshocks (Omori's law). The disadvantages lie in the need to arbitrarily choose the failure probabilities, and the corresponding large number of unknown parameters.

C.5.3 Symmetries and Scaling Properties

Heirarchy of Spatial and Temporal Scales: The presence of heirarchies of spatial and temporal scales is a recurring theme in simulations of this type. It is known, for example, that fault and crack systems within the earth are distributed in a scale invariant (fractal) manner over a wide range of scales. Moreover, the time intervals between characteristic earthquakes on this fractal system is known to form scale invarient set (Allegre et al., 1982, 1994 1996; Smalley et al., 1985; 1987). Changes in scaling behavior have been observed a length scales corresponding to the thickness of the earth's lithosphere (e.g., Rundle and Klein, 1995), but the basic physics is nevertheless similar over many length scales. It also is known that nucleation and critical phenomena, the physics that are now suspected to govern many earthquake-related phenomena, are associated with divergent length and time scales and long range correlations and coherence intervals (see, e.g., Rundle and Klein, 1995 for a literature review and discussion). For that reason, our philosophical approach to simulations will begin by focusing on the largest scale effects first, working down toward shorter scales as algorithms and numerical techniques improve. Moreover, our practical interest is limited primarily to the largest faults in a region, and to the largest earthquakes that may occur. Therefore, focussing on quasistatic interactions and long wavelenth and long period elastic wave interactions is the most logical approach. We plan to include smaller faults and events as a kind of background "noise" in the simulations, as discussed in the proposed work.

Renormalization Group Analysis: This topic is closely connected to the previous discussion of spatial and temporal scales. Many of the models that have been proposed (see e.g., the review by Rundle and Klein, 1995); Rundle et al. (1996); Klein et al. (1996); Ma (1976), Newman and Gabrielov (1991); Newman et al. (1994); Newman et al (1996); Stauffer and Aharony (1994); Turcotte (1994); are at some level amenable to analysis by Renormalization Group techniques. At the very least, an RG analysis be used to identify and analyze fixed points and scaling exponents, as well as to identify any applicable universality classes. Thus, these techniques can clarify the physics of the processes, and in particular identify the renormalization flows around the fixed point(s), which will help to provide physical insight into the influence of the scaling fields, and the structure of the potential energy surface. Friction models for which this could be particularly helpful include the hierarchical statistical models, the CA models, and the traveling density wave model, where there are expectations of the presence of fixed (critical) points.

For the most part, research during the past decade has focused on models of single faults embedded in an ideally elastic medium. Physics of failure processes, particularly the effects of heterogeneous friction on the fault surface, have been of interest. In this case, the interactions given by the elastic stress Green's function have a relatively simple spatial dependence. The more interesting case of a three dimensional array of faults is associated with a highly directionally dependent interaction. The physics of these systems will be extremely rich, inasmuch as it is quite possible for frustrated interactions and screening effects to arise from the interactions between the various faults. Also, relatively little work has been carried out using inelastic stress Green's functions. The current state of the art for three dimensional fault systems with inelastic interactions is still the work of Rundle (1988b).

Models of Single Faults

A major theoretical theme in recent years has been the development of increasingly sophisticated computer simulation techniques for understanding the dynamics of single faults. Recent observational work continues to confirm the importance of fault interactions (Bilham and Bodin, 1992; Hill et al., 1993; Sanders, 1993), which of course also exist between patches on the same fault. These interactions cause slip on the fault to self-organize, so that slip occurs in earthquakes obeying the Gutenberg-Richter relation, and not as random events (the existence of power law behavior is incompatible with truly random dynamics). Other work points out that the existence of scaling implies that small earthquakes are as important as large events in redistributing stress (Hanks, 1992), a fact that is clearly seen in simulations.

While a few authors use some version of continuum intereactions (stress Green's functions) in their models (Rice, 1993; Ben-Zion and Rice, 1993; Ben-Zion et al., 1993; Stuart and Tullis, 1995; Taylor et al, 1996; Ben-Zion, 1996), many more use some version of a nearest-neighbor slider block model. Rice (1993) in fact argues that nearest-neighbor models are inappropriate. However, it is known from extensive studies over many years on the statistical mechanics of spin dipole systems (Ma, 1985) that nearest-neighbor models such as Ising systems display much of the same physics as models with dipole (1/r3) interactions. This is a result of the large correlation lengths that appear near a critical point, in association with power law scaling like the Gutenberg-Richter relation. Examples of slider-block models include Carlson et al. (1991); Brown et al. (1991); Carlson (1991); Narkounskaia and Turcotte (1992); Narkounskaia and Turcotte (1992); Huang et al. (1992); Knopoff et al. (1992); Shaw et al. (1992); Vasconcelos et al. (1991); Olami et al. (1992); Rundle and Klein (1993); Ding and Lu (1993); de Sousa et al.(1993); Lomnitz-Adler (1993); Pepke et al. (1994); Rundle and Klein (1994).

Some of these papers attempt to put earthquakes into the same category as the self-organized criticality model for sandpiles that has no tuning parameter. Others show that earthquakes are more probably an example of self-organizing systems that can be tuned to approach a critical point, the different physics implying membership in a different universality class and different pattern formation mechanisms. These models typically use concepts from percolation theory to analyze the simulation data. Sahimi et al. (1993) use the percolation idea directly to show that large earthquakes are associated with the backbone of the percolation cluster of fractures forming the fault system in the earth. Others found chaos in the dynamics of their models (Narkounskaia et al., 1992), even in low order systems (Huang and Turcotte, 1992). Models can in principle also be used as testbeds to develop techniques for earthquake forecasting. Ward (1992), Dowla et al. (1992) and Pepke et al. (1994) suggested several methods for testing the predictability of synthetic earthquake sequences, using for example neural networks as well as pattern recognition.

Models of Interacting Fault Systems

While models of single faults can be used to compute surface deformation associated with failure processes (e.g., Stuart and Tullis, 1995), only much more realistic models of three dimensional interacting fault systems can approach the kind of detailed predictions necessary in using actual data. This category of model includes that of (Stuart ; 1986; Rundle, 1988b; Ward, 1992; and Ward and Goes, 1993).

In fact, the types of data from real fault systems to which the GEM simulations will be compared include all the same data that characterize natural fault systems (see, e.g., Kanamori,

Real and model faults in southern California (see Rundle, 1988 for details of model calculations)

Space-time diagram for model seismicity, and model friction, as a function of distance along model San Andreas fault. Dark portion in middle is creeping part along central San Andreas fault in Carrizo Plain. Friction parameters in Big Bend section were optimized to reproduce events similar to real data (see figure 3 below).

Comparison of real and simulated earthquakes. Slip as a function of distance along the fault is shown for several model and real events. note similarity between real 1857 event, and simulated event at model year 9707. Data from Sieh & Jahns (1984)

Time dependent deformation field associated with model event at model year 9707.

Comparison of real surface strain field in southern California and strain field calculated from simulations. Simulation strain field is calculated at a time corresponding to relationship between 1857 event and observed strain field observed today (e.g., about ~130 years later). Data from Savage et al. (1986)

1993). This includes, for example, 1) earthquake histories obtained from paleoseismic investigations, other geological data and historical records; 2) fault creep data obtained from creepmeters and strain meters; 3) the temporal development of slip on faults as determined by local, regional, and teleseismic data, including strong motion instruments; 4) spatial and temporal seismicity patterns, together with scaling statistics such as Gutenberg-Richter distributions, Omori laws, and foreshock and aftershock statistics; 5) surface deformation data obtained from historical geodetic techniques, land-based methods such as multiwavelength trilateration, and newer space-based methods such as Global Positioning System technology and Synthetic Aperature Radar Interferometry; 6) in situ data such as stress and dilatometer measurements, pore fluid pressure, and radon emission rates; and 7) earthquake source mechanism compilations.

As an example of the approach described above, and the synthetic data that may be computed, figures 1- 5 show results obtained by Rundle (1988) in a crude simulation of the San Andreas fault system using 80 fault segments embedded in an elastic lithosphere overlying a viscoelastic half space. This model used a modified CA algorithm for a friction law, in which changes in fault normal stress were included in the CA version of Amonton's law of friction. Additionally, the geometric configuration of the faults is held fixed in time, which for short time intervals is probably a good first approximation. The far field displacements rates are obtained from rigid plate motion models. The simulations clearly show that field data can be simulated to as realistic an extent as is desired for comparison with natural, observed data.

The collection of new observational and laboratory data is not the focus of this proposal, although there may be some secondary need to verify a laboratory result, or to add a new field observation to the data base. Instead, data is viewed here as a means for validating the simualtions. The GEM team expects, however, that recommendations for new data collection activities will emerge as a natural outgrowth of the simulations, and that an interesting new feedback loop will emerge between observationalists and modelers as a result of the GEM project.

Management of Existing Earthquake Data

Primary responsibility for earthquake data collection and archiving lies with the Southern California Earthquake Center (SCEC) at the University of Southern California, the Seismological Laboratory of the California Institute of Technology, and the Pasadena field office of the United States Geological Survey. Data in these archives are of several kinds, 1) broadband seismic data from the TERRASCOPE array owned and operated by CalTech, 2) continous and "campaign style" geodetic data obtained through a number of mechanisms such as trilateration, satellite laser ranging, VLBI, and GPS, 3) paleoseismic data obtained from historic and field trenching studies detailing the history of slip on the major faults of southern California, 4) near field strong motion accelerograms of recent earthquakes (last two decades), 5) field structural geologic studies and reflection/refraction data indicating the orientations and geometries of the major active faults, 6) other sparsely collected kinds of data including pore fluid pressure, in situ stress, and heat flow data. These data will be used, for example, to update the fault geometry models developed under the GEM project, and to update fault slip histories that can be used to validate various earthquake models developed under GEM.

Primary responsibility for interacting with elements of this data base will be given to Kanamori, Jordan, Stein, and Minster.

New Data: Stress Analysis of the Earth from Synthetic Aperature Radar Interferometry

A number of SAR missions are curently taking data over the southern California region. These include the C-band (5.8 cm) French ERS 1/2 satellites , and the L-band Japanese JERS satellite. These missions have already produced revolutionary images of the complete deformation fields associated with major earthquakes in the United States (e.g., Massonet et al., 1993) and Japan, and much more is on the way. Moreover, there are at this time several proposals before NASA to initiate a US SAR mission. A SAR instrument as currently envisioned by the investigator team led by Minster and upon which Rundle and Turcotte serve will provide the definitive technology for testing GEM models, and will lead to a "quantum jump" in our ability to monitor and perhaps forecast these events. As currently configured, the mission will have the advantages of 1) all weather imaging; 2) with Left/Right (L/R) radar capability we can image any spot on the earth in less than 2 days and probably 1 day; 3) capabilities including 24 day exact repeat, 3 swaths every 8 days, ascending/descending orbit repeat of 4 days, L/R repeat of 2 days. A major advantage of this mission will be extremely frequent repeats, so that fine space-time details of major events can be studied. Optimal use of these capabilities depends on having the ability to send instructions to the spacecraft, i.e. short response times, and the ability to download data and make interferograms rapidly.

For surface change missions to study earthquakes, deformation maps produced by SAR Interferometry will allow unprecedented levels of detail. Current models of deformation are very crude, and severely limited in detail due essentially to the Uniqueness Theorem of linear elasticity. This theorem asserts that displacements and forces everywhere on the boundary of an elastic material must be known in order to uniquely infer the forces and displacements within the medium that caused those surface effects. With GPS, only a few tens or at most, hundreds of points will ever be available. However, ISAR images will transform the field from "data poor" to "data rich", and it will be possible to study details of earthquakes as never before. In other words, we will be able to carry out "stress analysis of the earth", similar to the stress analysis methodologies employed by civil and mechanical engineers to study materials and structures. This change data will in fact be the most important dataset with which to test GEM-type models. The extremely fine details of these models demand an equally highly detailed space-time data set with which to test the models. For example, current ERS images taken over the Los Angeles basin clearly show deformation associated with fluid withdrawal-induced subsidence and other natural and man-induced processes (G.Peltzer, personal comm., 1996). Other SAR data sets recently compiled clearly show evidence of crustal pore fluid-induced deformation following the June 1992, Magnitude ~ 7 Landers earthquake. The extent of these processes must be known in order to test and modify the GEM suites of models.

Primary responsibility for interacting with the SAR data will be given to Minster, Turcotte and Rundle. Again, Minster is the team leader for the NASA Earth System Science Pathfinder proposal to orbit an L-band Interferometric Synthetic Aperature Radar in the near future, and Rundle and Turcotte are team members.

C.8 Overview of Computational Science Issues

The GEM Problem Solving Environment (PSE) requires several technology components which we discuss in the following. One major component is the detailed simulation modules for the variety of physics and numerical approaches discussed in the preceding sections. These modules will be set up in an object oriented framework using the Los Alamos POOMA system. The fast multipole, statistical mechanics and cellular automata subsystems will all need state of the art algorithms and parallel implementations. However the overall PSE will need a coarse integration to enable both distributed collaboration on the simulations and to link simulation, data analysis and visualization. In the following we discuss these various aspects of the GEM PSE, starting with the basic numerical formulation to set the scene.

Section C.5 presents an overview of the scientific framework for GEM. The basic degrees of freedom in GEM correspond to fault segments with the independent variables corresponding to the slip vector s(i,t or x,t) or more precisely the slip deficit f(i,t or x,t). Here the slip deficit represents the deviation from the average long term movement of the fault. In the original (Rundle 1988) calculation, there were N=80 fault segments which averaged 30 km (horizontal) by 18 km (vertical) in size. In a maximally realistic simulation, one would expect to need time and position varying segment sizes with 100 meters as a initial goal for size in critical regions. This indicates that one needs several (perhaps up to 100) million segments and so scopes the problem as of Grand Challenge proportions! In fact a 10m meter resolution (corresponding to local magnitude of about 0) is likely to be important with a corresponding two order of magnitude increase in the number of degrees of freedom.

In more recent simulations of long range single fault models (Rundle et al., 1996; Klein et al., 1996; Ferguson, 1996), networks of sites much larger than Rundle's original simulations have been modeled. Runs of 100 x 100 sites, with individual sites interacting with 212 -1 neighbors, have been carried out for 10,000 earthquakes at a time. These runs typically take of the order of several hours on a Pentium 150 MHz single processor desktop machine. Models with 256 x 256 sites have run for 10,000 earthquakes in which all sites interact with each other. Making use of a simple Fourier transform algorithm, these runs take several hours to complete on a CM-5 128 node processor, or current Boston University Power Challenge system. In real time, this would correspond to as much as several years of seismicity on a fault system the size of southern California.

In section C.5, the general framework for the problem is formulated as a set of Green's function equations for the stress and displacements at the fault segments. It is convenient to "think" about fault segments as "particles" interacting with a long range force which originates from a "double-couple", i.e. a pair of equal and opposite dipoles. Initially we will ignore time dependence but this approximation will be relaxed (in year 3) with the explicit inclusion of radiative time dependence. This formulation corresponds to a set of N ordinary differential equations -- one for each "particle" i (i.e. fault segment) -- where each right hand side has a sum of N terms corresponding to the long range force structure. At this level, we have a situation like conventional the O(N2) N-body problem but there are many important differences. For instance, the critical dynamics -- namely earthquakes -- are found by examining the stresses at each time step to see if the friction law implies that a slip will occur. As discussed earlier, there are many different proposals for the friction law and the computational system needs to be set up to be able to flexibly change and compare results from different friction laws. The analogies with statistical physics are seen by noting that earthquakes correspond to large scale correlations with up to perhaps a million 10-to-100 meter segments slipping together. As in critical phenomena, clustering occurs at all length scales and we need to examine computationally this effect. We find differences with the classical molecular dynamics N body problems not only in the dynamical criteria of importance but also in the dependence of the Greens function (i.e. "force" potential) on the independent variables. Thus we expect that a major algorithmic effort will be needed to examine succesful approaches to the classic N body problem and adapt them to GEM. Colorado and Syracuse will play a key role in this aspect of the project.

There are many important and distinct ways of tackling the basic GEM equations and these must all be supported by the GEM Problem Solving Environment. Note that a sophisticated fully operational fast multipole simulator does not mean that simpler approaches become irrelevant. Rather one still needs approximate methods such as cellular automata, as these provide important insight as to the fundamental mechanisms underlying earthquakes. One can use either direct time evolution or Monte Carlo methods and the analogies with statistical mechanics suggest the latter will be valuable although current simulations have not investigated this. However the GEM computational framework must support both possibilities. In the following three subsections, we describe the different computational issues in the basic GEM simulation system. The final subsection describes the overall system integration.

This section describes two important approaches that allow improved intuition and simplification by using formulations familiar from other fields. These need to be supported by the GEM PSE and will be a major focus of the initial simulations while the new multipole methods are being developed. In the "Slider-Block" approach one only keeps the "nearest-neighbor" interactions and this approximation can be generalized to a screened approximation where Green's function elements are ignored when the interacting segments are separated by an appropriate distance. With such a screened force, the computation becomes O(N) of course and not the O(N2) of traditional long range force case. In the "cellular-automata" approach (used in Rundle 1988), one uses a rather large time step ( 1 year in Rundle 88) and evolves the "lattice" of fault segments in discrete steps. In this approach, one sees a formulation very similar to that of statistical mechanics with an earthquake in GEM corresponding to a long range correlation-induced cluster in the statistical physics case. Clustering problems are rather hard to parallelize (we have experience with this from Swendsen Wang and other statistical mechanics problems as described in Fox et al., 1994) There are often not enough points in the cluster (earthquake) to allow much parallelism. In GEM this depends on total number of clustered points at each time step -- it must be either zero or very many to allow efficient parallelization. As earthquakes (clusters) are not uniform in physical domain (at given time step), one must use some sort of cyclic (block scattered) decomposition so that each processor participates roughly equally in any given earthquake scenario. The parallelism is hardest in the case of a screened or O(N) force computation. As in the long range force version, one can get parallelism in force computation, and more generally, the computation of forces (which is easy to parallelise) becomes a larger part of the computation time. The level of effort needed in parallel cellular automata methods will depend on experience gained. We do know how to parallelize the general fast multipole method and this may be the approach of choice on large systems with the CA methods used on smaller SMP's without explicit parallelism. There are some important software differences between CA and fast multipole methods. The former does not typically need sophisticated data structures and so approaches such as HPF have been effectively used in this problem class.

General N Body Problems: In an N-body simulation, the phase-space density distribution is represented by a large collection of "particles" which evolve in time according to some physically motivated force law. Direct implementation of these systems of equations is a trivial programming exercise. It is simply a double loop. It vectorizes well and it parallelizes easily and efficiently. Unfortunately, it has an asymptotic time complexity of O(N2). As described above, each of N left-hand-sides is a sum of N-1 right-hand-sides. The fact that the execution time scales as N2 precludes a direct implementation for values of N larger than a few tens of thousands, even on the fastest parallel supercomputers. Special purpose machines (such as GRAPE) have been succesfully used but these do not seem appropriate in an evolving field like GEM while we are still in the prototyping phase

General Multipole Methods: Several methods have been introduced which allow N-body simulations to be performed on arbitrary collections of bodies in time much less than O(N2), without imposition of a lattice (Appel, 1985; Barnes and Hut, 1986; Greengard and Rokhlin, 1987; Anderson, 1992). They have in common the use of a truncated expansion to approximate the contribution of many bodies with a single interaction. The resulting complexity is usually determined to be O(N) or O(N log N) which allows computations using orders of magnitude more particles. These methods represent a system of N bodies in a hierarchical manner by the use of a spatial tree data structure. Aggregations of bodies at various levels of detail form the internal nodes of the tree (cells). Making a good choice of which cells interact and which do not is critical to the success of these algorithms (Salmon and Warren, 1994). N-body simulation algorithms which use adaptive tree data structures are referred to as "treecodes" or "Fast Multipole" methods. By their nature, treecodes are inherently adaptive and are most applicable to dynamic problems with large density contrasts, while fast multipole methods have been mostly non-adaptive and applied to fairly uniform problems. It is likely that this distinction will blur in the future, however.

The fundamental approximation employed by a gravitational treecode may be stated as:

S j \F(G mj dij,| dij |3) » \F(G M di cm,| di cm |3) + higher order multipoles

where di cm = xi - xcm. Here xcm is the centroid of the particles that appear under the summation on the left-hand side. These Newtonian or Coulomb interactions are a special case of a more general formulation that can be applied to any pairwise interaction or Green's function. Even short-range interactions, i.e., those which can be adequately approximated as vanishing outsidthe GEM collaboration have established the fundamental mathematical framework for determining error bounds for use with arbitrary force laws (Salmon and Warren, 1994,Warren and Salmon, 1995). This error bound is essentially a precise statement of several intuitive ideas which determine the accuracy of an interaction. The multipole formalism is more accurate when the interaction is weak, when it is well-approximated by its lower derivatives, when the sources are distributed over a small region, when the field is evaluated near the center of the local expansion, when more terms in the multipole expansion are used, or when the truncated multipole moments are small.

Multipole Algorithms and Computational Approach: One of two tasks generally takes the majority of the time in a particle algorithm: (1.) Finding neighbor lists for short range interactions. (2.) Computing global sums for long-range interactions. These sub-problems are similar enough that a generic library interface can be constructed to handle all aspects of data management and parallel computation, while the physics-related tasks are relegated to a few user-defined functions. One of our primary computational research objectives will be to define and implement a series of generally useful, independent modules that can be used singly, or in combination, to domain decompose, load-balance, and provide a useful data access interface to particle and spatial data. In combination, these computational modules will provide the backbone for the implementation of our own (and hopefully others) algorithms. These algorithms imply very complex data structures and either direct message passing or parallel C++ (supported by POOMA) is the natural software model. In contrast , the direct N body methods are straightforward to implement in HPF.

Application of Fast Multipole Methods to GEM: Fast multipole methods have already been applied outside the astrophysics and molecular dynamics area. In particular the Caltech and Los Alamos groups have successfully used it for the vortex method in Computational Fluid Dynamics. There are several important differences between GEM and current fast multipole applications which we briefly discussed in section C.7.1. We believe a major contribution of this Grand Challenge will be an examination of the software and algorithmic issues in this area, as well as showing how good modular object oriented design can allow a common framework across many fields. An area of importance which is still not well understood in current applications will be use of spatially dependent time steps with smaller values needed in earthquake regions. An important difference between true particles and the GEM case is that in the latter, the fault positions are essentially fixed in space. Of course a major challenge in astrophysics is the time dependent gravitational clustering of particles. It is not clear if one can exploit this simplification in the case of GEM.

C.8.5 POOMA and the Overall Software Support

The Parallel Object Oriented Methods and Applications framework encapsulates the complexities of data distribution and interprocessor communication beneath a high-level programming interface which maps well onto the equations being solved. The computational scientist, using classes representing the relevant high-level physics abstractions, can focus on designing numerical algorithms, running simulations, and analyzing results. Computer scientists and experts on parallel computing can focus independently on the implementation layers beneath. This independence of programming interface and internal workings allows applications to be optimized without changing the high-level code.

Classes for representing fluid and electromagnetic fields discretized on meshes; classes for representing particles in monte-carlo, particle-in-cell, and molecular dynamics simulations; and the interactions of particles with fields on meshes have already been implemented in the POOMA framework. These all have parallel processing on distributed-memory machines or SMP's included in their design; the basic unit of parallelism--the virtual node--provides a mechanism for load balancing (there can be, and typically are, multiple virtual nodes per physical compute node). Usage and extension of these classes, coupled with current ASCI development toward adaptively refined meshes, will facilitate development of the specialized Seismic framework.

In somewhat more detail, Fast multipole techniques work by clustering the particles into regions of space defined by an octree. This data structure maps well onto the POOMA framework's virtual nodes (see below). Each v-node represents a region of space and manages the containers for whatever particles are in that region. As particles move from the domain of one v-node to another they are transferred to that v-node. If the other v-node is on another processor of a parallel computer, it is packed up, buffered with other particles moving to that processor, and communicated. This will employ technology already in use by the POOMA team.

Because of the object-oriented structure of the framework, the job of handling the particles on v-nodes is separated from how v-nodes are created and assigned to processors. This will allow us to use the techniques developed by Warren and Salmon for distributing an octree near-optimally on a massively parallel computer and efficiently coordinating the communication between v-nodes.

Much of the dipole dynamics will then be handled similarly to the molecular dynamics and particle-in-cell work that is currently being done with the framework. The main new element here will be the inclusion of the algorithms designed by Warren and Salmon. The first development atop the POOMA framework will entail casting the algorithms for:

1) Calculating the layout of an octree of v-nodes across the processors of a parallel computer

2) Coordinating the communication of messages between the v-nodes

3) Organizing the order of the computation and the replication of data so that the calculations may proceed in parallel as much as possible

With this infrastructure in place, capabilities specific to seismic modelling will be implemented including:

1) Simulations with adaptively refined meshes,

2) Particle-grid interactions,

3) Fast-time behaviour,

4) System subcycling, and

5) Inhomogeneous modelling.

Because the framework was designed to allow users to provide their own load balancing techniques, introducing these methods should be straightforward. Building this sort of simulation tool on top of the POOMA framework has the potential benefit of allowing us to easily do grid-based calculations at the same time as particle-based ones. This could be useful for visualization purposes or for working with hybrid algorithms that supplement or modify the dipole interactions based on PDE's solved on the grids. The important thing to note is that although sophisticated techniques are being developed and deployed to handle the parallel aspects of multipole in POOMA, the application developer is presented with an interface which represents a global physics structure (for example, a species of particles or field tensor). In this manner, the application developer is able to focus on the physics issues of a given model or simulation without having to manage the extra complexities of parallelism.

C.8.6 I/O, Visualization, and Data Display and Analysis

By designing I/O in terms of objects, the I/O interface for the computational scientist is simple and high-level. Since the high-level classes in the framework represent the physics and numerical abstractions of the problem domain--for example, vector fields discretized on a mesh of nodes--this interface is also a natural one. It also permits us to optimize the design of the underlying implementation details independently, without changing the high-level application code's interface to I/O. Our work to date has focussed on implementation layers using portable data formats such as Unidata's NetCDF (built atop special parallel I/O substrates). Having the end output in the form of collections of files in a standardized, machine-independent format makes it easier to use different tools for postprocessing such as generic reader modules for visualization systems. We will implement the I/O interface for the Seismic framework's classes along with the rest of the classes' interfaces, leveraging off our I/O development for ASCI and other applications. This will fit well into our ongoing collaboration with LANL ASCI scientists working on the general problem of parallel I/O and data sharing between large-scale simulations.

In Years 2 and 3 of the Proposal, it will be essential to build the overall infrastructure to allow the integration of data analysis and simulations for GEM. Giles at Boston and Syracuse will have major roles here. Although important, we expect that an approach using the Web will be very effective and require modest funding from the proposal. The general approach is described in two online papers (Fox et al, 1995, Fox et al., 1996) but we do not feel it is appropriate to go into detail as we believe that the Web is evolving so rapidly that detailed plans at this stage are unrealistic. However we can indicate some general characteristics of this systems integration problem. The Web commercial industry will supply the basic collaboration environment with digital video and audio support of videoconferencing and a rich variety of tools including common whiteboards, databases etc. We can use Java to develop high level visual interfaces to the GEM resources and allow user to specify the necessary linkage. Protypes of this can be seen at http://www.npac.syr.edu/projects/webbasedhpcc/index.html at the NPAC web site.

Approach: A Problem Solving Environment for the GEM Grand Challenge project will be constructed to accomodate a series of increasingly realistic and complex levels of simulations using novel parallel algorithms. Visualizations will be particularly important inasmuch as observational data sets must be overlain on simulations, and methods must be developed that will allow an "intuitive" approach to data analysis as well as more rigorous statistical methods of analysis. With each new level in simulation complexity, more spatial and temporal details will be included than was the case for prior simulations, as well as more realistic rheologies, more detailed friction models, and so forth. Thus, simpler models can be understood, evaluated and discussed before proceeding to the next level. This approach will be supported by the POOMA object-oriented framework developed at Los Alamos. The physics included in a given set of simulations will be discussed and decided upon by the investigator group, or large subsets of the investigators. Throughout this work, close coordination will be pursued with the Southern California Earthquake Center, which Rundle currently advises as Chair of the Advisory Council of SCEC.

C.9.1 Basic Model

Level 0 simulations will be built upon constructing algorithms for the following

model sets:

1) Major fault systems set in a "meanfield" background of smaller secondary faults,

with reference to the most recent observed fault geometries. Fault geometry

will not evolve in time initially. Faults will be predominantly vertical strike

slip faults for southern California

2) Linear rheologies (isothermal elastic and viscoelastic) for the

upper crust of the earth.

3) Simple algorithms for implementation of friction models, including Cellular

Automaton, Traveling Density Wave, and Rate and State friction

Note that this will still allow a moment release function to be computed for

simulated earthquakes.

4) Neglect of inertia in the sliding process.

5) Neglect of faults below a spatial scale to be determined

Most of these assumptions are the same or very similar to the 80 fault model as implemented in simulations carried out by Rundle (1988).

Level I and higher level simulations will include the following to varying

degrees. Note that all of these additional effects involve questions of science that to varying degrees are still fundamental research questions.

1) Evolving fault geometries: shear and tensional fracture problem

2) Different and more realistic frictional physics

3) Inclusion of the dramatically increasing number of smaller fault

systems as a kind of "noise" in the simulations

4) Inclusion of radiative effects: Inertia and elastic waves

5) Inclusion of other rheologic effects, including pore fluids,

6) Thermal effects, and possibly chemical processes, particularly as related

to stress corrosion

Scientific and Computational Problems: In all of the simulations and model activities, the GEM team will focus on a basic set of scientific and computational problems as listed below. These are generic problems which are applicable to all GEM type simulations, thus they will recurr repeatedly.

Generic GEM scientific problems include:

1) Developing the procedures and mechanisms for using a General Earthquake Model (GEM) methodology for use as a focus in organizing the data collection and interpretation activities in the field of earthquake science.

2) Evaluating whether forecasting of earthquakes at some level is at all possible, and how this depends on the length and time scales upon which the events are to be forecasted.

3) Developing a quantitative understanding of fault interaction and screening at various scales, and of the appearance of cooperative effects on members of fault systems, including earthquake triggering and self-organization phenomena.

4) Examining the roles played by quenched and annealed disorder, as well as the effects of random fluctuations, in determining the dynamical evolution of slip on faults embedded in the elastic and inelastic materials that comprise the earth's crust.

6) Evaluating the applicability of proposed nucleation mechanisms for describing precursory effects, aseismic slip, scaling, and seismicity patterns associated with earthquakes.

8) Elucidating the fundamental equations that describe the time-dependent nonlinear processes associated with earthquakes, fault creep,evolution of fault geometry and stress

transfer in the earth's crust.

Generic GEM computational problems involve:

1) Development of appropriate algorithms and approximations capable of simulations on parallel computers of an interacting fault system down to a 10 - 100 meter length scale.

2) Development of the Problem Solving Environment including advanced visualization techniques that allow both rigorous statistical comparisons with data, as well as more intuitive and visual comparisons. This will use commercial Web technology to support collaboration.

3) Development of a fast multipole or equivalent algorithms for including interactions between fault segments in a computationally efficient manner. These algorithms must allow numerical solution of very large numbers of nonlinear coupled equations using a geographicall distributed, heterogeneous computional environment based on a shared memory configuration. Other novel computational modes may be considered as appropriate.

4) Development of Object-Oriented techniques to allow rapid prototyping with a variety of different approaches and models. These will focus on the use and extension of Los Alamos's POOMA framework.

5) Evaluation of HPCC software methods and systems using benchmarking on a variety of different architectures.

C.9.2 Proposed Tasks and Research Foci by Year

Category I Tasks: Science: Earth Science & Statistical Physics Tasks

1) Collecting and setting up (in object oriented code) the quasistatic Green's functions (linear elastic, viscoelastic, and to a lesser extent poroelastic), so that spatial (and later temporal) multipole representations of them can be established

2) Using the existing data base to establish the basic model parameters, including major fault geometries, slip rates, earthquake histories, and seismic and aseismic activity on the faults. An important decision will be the number and nature of the blind faults (no surface expression) to be included.

3) Establishing the basic specifications for Geographic Information System-type overlays of simulation output with data

4) Analyzing and predicting the physical processes to be expected from the various meanfield versions of the model geometry, including screening and frustrated ordering effects

5) Consideration of the characteristics of the noise to which the faults are subjected due to unmodeled small scale faults, for example: a) geometric (representation of major faults can only be correct down to some scale), b) temporal (smaller faults can easily be created as new fractures in the medium, or unmodeled time-dependent processes in the crust), c) loading (space and time variations in the effective loading rate of the crust, whether real or induced). The nature of the noise, and whether it is correlated or uncorrelated, will have a major impact on the dynamics of the problem.

Year 2: Major Activities:

1) Testing and initial validation of simulated fault systems, by using initial comparisons with field data. In particular, evaluation of adaptive multigrid techniques for implementation of time evolution of failing cascades of fault segments during earthquakes, and initial comparison to source time functions of real events.

2) Evaluation and comparison of the effects expected from the basic friction laws discussed in section C.5, emphasizing in particular the Traveling Density Wave, Cellular Automaton, and Rate and State friction models.

3) Analytical and numerical analyses of excitatory or inhibitory interaction of fault segments, and how these interactions might lead to enhancement of seisimicity or screening effects, as well as the tendency to eventually modify the geometry of the fault system by means of fresh rock fractures

4) Pattern evaluation and analysis for simple fault systems using phase space reconstruction techniques, together with graph analysis for failure (avalanche) trees. Application of scaling and renomalization group analysis to random field approximation of fault loading and interaction.

5) Evaluation of effects of non-planar fault geometries, in particular, how to modify Green's functions, presumably using an adaptive multigrid approach, to account for fractal nature of fault surfaces, and analysis of effects of non planar structures on seismicity patterns.

1) Evaluate requirements for visualization of simulation results, emphasizing in particular Geographic Information Systems-type overlays of simulation results with observational data

2) Develop protocols for calibration and validation of full-up simulation capability, in terms of numerical benchmarking, scaling properties of models, and field data needed to evaluate models

3) Organize and sort field data used for validation, including seismicity, surface deformation from Southern California Integrated GPS Network and Synthetic Aperature Radar Interferometry, source process data fromo TERRASCOPE and strong motion instruments, paleoseismic observations of regional earthquakes, and other data including water level data, in situ stress, and so forth.

4) On the basis of simulations, discuss, evaluate, and develop new hypotheses to be tested with new data collection activities. Data collection to be undertaken by scientists at CalTech, Southern California Earthquake Center, UCSD, Cornell, MIT and Colorado, under other proposals to federal agencies. In other words, develop and promote collaborative modes of activities and scientific interactions as envisioned in GEM approach, focusing particularly on SCEC scientists.

5) Define requirements for future levels and further refinements of GEM simulation technology

6) Transfer technology to third party users as appropriate through WWW interfaces, and publications.

Category II Tasks: Computational Science

1) Develop a flexible object-oriented environment based on POOMA which will support current and expected GEM needs

2) Analyse the specification of GEM Information (geometry and degrees of freedom of faults) so that we can design appropriate data structures needed in databases and software for GEM.

3) Examine the base visualization and I/O requirements of GEM as well as collaboration support needed. Use this to design and develop prototype Problem Solving Environment based on Web Technologies

4) Develop prototype integrated systems for the better understood and simpler GEM computations such as the Cellular Automata approach.

5) Develop Parallel versions of the Statistical Mechanics and Cellular Automata codes and examine the difficult load balancing issues

6) Study and prototype the software needed to support the fast multipole method with the changes needed by GEM compared to better understood gravitational problems.

Many of the activities, including refinement of POOMA, will continue from the previous year. We also expect continued rapid change in available commercial Web Technologies and will incorporate them as appropriate.

1) Develop and use in simulations a simple brute force O(N2) Travelling Wave solution system with fixed (in time) spatial resolution (which can however vary with position in fault system) and fixed position independent time steps. Extend this to variable time and space resolutions and to the well known O(N2) parallel algorithms.

2)Test the first parallel multipole schemes with benchmarking on a variety of machines

3) Incorporate in multipole solver on an ongoing basis, the technologies (friction laws, multiresolution time steps etc.) tested in the simpler simulation systems. Continue to interact with POOMA to ensure that the same modules can be re-used in a variety of simulation approaches.

4) Integrate the simpler simulation approaches into a full operational Problem Solving Environment supporting distributed simulation, data analysis and visualization.

1) Develop Operational Parallel Fast Multipole System in terms of a full capabilities of Problem Solving Environment.

2) Investigate and produce prototypes of a full time dependent multipole method

3) Examine implications of new architecture developments in HPCC

John Rundle, University of Colorado: As Co-Principal Investigator, Rundle will work with Fox to lead, manage and coordinate the scientific activities described in this proposal. Both Rundle and Fox will collaborate closely with all of the investigators to develop and exploit the GEM simulation technology. As part of the management and coordination acitivities, a major workshop will be organized each year, to be held at a centrally located site in Boulder, Colorado, Santa Fe, NM, or Los Alamos, NM, which all of the investigators will attend to discuss and coordinate results, and to plan future work. Additionally, this will be an open meeting that any interested scientist may attend, and previous experience indicates that it may turn into a national workshop on earthquake modeling. There will also be smaller workshops in the New England area, and in California to strengthen coordination and collaborations among the investigators in those areas. Rundle, Jordan, and Minster will also coordinate activities with the Southern California Earthquake Center, and Minster, Rundle and Turcotte will coordinate activites with NASA efforts directed at developing and acquiring Interferometric Synthetic Aperature Radar data for the Southern California area. Rundle is currently chair of the Advisory Council of SCEC, Jordan is a Council member, and Minster is Vice Chair of the Board of Directors of SCEC. Minster, Rundle and Turcotte are heavily involved in team activities for planning an interferomentric SAR mission in the near future.

Year 1 Activities: Rundle's group will collaborate on parallelizing Green's functions, based on standard elastic (e.g., Okada, 1992), viscoelastic (Rundle, 1982a, 1988) and porous elastic Green's functions (Biot, 1962; Bowen and Wiese, 1969; Rice and Cleary, 1976; Cleary, 1977; Rundle, 1982b; and others). Rundle will also participate in load sharing calculations, and will work with Fox on multipole algorithms, by e.g., supplying multipolar terms from various Green's functions. In collaboration with Klein, Rundle's group will focus a major part of their activity on the evaluation of screening effects, in particular, how these will modify the multipole expansions used to propagate stress from on fault to another.

Year 2 Activities: Continue work on above, as well as collaborate extensively with Fox and others to parallelize 1) Traveling Density Wave; 2) Cellular Automaton; and 3) Rate and State friction laws. For all of the friction laws, this task will entail a careful analysis of the mechanics of implied space-time correlation of events, and the patterns that cascades of failing sites display during an earthquake, so that appropriate multigridding techniques can be developed. The primary consideration will be to devise means to regrid areas of faults as needed, to preserve as much detail of the space-time dynamics as possible.

Year 3 Activities: The development of Calibration/Validation & Visualization Software will be the major focus of activities in year 3. Rundle's group, in collaboration with Stein, Kanamori, Minster, Jordan and Turcotte, as well as Fox and the ACL computational group, will use data plus calculations, to design and implement criteria for the construction of visualization techniques. Rundle will coordinate the group's activities to carry out systematic comparisons to 1) statistics of seismicity; 2) surface deformation data from both Southern California Integrated GPS Network (SCIGN) and Interferometric Synthetic Aperature Radar (ISAR); 3) source process observations from both strong motion recordings and TERRASCOPE data as available; 4) paleoseismic observations of earthquakes on regional faults; as well as 5) other data such as water well levels (fluid pressure levels), in situ stress data, creepmeters, strainmeters, GPS campaigns, and so forth.

The requested funding for these activities will include a "geophysical postdoc", a "computer postdoc", a student, and and administrative assistant for help in organizational and communication issues. Emphasis will be on postdoctoral level workers who can help contribute to the success of GEM on the short time scales required.

Geoffrey Fox, NPAC Syracuse University: As Co-Principal Investigator, Fox will work with Rundle to coordinate the group activities with especial responsibity for the computational science tasks. This will draw on Fox's extensive experience with such interdisciplinary work both at Syracuse and Caltech.

Year 1: A critical task will be to understand how the multi resolution multipole approach should be integrated into the GEM computational environment. Fox's group will work with both Rundle and the Los Alamos/Caltech group who will be largely responsible for base software. We will also be responsible for systems integration of the various tools for different models and the development using Web Technologies of a customized collaborative system. We will not develop any base Web technologies but use existing commercial database and collaboration capabilities with which we have extensive experience from other NPAC projects.

Year 2: The focus will be on implementation and integration of the full system with multipole and cellular automata capability. We will also study implications of time dependent multipole methods from both software and algorithmic point of view.

Year 3: The focus will involve substantial use of the full GEM analysis, and the system integration issues will involve interacting with all investigators including the data and simulation visualization. The requested staffing includes some funding for Fox; a postdoc and student who will work on the computational science with software and parallel algorithm experience. These activities will include appropriate benchmarking and MPP architecture evaluations so that we can draw lessons on MPP Systems from the GEM application. There is funding for technical support (.33 FTE) for supporting the Web aspects of the collaborative environment.

Advanced Computing Laboratory, Los Alamos National Laboratory (John V. W. Reynders, Timothy J. Williams, Steve R. Karmesin, and Julian C. Cummings, Michael Warren) The work to be carried out by this group divides into three distinct areas. These efforts will continue throughout the three year period of the proposal.

Object-Oriented Framework Development Activities: We propose to design and implement an object-oriented Seismic framework of layered class libraries for high-performance predictive seismic simulation. This effort will build upon the Parallel Object-Oriented Methods and Applications (POOMA) framework at the Advanced Computing Laboratory.

Fast-Multipole Implementation Activities: The ACL group will focus a large share of their activies on implementing the Fast Multipole computational approach within the POOMA framework. Warren in particular will assis