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 assist in the analysis and development of time and space adaptive fast multipole algorithms for GEM. Initial analysis will include the determination of multipole error bounds for the particular Green's functions involved in the interactions, based on the techniques developed in Salmon and Warren (Salmon and Warren, 1994). Software development will primarily focus on extensions to the fast multipole method to accomodate features which are particular to GEM, based upon the work described in Warren and Salmon (Warren and Salmon, 1995). Techniques to allow an efficient parallel implementation of a multipole method allowing variable time steps will be developed, Additionally, advanced load-balancing techniques to accommodate the varying time and space resolutions in the problem will be investigated. This work will begin in the first funded year of the proposal and will continue throughout the project.

I/O, Visualization, and Data Analysis Activities: In addition to providing a high-level I/O interface, our object-oriented framework will provide mechanisms for run-time visualization of simulation data. Visualization specialists at the ACL are interested in examining the Seismic framework for applications of their parallel rendering algorithms. We propose to work with them through this avenue on runtime visualization of simulation data in the Seismic framework, as well as through the avenue of portable output files for postprocessing; we have already made progress with the POOMA framework on postprocessing using our current parallel I/O implementation.

Both for postprocessing (via I/O mechanisms) and runtime data analysis, aside from visualization, the framework approach has merits we propose to exploit. We propose to develop analysis methods, both specific for the seismic simulation data and general for a wider class of scientific data. In the same way that the framework's high-level programming interface maps well onto solving equations numerically, it maps well onto data analysis and reduction algorithms such as spectral transforms and spatial averages. This work will leverage off other POOMA developments, such as the fast Fourier transforms developed for particle-in-cell simulations.

James Crutchfield, Santa Fe Institute: Crutchfield will collaborate with Klein, Rundle and Fox to extend recent results from dynamical systems theory to the GEM simulations. Crutchfield's new approach to computation in dynamical systems --- computational mechanics (Crutchfield, 1992; Crutchfield and Young, 1988) --- suggests several novel methods to analyze and visualize earthquake dynamics. Crutchfield's work will focus on both computational methods and algorithm development for simulation and analysis of earthquake modeling and prediction. One goal is to develop a range of visualization methods that parallel and complement the simulation and analysis methods.

In the analysis of nonlinear processes work to date has focused on statistics such as fractal dimension, metric entropy, and related scaling properties. Crutchfield recently introduced a new, but complementary observable---the statistical complexity---which is a measure of the amount of memory used by a process to produce its behavior.(Crutchfield, 1992; Crutchfield and Young, 1988) It should be quite helpful in characterizing new aspects of earthquake dynamics and its prediction. The estimation of this new observable, and several allied statistics, will be provided by the appropriate software libraries.

Crutchfield will also adapt epsilon-machine reconstruction, an automated modeling method that can reconstruct continuous-state difference/differential equation models from data series. The goal is to implement several methods for epsilon-machine reconstruction with an eye toward allowing users to select which is most appropriate for the given field data or simulation data. Obviously, these kinds of data have different characteristics. For example, data from simulations can be much more extensive, assuming that one invests the computation effort to produce it. This is typically a benefit, since statistical estimates can be greatly improved. But making use of such data sets brings up issues of optimally using large amounts of data.(Crutchfield, 1992; Crutchfield and Young, 1988; Crutchfield and McNamara, 1987; Crutchfield and Hanson, 1993)

There are several novel methods that Crutchfield and Hanson have developed, based on their computational mechanics analysis of spatial systems, that will be adapted to visualize the dynamics of earthquake models.(Crutchfield and Hanson, 1993; Hanson and Crutchfield, 1992) There are three goals in this.

The first is to adapt the nonlinear spatial filtering methods they've developed to track localized disturbances. These methods allow for a great reduction in the complexity of visualizing the temporal development of patterns in one and two spatial dimensions.

The second goal is to develop software tools that will allow analysts to build up the attractor-basin portrait for the earthquake models. This will give a way to summarize the spatio-temporal dynamics in terms of their qualitative dynamics. The GEM model, being spatially-extended collection of relaxation oscillators, poses interesting problems to conventional dynamical systems modeling methods. Several of these have been addressed by Crutchfield's recent work on the attractor-basin portraits of spatial processes.(Crutchfield, 1988; Crutchfield and Kaneko, 1988) The question now is how this new perspective can elucidate the underlying mechanics particular to the GEM model.

The third goal is to merge wavelet analysis methods, which are appropriate for localized space-time events such as earthquakes, with computational mechanics. The hope is that the combined method will allow for a function basis that is not susceptible, as are Fourier methods, to the artifacts of multiple time-scale bursting dynamics. Access to the new method will be provided via software and graphical display libraries that are integrated into the overall simulation package.

These tasks will begin in year 1 of the proposed work and continue on through temination of the project. In particular, simulations developed in years 2 and 3 will be used as they become available.

Roscoe Giles, Boston University: In Year 1 of the proposal, Giles will work with Klein and Rundle to develop and extend theoretical models for earthquakes, and upon analysis methods. His particular emphasis will be on the computational realization of these models on high performance parallel computers, including the distributed shared memory architectures available at Boston University's Center for Computational Science. An important goal is to be able to begin preliminary evaluations of methods to validate the models using the various southern California databases. Giles will work with York, Klein, Stein, Kanamori, Minster, Jordan and Rundle to comparing model results to the SCEC and USGS databases. Giles will also begin development of a common WWW/JAVA user interface for internal dissemination of all our work within the collaboration.

In Year 2, Giles will work with Klein and York on incorporating additional features such as surface roughness and noise into the models. He will work toward the integration of statistical and CA simulations into a common object oriented framework. He will also examine the scalability of 3D models.

In Year 3, Giles will continue the activities described above, and will in addition focus on refining the techniuqes needed for comparison with observational data and for refinement of the model and WWW interface.

Thomas Jordan, Massachusetts Institute of Technology: Jordan's research will focus on the phenomena of "slow" earthquakes and their implications for earthquake nucleation and fault friction in the context of GEM simulators. In this activity, he will collaborate in particular with Kanamori and Rundle. Slow earthquakes -those events whose characteristic velocities are small compared to ordinary events- have been observed in a variety of tectonic environments, including southern California (Kanamori, 1989; Kanamori and Hauksson, 1992) and they are very common on the oceanic ridge-transform system (Kanamori and Stewart, 1974; Okal and Stewart, 1982; IhmlŽ and Jordan, 1994). Slow events are especially interesting because they appear to precede and perhaps initiate some large earthquakes (Kanamori and Cipar, 1974); Cifuentes and Silver, 1989; IhmlŽ et al., 1993; IhmlŽ and Jordan, 1994). The questions to be addressed are as follows: (1) Can slow earthquakes be generated on faults governed by laboratory-derived friction laws or by other friction laws used in GEM? (2) If so, over what range of parameters is this behavior observed? (3) If not, what other types of instabilities might be capable of generating such ruptures and how could they be incorporated in GEM? (4) What role do slow events play in the earthquake initiation process? Jordan's work will tie in with that of Marone and other investigators to assess the role of friction laws in GEM.

Hiroo Kanamori, California Institute of Technology: Kanamori will work with Rundle, Jordan, Turcotte and Marone to investigate crustal processes leading up to fault instability. In particular, he will investigate processes that involve fluids, such as rectified diffusion. Triggering of earthquakes by seismic waves generated by earthqukes from different segments of fault will play an important role in controlling the temporal and space variations of seismicity. In particular, strong shaking caused by a major earthquake will influence seismic activities in the adjacent areas. He will investigate several examples of such triggered seismicity, and interpret them in terms of the numerical simulations to determine the existence of seismicity patterns.

William Klein, Boston University: In the first year of the grant Klein will, in conjunction with Rundle and Giles, be involved with the theoretical analysis of fault interactions as described in section ?? of the proposal. This will involve applications of renormalization group techniques as well as ideas such as Supersymmetry to the long range interaction slider block models as well as the Traveling Density Wave model. In addition, Klein will be involved with York, Rundle, Giles, Reynders, Warren and Krogh in building models of fault interaction based on an analysis of data from SCEC databases.

In the second year Klein will continue with the theoretical effort while also working with Giles and York to investigate the effect of adding noise, the way plates are updated and surface roughness to the model. These investigations will be done in parallel with a theoretical investigation of the same phenomena.

In the third year of the grant Klein will be involved with the evaluation of the model data and its comparison with measurements of similar parameters on real faults. In addition he will be part of a subgroup whose task is to use the results of the theoretical and modeling effort to suggest new measurements for the seismologists and experiments for laboratory scientists.

Karen E. Carter Krogh , Los Alamos National Laboratory: Carter Krogh will be involved in the initialization of the proposed earthquake models by incorporating actual field data which encompass the morphology, geometry and kinematics of active faults in the study area. Specifically, a three-dimensional model of seismically significant faults will be generated in collaboration with Rundle, Kanamori, Turcotte and Jordan by incorporating both observational and derived data. Data will include both surface and subsurface faults. Whereas surface faults are visible and measurable, the character of subsurface faults will be modeled through the use of geophysical data. Working with Rundle and Minster, Krogh will correlate information such as seismic, GPS and field data, and match it to fault surfaces. Once the computational model has been constructed, simulations on that model will be constrained, in part, by local and regional stresses. Local stresses are influenced by the spatial distribution and geometries of seismically significant faults. Additionally, the geometry and morphology of secondary and minor faults will be analyzed as indications for local stress states. Geologic methods will be developed to enable the incorporation of such unmodeled minor faults as background "noise" for the major faults in the simulations. Regional stresses, that are invariant over the proposed crustal scale for short time steps, have been determined by previous studies and will be included. Similarly, in collaboration with Minster and Jordan, Krogh will document kinematic information on active faults and will used it to constrain the tectonic framework within which the computational earthquake models will function. Results from the model will be validated by comparing simulated data to observational data. For example, new surface ruptures over the past decade can be compared with those simulated, for that time step in this model. This work will begin in the first funded year of the proposal and will continue through the analysis and validation phases of the project.

Chris Marone, Massachusetts Institute of Technology: Marone will work to incorporate and test laboratory-based friction laws within GEM and to study the role of friction in simulations of earthquake physics. The work will fall into two main areas: 1) incorporation of existing rock friction laws into computational models of earthquakes, and 2) investigation of the problem of scaling laboratory results to earthquake faults and fault systems. Traditionally, cellular automata and spring-block models of earthquakes have used extremely simplifed friction laws. This is justified by the computational efficiency of using such laws. However, in some cases these laws omit critical features of rock friction. Thus, an important question is, what effect does this have on model output, and what other parameters of the simulation are impacted by the use of a given friction law? Because of the relative simplicity of the current generation of GEM simulators it is important to begin testing realistic friction laws and calibrating the effects of simple friction laws. In particular, Marone will work with Minster and Shkoller to better connect GEM results with laboratory-derived friction laws and related earthquake modeling studies. Marone will also study the scaling problem via laboratory experiments and numerical studies. The most important scaling parameter is the critical slip distance (Marone and Kilgore, 1993), which enters in either slip-weakening or rate and state friction laws. Marone will use laboratory-based scaling models to investigate earthquake nucleation (Roy and Marone, 1996) and will collaborate with Jordan to study the phenomena of slow earthquakes in the context of rate and state friction laws and GEM simulators.

J.-Bernard Minster, University of California, San Diego: Minster, together with a postdoctoral fellow (S. Shkoller), will work to formulate rate and state dependent friction laws in the form of partial differential equations that are compatible with constraints imposed by a variety of theorems relating to the dynamics of nonlinear dissipative systems. Recently, they have shown that the dynamical system describing 1D frictional sliding, using the generally accepted class of state and rate-dependent friction laws, displays deterministic chaos with attracting sets of fractal dimension. This result is largely due to the recent discovery that hyperbolic attractors persist in a finite dimensional phase space under the small perturbation of numerical time integration (Pliss and Sell, 1991). Since a physically-based law for sliding of the earth's faults will require a complete 3D model described by a dissipative partial differential equation in infinite dimensional phase space, it is of the utmost importance to generalize the perturbation result to overflowing invariant manifolds contained in this infinite- rather than finite-dimensional setting. Minster and Shkoller will thus work with Marone from MIT to reformulate these friction laws for GEM simulations, making use of constraints imposed both by laboratory data and theory. In addition to these activities, Minster will also work on merging current and anticipated future SAR Interferometry data with the GEM simulations as described in section C.7 on data validation issues.

CACR Group (John Salmon & Paul Messina), California Institute of Technology: As part of the GEM project, John Salmon will work with Earth Scientists (Rundle, Kanamori, Jordan, Minster, Turcotte) to implement parallel "Fast" algorithms for the evaluation of stress fields from fault displacements. Initially, in the time-independent formulation of the Traveling Density Wave Model, the stress field induced by a fault is represented by a Green's function with a long-range component, i.e., a field for which it is not mathematically permissible to neglect interactions beyond some radius. Salmon's previous work in astrophysics and vortex dynamics uses "Fast" summation methods and oct-tree data structures to evaluate similar Green's functions for Newtonian gravity and the Biot-Savart law. Warren and Salmon have developed a library which separates the "physics" from the "computer science", and which allows one to alter the interactions, error tolerances, etc. without worrying about the underlying parallel computational mechanisms (Warren and Salmon, 1995). Salmon's primary task will be to work with the Earth Scientists to allow their codes to use this library. This will have two important consequences: the codes will run in O(N log N) time or O(N) time, and they will run efficiently on a wide variety of shared and distributed memory parallel computers. Further refinement of the physical model will require introduction of spatial and/or temporal dependences into the Green's function, modeling non-uniform physical characteristics of the earth and the radiative time-dependence of the stress field, respectively. Here it is likely that the object-oriented framework supplied by POOMA will be an important asset. It will then be necessary to either create an interface between POOMA and the tree library, or to reimplement some of the "Fast" methods supplied by the tree library in POOMA itself. Salmon's role will be to work closely with both the POOMA team and the Earth Scientists to find the most efficient methods of implementing these advanced techniques.

Ross Stein, United States Geological Survey: Large spatial datasets provide critical input for testing models. Such data includes three-dimensional displacement vectors for hundreds of points on the earth surface (Global Positioning System and Very Long Baseline Interferometry); baseline changes between a radar satellite and millions of points spaced 30 m apart on the earth's surface (Synthetic Aperture Radar); and positions, geometries, and fault-slip information on the several thousand active faults and folds longer than 8 km in the western U.S. Such data must be easily interrogated and changed as errors are uncovered and modifications are made. Model output includes calculations of stress and stress change, fault friction, fault slip and earth displacement at the earth's surface and also at depth. Output datasets will exhibit a strong time-dependence and thus will be frequently updated, compounding the handling problems. These very large spatial and temporal datasets must nevertheless be accessible to develop physical insight. The visualization needs of this project will thus overwhelm existing desktop-type Geographic Information System (GIS) codes such as ArcInfo or MapInfo, and matrix manipulation and display codes such as Transform and Dicer. In addition, remote collaboration by the many investigators on this project means that the visualization tools must be client-server based.

Stein will work with Reynders, Williams, and the LANL/ACL group, and in addition with Giles and York to create a database with all known active faults longer than 8 km in the western U.S. Distributed deformation at fault bends, break, and ends will be implemented as well to avoid singularities and artificial instabilities. Stein will also work with Rundle and Minster to develop an inventory of historical earthquake sources with slip distributions (where known), slip histories (where available), and the net slip distribution. Procedures to download and display continuous and survey-mode GPS and other surface displacement data will then be incorporated. Output datasets must be converted to 3D animations in which one can control and adjust the time, and three-dimensional spatial range of the output easily for inspection.

Donald Turcotte, Cornell University: Turcotte will focus on two areas in collaboration with Rundle: 1) relating the stresses driving slip on the simulated fault systems to the underlying mantle dynamics, and more extensively, 2) validating simulations with seismicity data and Interferometric Synthetic Aperature Radar data. With respect to the first activity, Turcotte will draw on his long association with the mantle convection and dynamics community to assess the nature and patterns of crustal stress in the simulations, and estimate the extent to which these accurately represent the magnitude and orientation of stresses arising from mantle flow.

With respect to the second area, a primary data set is the spatial and temporal distribution of seismicity. High resolution data sets are clearly preferred and for this reason, comparisons will emphasize regions where many earthquakes occur and where seismic networks archive quality data over a relatively wide range of magnitudes. For these reasons, aftershock sequences are particularly attractive for study. Seismic networks in northern and southern California also give excellent quality data with particularly high resolution data available for the Parkfield section of the San Andreas fault in Central California. Although these areas have the highest priority, it will also be important to compare patterns of seismicity there with those observed in other areas of the United States, including the Pacific Northwest, the Wasatch Fault, the New Madrid seismic zone, and seismic zones in the eastern United States.

While seismic data are a primary data source, the great improvements in geodetic observations are providing complementary data. GPS data have provided sub-centimeter resolution on occupied geodetic monuments for five to ten years. However, the density of monuments is still too sparse to provide the desired resolution. Recent studies using radar interferometry have yielded sub-centimeter accuracy with a 50 meter (pixel) resolution in many regions. With monthly or even weekly repeat passes, this data can provide high resolution, extremely accurate spatial coverage of stress and strain release, coseismically, aseismically, or with continuous and discontinuous creep. In particular, Turcotte, Rundle and Minster are on the NASA ESSP pathfinder SAR team, and they will plan for the integration of that data into GEM models in the future.

Bryant York, Northeastern University: In the first year York will focus on the development of parallel/distributed algorithms for the extension to three dimensions of cellular automata (CA) based earthquake models, such as the well-known "slider block" model. Initially he will consider the problem of describing 3D fault geometries for multiple faults in a manner which allows them to be easily distributed and computed with minimum communications delays. The two key issues here are: (1) how to define and implement variable-resolution gridding, and (2) whether or not it is feasible to treat the 3D model as a network of 2D slider-block models in which distant blocks are connected through an elastic medium. In order to support this work, York will begin a careful analysis of data from several SCEC databases in conjunction with Rundle and others. York will be on sabbatical leave during academic year 1997 - 1998 and plans to spend that leave at the Boston University Center for Computational Science. This sabbatical will afford York substantial time to work closely with Giles, Reynders, Williams, Warren, Salmon, Klein, Fox and Rundle on several aspects of the GEM project which complement the work on the development of object-oriented modeling environments for earthquake physics put forth in the earlier proposal.

In the second year York will work with Fox, Giles and Klein on the addition of several novel features to CA-based earthquake models. In particular, we will look closely at the incorporation of surface roughness and fractal dimension parameters as well as the systematic incorporation of noise. We will look into the "plate update" scheme and investigate model dependence on the order of block evaluation. We will be particularly concerned with parallel block update, correlated block movement, and continuous deformation models for connected blocks. A careful sensitivity analysis will be done of all the parameters and we will begin to analyze the scalability characteristics of the enhanced 3D models.

In the third year York will work with Giles, Klein and Stein to incorporate time renormalization in a generic way into the GEM modeling system and to develop tools for evaluating model performance against recorded data.

Salary support for York is requested for 2 months of summer each year and 1/2 salary for academic year 1997-98 while York is on sabbatical. The sabbatical salary will be spread across two grant years. No graduate student support is requested as York has a fellowship student working with him. Request for undergraduate programming support was included in the earlier proposal.

Simulations will be carried out in a variety of computational environments. These are described below.

The computing hardware in the LANL Advanced Computing Laboratory is about to begin a multi-year, phased upgrade. We will focus on clustered SMP technology. In the broadest terms, one might describe an initial SMP-cluster arriving in Fall, 1996, of about the same capacity as our 1024 node TMC CM-5 but with greater efficiencies; a 5-fold technology refresh in Fall, 1997, and a final configuration in Winter, 1998.

October 1996 October 1997 October 1998

---------------- ---------------- ----------------

100 Gflops 600 Gflops 1000 Gflops

30 Gbytes memory 150 Gbytes memory 250 Gbytes memory

600 Gbytes disk 3000 Gbytes disk 5000 Gbytes disk

multiple HIPPIs multiple HIPPIs multiple HIPPI 6400s

The visualization environment for the Grand Challenges will consist of state of the art high-end SGI servers for hardware assisted visualization. At the start of the program (Oct 1, 1996), our graphics server facilities will consist of a dual headed Onyx with 2 RE-2 graphics engines, 8 CPUS, 24GB/disk and a 2 processor R10K with an Infinite Reality graphics engine. We expect to phase out the Onyx and bring on-line a minimum of 3 Infinite Reality engines over time. These will be connected via striped HIPPI with the SMP-cluster forming the main compute capabilities for the ACL. In addition to the high-end graphics hardware, the ACL provides video laser-disk assisted animation capabilities for generation of video tapes (SVHS/VHS/Hi-8/8). We encourage participants to leverage the visualization environment in the ACL with mid-range SGI desktops such as the Maximum Impact Indigo^2 with a full complement of memory.

The Northeast Parallel Architectures Center (NPAC) is a research and development center focusing on high performance computing and communications. Established at Syracuse University in 1987, the center is directed by Geoffrey C. Fox, a pioneer in the development and application of parallel computers, who assumed NPAC leadership in 1990. NPAC has a strong emphasis on computational science applications where one seeks to obtain "real solutions to real problems." NPAC's InfoMall technology transfer program puts high-performance computing and communications (HPCC) to work in industry. Other major projects include research and development in the areas of parallel languages and parallel compilers, including work on High Performance Fortran (HPF) -- a standardized parallel version of Fortran; distributed computing; parallel database technology; integration of relational databases with the NII; parallel algorithms; distributed and cluster computing; and networked digital multimedia applications providing Information, Video, Imagery, and Simulation on demand.

NPAC's role in this collaboration will focus on appropriate use of the best HPCC technologies both the Southern California models, as well as in GEM, such as MPI and HPF. This will naturally use the expertise of the NSF Science and Technology center CRPC -- the Center for Research in Parallel Computation -- led by Kennedy at Rice where Fox is a founding member. The Southern Calfornia model and GEM are both natural projects in which to exploit World Wide Web technology where NPAC is a pioneer for expoiting both databases (http://asknpac.npac.syr.edu) and distributed computing (http://www.npac.syr.edu/factoring.html). The Web naturally links the data and simulation aspects of GEM and this integration is at the heart of a joint MARINER-NPAC project WebWork where Fox is collaborating with Marina Chen -- the head of the computer science department at Boston University.

MARINER:

The MARINER (Metacenter-Affiliated Resource in the New England Region) project (established at Boston University in October 1995) provides a focal point for HPCC/NII related activities in New England. As a Metacenter Regional Affiliate, its goal is to act as a catalyst for wide diffusion of technologies and as a channel to the Metacenter. MARINER will draw on the University's advanced computing and communications facilities and build on its broad experience in computational science, scalable parallel computing and wide area networking in order to provide education and training to a wider audience as well as to facilitate computational science research and education.

Computational facilities that will be available for this project are the supercomputing facilities associated with the Center for Computational Science at Boston University. At present these facilities include a 38 processor SGI POWER CHALLENGEarray with a peak performance of 13 Gflops, high performance networks and associated visualization peripherals, and a laboratory of networked SGI graphics workstations. The POWER CHALLENGE array represents the first phase equipment of Boston University's third generation supercomputing facility. The equipment for the second phase will consist of SGI's next-generation, highly-scalable supercomputing system with about 25 Gflops of peak performance. The first phase systems will continue in operation, so that our total computational power will increase to about 38 Gflops at that time. About 12 months later, the third and final phase will consist of an upgrade to the second phase system, bringing it to a peak power of about 65 Gflops, or a total of more than 75 Gflops including the first phase systems. The computing environment also includes numerous high performance workstations, HiPPI and FDDI high performance networks in addition to universal Ethernet connectivity and a 10Mb per second connection to the NSF backbone. In the near future we plan to provide ATM connectivity by routing through the SGI Power Challenge machine to a Fore Systems ATM switch. ATM will also be brought to several Indigo-2 workstations serving scientific computing at the University.

G.1 Summary of Past and Present Collaborations

The challenges posed by the GEM project are at least equal to those found by atmospheric scientists working on GCM models, because the physical processes are less susceptible to direct observation. We have therefore brought together an extremely high quality team of some of the best theoretical and observational(experimental) geophysicists, condensed matter physicists and computational scientists in the country. These scientists have all collaborated in the past in various subsets.

Rundle and Klein have collaborated since 1987 on a series of papers whose subject is the physics of earthquakes, and the relation between slip events and nucleation phenomena. The most recent of these papers was published in Phys. Rev. Lett., 75, 1658-1661, 1995. Rundle has also collaborated with Turcotte, Kanamori, Minster, and Stein at various times, and has published with all of them on the subject of stress redistribution in earthquakes; source processes of earthquakes; and geodetic tests of earthquake models. Crutchfield and Rundle are currently developing pattern recognition algorithms for seismicity distributions. Jordan and Minster have a long history of collaboration on plate tectonic and other research dating back 25 years to the time they were graduate students together at Caltech. Turcotte and Allegre are long time collaborators in reseach dealing with mantle convection and the applicability of Renormalization Group techniques to scaling problems in the earth sciences.

Fox has worked in many interdisciplinary projects during his years at Caltech, including neural network activities with Hopfield (unfunded collaborator with Rundle), as well as several other data analysis and simulation projects with geophysicists within the Seismological Laboratory. These are described in a recent book "Parallel Computing Works!" Fox is also collaborating with Marina Chen, head of the computer science department at Boston University. Fox has successfully integrated HPCC technologies into many applications. The book Parallel Computing Works describes some 50 examples from his time at Caltech where he collaborated with the Geophysics group there on seismic (Clayton) and mantle (Hager before he left) simulations on parallel computers.

Fox has also collaborated with the computer and computational scientists at the Advanced Computing Laboratory at Los Alamos National Laboratory. In particular, two of his students (Warren & Salmon) are currently funded by ACL, either directly (Warren) or indirectly (Salmon). These former students played major roles in developing the fast multipole algorithms. Rundle has known Krogh for over a year. Krogh will work closely with the ACL group and will act as an in-house interface between the earth science community and the ACL lab.

York, Klein and Giles have all collaborated together in the past in carrying large scale supercomputer simulations of nonlinear problems. York and his students are curently collaborating with Rundle to port Rundle's 1988 model to the SGI PCA in Fortran 77. In addition, they are implementing primitive surface patch models for use in developing sophisticated fault geometries.

G.2 Investigator Vitae

1. John B. Rundle, University of Colorado

Computer Equipment Justification:

As a part of the University of Colorado budget, the Principal Investigator is requesting funds to purchase a fully configured Maximum Impact Indigo2 SGI desktop workstation, together with two P6 Intel machines each with 256Mb memory and two 200 MHz processors. We also plan to purchase a 10Gb disk for each workstation, a tape drive backup for the Indigo, a CD ROM reader, and a new laser printer. The rationale for the SGI Indigo is that no such graphics engine exists here at CIRES, and there is a critical need for visualization facilities that will be compatible with ACL at Los Alamos, NPAC at Syracuse, MARINER at Boston University, and so forth. The two 200 MHz P6 Intel machines will be a substantial upgrade from the SPARC 512 90 MHz machine and the two SPARC 5 90 MHz machines currently available in the Geophysical Simulation Laboratory at CIRES. The P6 machines will be used for code development and validation using LINUX software, as well as for data display and simple visualization (as compared to the capabilities of the SGI Indigo). A state-of-the-art color dye printer is needed because the current CIRES color dye printer is 5 years old, is very slow in part because it has only a serial connection port, hence a very low bandwidth, and it is heavily overused. Thus the current facility color printer will not accomodate processing of the large volume of color images the GEM project will generate, without making an unacceptably large impact on the current CIRES users.

Justification for Administrative Assistant

The management plan (section C.2.4) for administering the GEM Grand Challenge calls for an executive committee, a full technical committee, and an external advisory board. These committees will meet regularly, in Boulder and elsewhere, to assure that the proposed work is being carried out in a responsible and expeditious manner. In addition, it is anticipated that a series of major annual reports, as well as other technical reports at irregular intervals, will be required by the sponsoring agency. These reports will be in addition to the normal technical papers published in the peer reviewed literature. Finally, long experience, both by the PI and with the organizations such as the Southern California Earthquake Center, clearly shows that the nature of the work (implications for earthquake prediction/forecasting) will attract widespread public interest. There will thus be a considerable need for substantial interaction with the public, California Local & State, and Federal Government agencies. None of these tasks would be eligible for coverage by the usual ICR resources, thus there is a need for an assistant to handle these duties.