Julian J Bommer

Norman A Abrahamson

Fleur O Strasser

Alain Pecker

Pierre-Yves Bard

Hilmar Bungum

Fabrice Cotton

Donat Fäh

Fabio Sabetta

Frank Scherbaum

The challenge of defining upper bounds on earthquake ground motions

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Recent studies to assess very long-term seismic hazard in the United States and in Europe have brought the issue of upper limits on earthquake ground motions into the arena of prob lems requiring attention from the engineering seismological community. Few engineering projects are considered suffi ciently critical to warrant the use of annual frequencies of exceedance so low that ground-motion estimates may become unphysical if limiting factors are not considered, but for nuclear waste repositories, for example, the issue is of great importance. The definition of upper bounds on earth quake ground motions also presents an exciting challenge for researchers in the area of seismic hazard assessment. This paper looks briefly at historical work on maximum values of ground-motion amplitudes before illustrating why this is an important issue for hazard assessments at very long return periods. The paper then discusses the factors that control the extreme values of motion, both in terms of generating higher amplitude bedrock motions and of limiting the values of motion at the ground surface. Possible channels of research that could be explored in the quest to define maximum pos sible ground motions are also discussed.

HISTORICAL BACKGROUND

In the period between the recording of the first strong motion accelerograms in the Long Beach earthquake of March 1933 and the end of the 1960's, a number of studies were published proposing possible upper limits on earth- LGIT, Universite Joseph Fourier, France 7. Institut fiir Geophysik, ETH-Honggerberg, Switzerland 8. Servizio Sismica Nazionale, Italy 9. lnst. fii r Geowissenschaften, Universitat Potsdam, Germany I 0. Studer Engineering, Switzerland quake ground-motion amplitudes. Some studies were purely empirical and influenced to a large extent by the El Centro recording of the 1940 Imperial Valley earthquake: Hausner (1965) proposed that peak ground acceleration (PGA) would not exceed 0.5 g; Newmark (1965) proposed a limit in the range 0.5-0.6 g on PGA and between 76 and 91 cm/s on peak ground velocity (PGV); and Newmark and Hall (1969) proposed a limit of 0.75 g on PGA and agreed with 91 cm/s as the limit on PGV. Newmark and Rosenblueth (1971) referred to the estimate made by Housner (1965) and argued that the upper limit must be higher, at least 1.0 g and possibly 1.5 g. Their argument for this latter value is based on the fact that surface accelerations in the vertical direction exceeding 1.0 g had been inferred from observed effects in many earth quakes, notably the 1897 Assam earthquake; the estimate of 1.5 g for the limit on the horizontal acceleration is then inferred from the rule of thumb that vertical accelerations are generally of the order of two thirds of those in the horizontal direction. A recent study by Anderson Q. G. Anderson, per sonal communication, 2003) using more than 3,000 seismo grams from the Guerrero (Mexico) network found good agreement in general between the distribution of the V/H ratio and a lognormal distribution with mean of 0.67, but that there is a deviation from the lognormal shape in the upper 15 percent of the cumulative distribution function, where more high ratios are in the data than the lognormal dis tribution would predict. Therefore, the use of the 2/3 ratio to infer maximum horizontal acceleration from the maximum vertical acceleration might not be appropriate.

Other studies used simple models of slip on a fault, which were essentially rock mechanics solutions, such as Ambraseys and Hendron (1967), who estimated maximum values of PGV in rock in the range of 90 to 120 cm/s. Ambraseys (1969) later revised the estimate to include an upper limit of 150 cm/s. Hanks and Johnson (1976) subse quently combined a dynamic faulting model with the limit ing strength of rock to estimate a maximum PGA of 0. 75 g based on average rock strength; considering regions of higher stress in areas of greater rock strength, they estimated a more likely upper bound to be 1.8 g. McGarr (1982) per formed similar analyses for inhomogeneous faulting and Julian Bommer 1 , Norman Abrahamson 2 , Fleur Strasser 1 , Alain Pecker 3 , Pierre-Yves Bard 4 , Hilmar Bungum 5 , Fabrice Cotton 6 , Donat Fäh 7 , Fabio Sabetta 8 , Frank Scherbaum 9 , and Jost Studer 10 related the maximum ground motions to the tectonic regime, leading to maximum PGA values of 0.4 g for exten sional regimes, 2.0 g for compressional regimes, and 0.7 g for pure strike-slip.

The studies described in the previous paragraph were all based on the maximum possible strength of radiation from the seismic source. Ambraseys (1970) pointed out that for the horizontal components of motion, nonlinearity and the lim ited shear strength of soil deposits control the maximum accelerations that can be transmitted to the surface, leading to estimates of maximum PGA on normally consolidated clays of 0.10-0.15 g, 0.25-0.35 g for highly plastic deposits, and 0.50-0.60 g for saturated sandy clays and medium dense sands (Ambraseys, 1974). Similar values have been suggested by Dowrick (1987) as summarized in Table 1. This range of values was later supported by both empirical data and theo retical nonlinear models of soil response (Mohammadioun and Peeker, 1984).

Following the San Fernando earthquake in February 1971, which more than doubled the databank of strong motion accelerograms available at the time, attention shifted from consideration of upper bounds to the derivation of empirical curves through regression analysis, although there have been a few excellent studies of extreme ground motions (e.g., Oglesby and Archuleta, 1997). As the database of strong-motion records has continued to grow at ever-increas ing rates, with expanding accelerograph networks throughout the seismically active areas of the world, the number of ground-motion prediction equations has grown in propor tion (e.g., Douglas, 2003). Common to nearly all of these empirical equations, often referred to inappropriately as attenuation relations, is the assumption of a lognormal distri bution of the residuals, resulting in the modeling of the alea tory variability in the ground-motion predictions as a zero mean Gaussian distribution characterized solely by its stan dard deviation, i.e., ground motions are formally considered to be unbounded.

At the same time that empirical equations started to be derived in large numbers during the 1970's, probabilistic seis mic hazard analysis (PSHA) was becoming widely adopted in engineering practice, following the presentation of the funda mental concepts by Cornell (1968). Although not included in the original formulation by Cornell (1968), the contribu tion of the scatter in the ground-motion prediction equations was incorporated into the calculations of annual frequencies of exceedance in the first widely available computer program for performing PSHA (McGuire, 1976). The integration across the lognormal scatter in the ground-motion prediction equations has now become a standard and fundamental ele ment of PSHA (e.g., Bender, 1984). Recent engineering projects have shown the limitations of the lognormal formal ism by extending PSHA to probability levels previously not explored, and thus given rise to the need to return to the issue of upper bounds on earthquake ground motions.

WHY UPPER BOUNDS NEED TO BE DEFINED

Upper limits on earthquake ground motions have recently been identified as the "missing piece" from seismic hazard assessment, for both deterministic and probabilistic approaches [START_REF] Bommer | Deterministic vs. probabilistic seismic hazard assessment: an exaggerated and obstructive dichotomy[END_REF]. Deterministic seismic hazard assessment (DSHA) is often interpreted to define the worst case ground motion. DSHA should therefore be based on the least favorable combination of earthquake source characteris tics and location, and the strongest ground motion that could be generated by this scenario. In practice, DSHA generally uses the logarithmic mean or mean-plus-one-standard-devia tion level of ground motion from predictive equations (e.g., [START_REF] Krinitzsky | How to obtain earthquake ground motions for engineering design[END_REF], which will generally be significantly below the worst-case scenario [START_REF] Bommer | Uncertainty about the uncertainty in seismic haz ard analysis[END_REF]. If DSHA is to be used to define the maximum earthquake loading to which a structure may be subjected, then an estimate of the upper limit on the ground motion that a particular scenario could generate is needed. Brune (1999) observed that PSHA using ground-motion prediction equations with untruncated lognormal scatter may overestimate ground motions with very long return periods. This inference was based on observations of the stability of precariously balanced rocks in the Mojave Desert. The need to define upper limits on ground motions in PSHA becomes clearly apparent only when ground motions are calculated for very low annual frequencies of exceedance. For very long return periods, the hazard estimates are driven by the tails of the untruncated Gaussian distribution of the logarithmic residuals (Anderson and Brune, 1999;[START_REF] Abrahamson | State of the practice of seismic hazard eval uation[END_REF]. The effect of truncating the distribution at different levels above the mean is illustrated in Figure 1. This figure indicates Peak Ground Acceleration (g) Peak Ground Acceleration (g)

Q) 0 c 0.01 � 0.001 Q) Q) 0 X w - 0 >. 0 c Q) ::I C"
Q) 0 c 0.01 � 0.001 Q) Q) 0 � 10-4
.A. that for the situation analyzed therein, at the 1 o-4 level which has been widely used as the basis for seismic safety analyses in nuclear installations in the past-the difference in the resulting hazard between truncating at 3 sigmas or 6 sig mas is almost negligible. Figure 1 also clearly shows that at annual frequencies of exceedance of the order of 1 o-7 or 1 o-S ' the issue of whether the truncation level is at 3, 4, 5, or 6 sig mas becomes a controlling factor on the computed hazard.

The hazard curves shown in the figure are for one particular configuration and the sensitivity to the truncation level may not always be so high, particularly for sites closer to seismic source zones. Moreover, the two sets of curves indicate that sensitivity to the truncation level is also related to the under lying seismicity rates.

For most engineering projects, such long return periods are of no relevance, but for the rare situations where the risk analysis must consider such extreme cases, the definition of upper bounds becomes a necessity. The PSHA carried out for the nuclear waste repository at Yucca Mountain in the Nevada desert (Stepp et al., 2001) considered ground motions for annual frequencies of exceedance as low as 10-8 , and as a result of not using truncations, extremely-and pos-sibly unphysically-high levels of ground motion were com puted. For example, for a hypothetical surface site near Yucca Mountain, the PSHA showed the mean PGV at the 10-7 and 1 o-s values as being of the order of 6.5 m/s and about 13 m/s respectively (L. Reiter, personal communication, 2003). Additional modifications, particularly those related to scaling real earthquake records, increased some of the ground motions even further. For example, some of the scaled ground motions at the 10-7 level reached PGA values as high as 20 g and PGV values of almost 18 m/s [START_REF] Corradini | Letter from chairman of the US Nuclear Waste Technical Review Board to the Director of the Office of Civilian Radioactive Waste Management[END_REF].

The Yucca Mountain PSHA has been one of the first comprehensive applications of the SSHAC Level 4 proce dures (Budnitz et al., 1997) for expert elicitation in seismic hazard assessment. A more recent application has been the PEGASOS project to perform a comprehensive seismic haz ard assessment for nuclear power plants in Switzerland (Abra hamson et al., 2002), which has considered annual frequencies of exceedance as low as 10-7 . In large part due to the outcome of the Yucca Mountain PSHA, the PEGASOS project has required experts in both the ground motion and site response subprojects-the latter being a feature not included in the Yucca Mountain study-to specifY bounding values on the ground motions.

Due to the sensitivity of the computed hazard to the truncation level, the critical point in the definition of this level is that it truly represents the boundary of physically acceptable ground motions and excludes only unphysical val ues. [START_REF] Romeo | Probabilistic versus deterministic hazard analysis: An integrated approach for siting problems[END_REF] proposed an upper bound at two standard deviations above the logarithmic mean of the prediction equation for PGA on the basis that "stronger motions are considered to be unlikely." Since two standard deviations above the mean corresponds to the 97.7th percen tile, it is not disputed that higher levels are indeed unlikely, but this does not mean that they are impossible. As Agathon pointed out, "It is a part of probability that many improbable things will happen"; the obvious corollary to this statement, in the context of upper bounds, is that impossible things will not happen. For an equation with homoscedastic scatter (i.e., constant sigma for all magnitude and distance combina tions), the upper bound will generally lie at least three stan dard deviations above the mean (Figure 2). To distinguish between improbable and impossible levels of ground motion can be achieved only if the physical processes controlling ground motions are identified and the interactions between these processes assessed.

FACTORS DRIVING AND LIMITING EXTREME GROUND MOTIONS

The maximum ground motions that can be experienced at the ground surface are controlled by three factors: the most intense seismic radiation that can emanate from the source of the earthquake; the interaction of radiation from different parts of the source and from different travel paths; and the limits on the strongest motion that can be transmitted to the surface by shallow geological materials. The maximum amplitude of seismic radiation from the earthquake source, for a given seismic moment, is controlled by the total energy release and the rate of energy release, which are dependent on factors describing the mechanics of the rupture process, such as the magnitude of the slip, its velocity (often approximated as a function of the source rise time and the final slip), and the velocity of rupture propaga tion. Although the average values of these quantities generally provide a good first-order approximation, it should be kept in mind that they possess a high degree of spatial variability, resulting from heterogeneities in material properties and stress conditions across the fault plane. This needs to be taken into account since in some instances it is the rate of change of these quantities rather than their absolute values that will influence the level of ground motion. In particular, high-fre quency ground motion results from abrupt changes in rup ture velocity (Madariaga, 1977), whereas lower-frequency motions are influenced more by the actual value of the rup ture velocity. While it is fully acknowledged that slip, rise time, rupture velocity, and their respective rates of change are highly interdependent variables, they are presented separately in the following discussion for the sake of clarity.

For larger events, the energy release is often concentrated in small zones of the fault plane, called asperities (Aki, 1984). Asperities are characterized by having much larger slip than the average over the entire rupture plane; the asperity slip contrast, defined as the ratio of average slip on the asperity to the average overall slip (Somerville et al., 1999), provides a first-order estimate of the relative strength of the asperity. The amplitude of the ground motions generated by an asperity can be expected to increase with asperity size (relative to the rupture area) and asperity slip contrast. The constraints of geometry and energy conservation, however, imply that both these quantities are bounded and moreover that their maxi mum possible values are inversely correlated. This inverse correlation can be observed in practice: In the database of Somerville et al. (1999), the largest asperity slip contrast (3.42) corresponds to the 1984 Morgan Hill earthquake, which has one of the lowest ratios of asperity area to total rup ture (0.14 compared to the average of 0.22). Conversely, the largest relative size of asperities (0.4) is identified for the 1983 Borah Peak earthquake, which has an asperity slip contrast of 1.62, well below the sample average of 2.01.

As discussed previously, details of the slip distribution (e.g., the slip contrast between the edge of the asperity and the surrounding region) might also be important, in particular for the generation of high-frequency motions. The impor tance of these details decreases with increasing distance from the source, however. For sites located far enough away from the source, the effect of the source on ground motions can be satisfactorily estimated using the average value of slip velocity. While the slip distribution reflects spatial variations in the density of energy release, the distribution of the slip velocity gives an indication of the rate at which this energy is released.

The velocity of fault rupture will also play a role in con trolling the most extreme ground motions that can be gener-ated, since rupture velocity affects the corner frequency of the radiated body waves. In a study of the 1906 San Francisco earthquake, Boore (1977) found that a change of rupture velocity from 2 to 3 km/s led to a fourfold increase in com puted ground motions at periods of about 5 seconds. Boore and Joyner ( 1978) investigated how the introduction of inco herence into the smooth propagation model of Boore (1977) affected the outcome and conclude that the sensitivity of near-field ground motions to rupture velocity and azimuth is preserved as long as the mean rupture velocity is the same. It should however be kept in mind that the figure quoted above is the result of a modeling exercise, rather than an observed quantity.

If the spatial variability of rupture velocity across the fault plane is considered, the definition of the maximum per missible values of this quantity becomes even more of an issue. [START_REF] Das | Spontaneous complex eanhquake rupture propagation[END_REF] discussed the development of proposals regarding maximum permissible rupture speeds, addressing in particular the issue of whether rupture velocity can exceed shear-wave velocity. The answer to this problem is to a large extent dependent on the assumptions of the model that is used. Day (1982) used a finite-difference method to study crack propagation in a 3D continuous medium. He found that in cases where supershear velocity is predicted by a uni form prestress model, the introduction of stress heterogene ities is sufficient to reduce the average rupture velocity to less than the shear-wave velocity, while local supershear rupture velocities can still occur in regions of high prestress. Day (1982) concluded that rupture models including extensive segments of supershear propagation should not be considered unphysical, even when the average rupture velocity can be reliably determined to be subshear.

Although the issue of supershear rupture velocity has been considered for some time (e.g., Murray, 1973;Archu leta, 1984), it is still a matter of controversy. One reason for this is that the complexity of the inversion problem makes it difficult to know to what extent the computed variables trade off with each other. For example, in the kinematic model of the 1979 Imperial Valley event presented by Archuleta (1984), the zones corresponding to supershear rupture veloc ities at the bottom of the fault plane are associated with high slip and short rise time. Another consequence of the com plexity of the inversion problem is that it is highly sensitive to the quality of the data used to calibrate it, as illustrated by the recent Kocaeli event. Bouchon et al. (2002) found that the central segment of the fault broke at the supershear speed of about 4.8 km/s. In the same issue, [START_REF] Sekiguchi | Rupture process of the 1 999 Kocaeli, Turkey, earthquake estimated from strong-motion wave forms[END_REF] examined both the hypothesis of supershear rupture and the alternative explanation of a P-wave-triggered asperity to account for the anomalously short S-P time at the Sakarya (SKR) station and concluded that the latter is the preferable interpretation. Therefore, although new evidence in favor of extended supershear rupture (e.g., [START_REF] Bouchon | Observations of long supershear rupture during the magnitude 8 . 1 Kunlunshan earthquake[END_REF] should not be ignored, the issue of the maximum proportion of the fault plane that can undergo supershear rupture still requires further investigation.

In the discussion above, the focus has been on the mechanical aspects of seismic rupture. As a consequence, the spatiotemporal characteristics of energy release are considered from the point of view of an observer located on the source. The advantage of this viewpoint is that it gives a detailed insight into the source processes that will control the ground motion leaving the source. Little information is gained, how ever, on how these source effects interact with propagation effects to ultimately produce extreme ground motions.

An alternative approach is to look at the problem from the receiver point of view, i.e., examine the temporal distribu tion of seismic radiation arrivals at a given geographical loca tion. From this perspective, the amplitude of ground motion at the location of interest at a given moment in time is simply the algebraic sum of the amplitudes of all the waves arriving at this location at that instant. In most cases, the maximal amplitudes of different wave trains will arrive distributed in time; it can also happen, however, that two or more wave trains reach their peak value simultaneously at a given loca tion. This phenomenon is known as constructive interfer ence. Constructive interference can be an isolated occurrence in time, which will lead to a single spike in the amplitude of the motion, without changing significantly the amount of energy that arrives at the location of interest. In the particular case of the arrival of two or more coherent waves, i.e., waves that travel with a constant phase difference, all the peaks and troughs will arrive simultaneously, leading to a significant increase in energy. This phenomenon affects only locations satisfying a number of geometric constraints defined by the location of the elementary sources emitting the waves.

Elementary sources can be of two types: separate areas of the rupture plane (possibly triggered at different moments in time with different rupture velocities) and secondary sources such as reflection or refraction interfaces. As a result, numer ous scenarios can lead to increased ground motions due to constructive interference along the propagation path. Close to the source, constructive interference is mainly the result of simultaneous arrival of contributions from different parts of the rupture plane, as a consequence of the spatial variability of rupture velocity across the fault plane discussed previously.

Examples include forward directivity (e.g., Somerville et al., 1997) and focusing of energy toward a station located along the axis of symmetry of the rupture (e.g., Oglesby and Archu leta, 1997). Also, for an extended source, late arrivals from the part of the fault where rupture initiates can reach the site at the same time as early arrivals from a part of the fault that ruptures later (e.g., [START_REF] Anderson | Expected shape of regressions for ground motion parameters on rock[END_REF]. Farther along the prop agation path, constructive interference is essentially the result of particular geometrical conditions leading to reflection or refraction of the waves in a preferential direction and thus to focusing of energy. Examples include increased amplitudes at the tip of the wedge above the hanging wall (e.g., Oglesby et al., 1998) and topographic effects. The latter can be due to features of the surface topography, such as hills (e.g., Bouchon and Barker, 1996) and canyons (e.g., Boore, 1973). Subsur face topography can also be the cause of constructive interfer-ence; the most prominent example is the interference of surface waves generated at the edge of deep sedimentary basins with direct arrivals, as observed for instance during the 1995 Kobe event. Also, the presence of local heterogeneities in the substratum can result in focusing of energy toward a particular area, as was the case for Santa Monica during the Northridge event (Gao et al., 1996).

As seismic waves propagate to the Earth's surface, other factors act to limit the maximum amplitude of the motion. These factors are associated with the failure of surface materi als, which are usually weaker than the underlying rock, under the loading conditions generated by the passage of seismic radiation. The principle is similar to that of a fuse: Once fail ure is reached at a given depth within the soil profile, the inci dent motion is filtered and no motion larger than the motion reached at the failure stage can be transmitted to the upper strata.

Following a postulate by Schnabel et al. (1972), it is gen erally assumed that the strong part of horiwntal motion in soil deposits is caused mainly by the vertical propagation of SH waves, while the strong part of vertical motion is caused by the vertical propagation of P waves. This simplified repre sentation requires two conditions to be met: First, the soil profile must be almost horizontally layered, a reasonable assumption in view of the mechanics of soil deposition and weathering, and second, the propagation must be vertical. This latter assumption is usually justified by the fact that the stiffness of surface materials decreases toward the surface, causing the wave path to undergo successive refractions and thus become nearly vertical.

Under these assumptions, any soil element within the profile will be submitted in the horizontal direction to cyclic simple shear strains and in the vertical direction to con strained one-dimensional compressional-extensional strains. For these simplified stress paths, constitutive soil modeling (e.g., Prevost, 1978) predicts that a failure condition can be reached only for the simple shear condition. Therefore, it can be anticipated that the horizontal ground motion is limited by the soil strength, whereas the vertical one is fully transmit ted to the ground surface. Such anisotropic behavior has been observed in the field: Aguirre and lrikura (1997) examined accelerograms recorded during the 1995 Hyogo-ken Nanbu earthquake on the vertical array at Port Island, Kobe, where liquefaction occurred. They found that the recorded horiwn tal peak acceleration was only about 25% of the value expected from linear theory, while vertical peak acceleration was close to predictions from linear theory.

As Beresnev et al. (2002) point out, however, the discus sion above holds only if propagation is strictly vertical. This stems from the fact that even a small nonzero incidence angle would significantly increase the contribution from SV waves to vertical ground motion, because they are associated with larger amplitudes than P waves. A spectral study of the verti cal motions recorded during five significant recent events in California yielded the result that for frequencies up to 10 Hz, the vertical ground motion is dominated by SV waves rather than P waves. If this threshold is to be considered robust, the validity of the conclusions reached in the previous paragraph is limited to PGA; for lower frequencies, all components of motion could be expected to be limited by soil strength.

ASSESSMENT OF UPPER BOUNDS

For application in seismic hazard assessment, upper bounds need to be expressed not as single extreme values but rather as functions of the explanatory variables used in ground-motion prediction equations, e.g., magnitude, rupture mechanism, distance, and site condition. A promising approach may be to derive equations to predict maximum ground motions in hard bedrock in terms of magnitude, style of faulting, and distance, and then separately to define the level of this motion that can actually be transmitted through the soil deposits [START_REF] Peeker | An estimate of maximum ground motion[END_REF].

Empirical Data

The first and most obvious tool for exploring upper bounds on ground motions is the ever-increasing databank of strong motion accelerograms. Some proposals have been put for ward to define the upper limits simply in terms of a number of standard deviations above the median values obtained from ground-motion prediction equations, these estimates ranging from 20"to 4u(e.g., Bernreuter etal, 1989;Anderson and Brune, 1999;[START_REF] Romeo | Probabilistic versus deterministic hazard analysis: An integrated approach for siting problems[END_REF]. These deci sions have been based mainly on the representation of the residuals in quantile-quantile plots, which show a deviation from the best-fitting lognormal distribution above a certain level, generally around 20" (Figure 2). The problem with such an approach is that the data at these percentiles are very sparse and the true nature of the distribution of the residuals is dif ficult to determine. For a lognormal distribution of ground motion amplitudes, three logarithmic standard deviations above the logarithmic mean-a level below which few would propose that the physical maximum lies--corresponds to the level of motion that will be expected to occur on average once in every 1,000 occurrences of a given combination of magni tude, distance, and other explanatory variables in the equa tion. In statistical terms, this level of motion corresponds to the 99.87th percentile. To sample the corresponding part of the tail of the distribution reliably, particularly if the model is not homoscedastic and it is assumed that the standard devia tion varies, would require several thousand data points for each combination of magnitude, distance, and site classifica tion.

One possibility is to extrapolate from the distribution of the empirical data, in a similar way to techniques based on extreme value distributions used to infer maximum earth quake magnitudes (e.g., [START_REF] Kijko | Statistical estimation of maximum regional earth quake magnitude[END_REF]. Restrepo-Velez and [START_REF] Bommer | Uncertainty about the uncertainty in seismic haz ard analysis[END_REF] explored an alternative of fitting the data to a truncated lognormal distribution (the Upper Limit Lognor mal, or ULLN, distribution) instead of the traditional unbounded lognormal distribution. This requires the intro duction of an additional parameter, namely the upper bound a. A practical problem is that existing algorithms to fit ULLN (Bezdek and Solomon, 1983) are highly sensitive to the initial choice of a. When the ULLN is fitted to the residuals of a strong-motion data set used for deriving a ground-motion prediction equation, using an approach that determines a as part of the solution, the results generally yielded values corre sponding to about six standard deviations above the mean, which cannot be verified with the existing data. A more fun damental limitation, which is shared by all statistical approaches, is the implicit assumption that the underlying data set is a representative sample of the population of all physically possible ground motions. The limited amount of data constituting the tails of the distribution makes it difficult to examine the validity of this assumption.

The nature of the scatter in ground-motion prediction equations causes other problems in estimating extreme ground motions because it simply lumps together the influ ence of all explanatory variables not explicitly included in the model. For example, some of the largest ground motions in a data set may be due to site-specific conditions such as surface topography or subsurface structure. For a flat site, should the aleatory variability caused by severe topographic conditions be included in the predicted ground motions?

All of this, however, is not to suggest that the empirical data serve no purpose in the search for upper bounds on ground motions. For a start, the maximum recorded motions provide lower bounds on the upper limits; in other words, the upper bounds must be at least as great as the recorded max ima for the same combination of magnitude, style of faulting, distance, and site classification (and possibly other explana tory variables as well). Quite naturally, the sense of what con stitutes a "reasonable" upper bound is conditioned by the difference between the proposed limiting value and the max imum recorded at the time. Figure 3 shows how estimates of limiting values on PGA and PGV have mirrored the growth of the largest recorded values, although a more optimistic interpretation would be that the proposed maximum motions have anticipated values subsequently recorded. Most engineering seismologists would agree that the estimates of PGA obtained in the Yucca Mountain seismic hazard study are too high, being an order of magnitude greater than the largest recorded horizontal accelerations, but no consensus exists on how much they would need to be reduced to be con sidered physically realizable. Should an earthquake produce an accelerogram with a peak of, say, 3 g in the near future, responses to this question would probably be modified.

In passing it is important to note that the maximum ground motions shown in Figure 3 correspond to those recorded on strong-motion accelerographs, which are there fore considered to be in the frequency range of interest to engineering applications. There are examples of much higher ground accelerations recorded from rock bursts in mines. In particular, McGarr et al (1981) report that during a 0. 72 ML tremor, a value of 7.7 g corresponding to a single major S wave arrival was recorded on a transverse component. The corresponding frequency is not given explicitly, but it can be ground conditions, it has been assumed that the estimates refer to rock sites); the maximum PGV estimate of 300 cm/s is taken from Esteva (1 970), who pro posed this value for the near-source saturation of peak ground velocity. In the upper graph for PGA, the stars correspond to proposed limiting values on dif ferent deposits by different studies, as summarized in Table 1 .

inferred from the paper that it was probably of the order of 100 Hz or more.

The real value of the empirical data is in the insight it can provide into both the factors that create extreme levels of motion and those that limit the maximum surface ampli tudes. Careful examination of the largest, and indeed the low est, recorded amplitudes of motion can provide insight into the factors that tend to produce the strongest motions (e.g., Anderson et al., 2002), although it is important always to keep in mind that often a number of factors will be acting in unison to create the higher recorded amplitudes and the temptation must be resisted to assign the cause of the extreme motions to the current fashion, be it site amplification, Moho bounce, rupture directivity, or buried ruptures. Inversion of strong-motion recordings enables detailed pictures of the fault rupture process to be obtained and compared, such as has been done by Somerville et al. (1999) for a number of earthquakes producing significant accelerograms. All of this information is vital input to the only current possibility for identifying upper bounds of rock motions: estimates based on seismological source models.

Source Models

In the absence of large numbers of recordings from very dense accelerograph networks triggered by many earthquakes, the best possibility of constraining the likely maximum levels that ground-motion amplitudes can reach is through the use of models to generate synthetic accelerograms. The usefulness of any particular model in achieving this objective depends on the degree to which it explicitly includes each of the factors that tends to increase or to limit the most severe motions. Even models that do not include all factors may still provide useful insight into the possible ranges of upper bound motions provided their results are carefully interpreted with respect to the effect of any missing factors.

The theoretical modeling of ground motions is based on two fundamental results: Burridge and Knopoff (1964) estab lished the equivalence between a double-couple tensor and the shear-force system acting on a fault plane in the so-called representation theorem, and Aki (1967) found that the far field radiation of seismic sources over a wide range of magni tudes could be approximated using an ol model for the Fou rier displacement spectrum.

The first theoretical models presented for ground motion prediction were simple deterministic models. Haskell (1964) considered the problem of a propagating rupture on a rectangular fault plane, assuming constant rupture velocity along strike and instantaneous propagation along dip. Savage (1966) presented two models to calculate the far-field radia tion of an elliptical fault, rupture being assumed to originate at one focus of the ellipse and then spread radially. In the first model, the slip is considered homogeneous across the fault surface, whereas in the second model it is variable and con strained to zero on the fault edge. Unlike Haskell's model, Savage's models include nucleation and rupture arrest fea tures, leading to the simulated body waves exhibiting two stopping phases. Brune (1970Brune ( , 1971) considered a shear-stress pulse on a circular fault and provided an estimate of the far field spectrum that follows the ol model and depends only on M0, the stress parameter measuring the average strength of the source (within this model), and wave propagation charac teristics.

The deterministic character of the models above implies that they can predict only coherent ground motions. This limits their applicability for high-frequency ground motions, which are strongly incoherent. Hanks and McGuire (1981) showed that high-frequency ground motions could be repre sented as finite-duration band-limited Gaussian white noise. Boore (1983) combined this result with Brune's ol model to obtain a widely used stochastic point-source model. The numeric implementation of this model, SMSIM (e.g., [START_REF] Boore | SMSIM Fortran Programs for Simulating Ground Motions from Earthquakes: Version 2. 0: A Revision ofOFR 96-80-A[END_REF][START_REF] Boore | Prediction of ground motion using the stochastic method[END_REF] also incorporates a theoretical propagation func tion accounting for path and site characteristics via a series of multiplicative filters. Beresnev and Atkinson (1998) extended the method to a finite source by dividing the fault plane into subfaults, each represented as a stochastic point-source, and summing the contributions of the subevents with appropriate time delays.

Zeng et al. (1994) present a composite source model in which a great number of subevents of various sizes are ran domly distributed across the fault plane; the subevent sizes are assumed to follow a fractal distribution (i.e., a power-law distribution with a fractional exponent) and are allowed to overlap. The radiation of each subevent is assumed to follow the shape of the Brune (1970) model. The contributions from the subevents are summed with appropriate time delays. This yields a complex source-time function, which is con volved with a synthetic Green's function obtained using the Zeng and Anderson (1994) algorithm to produce synthetic seismograms.

The models above are all kinematic, i.e., the slip function is imposed. An alternative approach is to impose the stress conditions and derive the slip function from these, as is done in dynamic models. The distinct advantage offered by dynamic modeling is the fact that the limiting strength of the rock, in terms of the maximum sustainable stress, can be specified, thus allowing physical constraint on the radiated motions. In this respect, the upper bounds obtained from the kinematic source models are likely to be overestimates. Dynamic models are based on dislocation theory (£shelby, 1957) and more precisely the propagation of shear cracks. Madariaga {1976) introduced a very simple model based on a circular crack. This simple model was later revised to exclude stress singularities at the crack tip (Madariaga, 1977) and then evolved to a finite-fault model considering asperities (Madariaga, 1979), and more recently to 2D and 3D crack models (Madariaga et al., 1998; Tada and Madariaga, 2001) Aki and coworkers (e.g., Das and Aki, 1977;Aki, 1979;Papageorgiou and Aki, 1983) explored the topic of spatial stress heterogeneity across the rupture surface using barriers and asperities. Although this model yields good results in backward modeling, it requires input parameters that are dif ficult to predict for future events, such as the locations of asperities and barriers. This limitation is related mainly to the entirely deterministic character of the model and somewhat alleviated by the latest generation of dynamic models [START_REF] Guatteri | Strong ground motion prediction from stochastic-dynamic source models[END_REF] that present a stochastic-dynamic model in which the characteristics of slip are allowed to vary across the fault plane, although the resulting model is so far applicable only for low frequencies.

Whichever source model is chosen, the exploration of upper bounds requires, for each magnitude, mechanism, and source-to-site distance combination, a large number of simu lations to be run, varying each of the variables such as the location, size, and strength of asperities, rupture velocity, and nucleation point. This task is onerous in terms of computa tional effort and time; the real challenge is to ensure that all of the combinations of rupture parameters are actually phys ical. There is little doubt that very extreme motions could be computed from a model including a large asperity, occupying about one third of the rupture plane, possessing a slip con trast of 3.0, and exhibiting variations of rupture velocity up to supershear and a short rise time; there is also little doubt that such a situation is unlikely to be realizable. Herein, once again, lies the heart of the problem: the need to identifY source characteristics that are not only unlikely (very small probability of realization) but actually impossible (zero prob ability).

Another reason that the process of searching for extreme values is computationally intensive and very time-consuming is that for a given earthquake magnitude and source-to-site distance, a large number of locations around the fault rupture will need to be considered. With relative ease, for a given velocity model, the least favorable locations could be identi fied, probably controlled by proximity to the major asperity and in the forward directivity position; lateral velocity varia tions would make this task considerably more difficult. How ever, there is also value in exploring the spatial variation of extreme motions around the fault rupture, since this can shed light on how likely the highest ground-motion recordings, invariably from sparse instrument networks, are actually to come close to the upper bounds, even if the rupture parame ters represented a worst-case scenario.

Strength of Surface Materials

Several possibilities exist to account for the influence of the strength of surface materials on the maximum surface ground motion. The first, and most natural, one consists of running site-response analyses with increasing amplitudes of the input bedrock motion (e.g., Mohammadioun and Peeker, 1984). As the input motion increases, the earthquake-induced shear stresses increase within the soil profile, eventually reaching the soil strength at some depth; at that stage, the incoming motion is no longer transmitted to the upper layers and the surface ground motion is bounded to the value reached before saturation of the soil strength. Although attractive, the method suffers some drawbacks: It is time-consuming and requires a genuinely nonlinear algorithm (equivalent linear models are of limited usefulness for such calculations), and the results are very much dependent on the soil constitutive model and the accuracy of the parameters entering these models. On the other hand, using any finite element model, great flexibility is offered to the analyst to model strongly het erogeneous profiles and even two-dimensional geometries. Given the heavy burden put on the analyst for implementa tion of the above approach, several attempts have been made to derive simplified estimates of the maximum surface ground motion. So far, these estimates are limited to the eval uation of the maximum peak ground acceleration. Betbeder-Matibet (1993) has developed an expression for the maximum acceleration at the surface of a uniform soil layer overlying a rock half-space. This solution is based on the relationship between the maximum surface acceleration and the average soil column acceleration; this is obtained by com putation of a fundamental "nonlinear" mode shape of the soil column. Although simple, this solution suffers from some limitations: The soil profile is assumed homogeneous with depth, which is seldom achieved in natural soil deposits, and the maximum ground surface acceleration does not depend on the frequency content of the input base motion. To over come some of these limitations, [START_REF] Peeker | An estimate of maximum ground motion[END_REF] has developed a simple but robust method for estimating the limiting value on PGA by deriving a solution for a heterogeneous soil profile in which the shear modulus varies with depth.

The proposed method is supported by observations in California available at the ROSRINE (Nigbor eta!., 1998) Web site (http:/ /geoinfo.usc.edujrosrine). The method is quite simple, relying mainly on the estimates of two quantities, the depth-dependent shear strength and the yield shear strain of the soil. Of these two quantities, the former can be reliably measured under static conditions using well defined proce dures that may vary according to the soil type, along with the appropriate strength criterion, and variations in the latter do not affect the outcome significantly. This solution involves a number of issues that deserve further investigation, however, such as the margin of error introduced by using static mea surement procedures to estimate dynamic quantities and by assuming a similar behavior at small and large strains.

The possibilities for investigating the limitations on earth quake ground motions imposed by soil strength through labo ratory investigations may be limited. Since the strength and stiffness of soils are pressure-dependent, shake-table tests are limited by the fact that the samples are generally in a 1 g envi ronment and too small to reproduce realistic in situ stresses. Some geotechnical centrifuges, provided that they are able to reproduce realistic seismic ground motions as the input to the soil samples, may be able to provide some useful results.

The existence of a well defined physical limiting mecha nism for ground accelerations on soil is encouraging, but it should be kept in mind that any method trying to capture the resulting limits requires bedrock motion as input, and there fore the problem cannot be entirely decoupled from the esti mation of extreme bedrock motions. The limiting values on other parameters, such as peak ground velocity, which may be controlled by different mechanisms, also need to be explored.

A Framework for Comprehensive Definition of Upper Bounds

As has been mentioned previously, the needs of seismic hazard assessment for critical facilities are not met by simply defining the maximum level of ground motion that could ever be achieved. Upper bounds need to be defined, possibly on a regional basis (and certainly separate definitions will be needed for crustal and for subduction earthquakes) as limiting surfaces defined by at least the following parameters: magni tude, style of faulting, depth of faulting, source-to-site dis tance, and site conditions, which will almost definitely need to consider both the strength of the materials in the upper tens of meters at the site and the structure over a depth of sev eral kilometers. The simple procedure illustrated in Figure 1, of truncating at a specified number of logarithmic standard deviations above the logarithmic mean of the ground-motion prediction equation, is unlikely to be an adequate solution. Among the reasons for this is the fact that in a logic-tree for mulation truncating different attenuation relations in this way will not define a unique surface of bounding values.

Another important point to be borne in mind is that to obtain the strongest possible ground motions at a particular distance from the fault rupture using seismological source models, the high spatial variability of ground motions needs to be considered. Oglesby and Archuleta (1997) argued that the very high amplitudes recorded at the Cape Mendocino station during the 1992 Petrolia earthquake were the result of its location at a specific point with respect to the radiation from the main fault asperity.

In the preceding discussions a clear distinction has been made between the effects that contribute to creating and lim iting extreme ground motions; although the separation is use ful from an analytical point of view, the final model for upper bounds will need to consider the interaction between the two, since the nature of site response is dependent upon the base excitation. For the model of [START_REF] Peeker | An estimate of maximum ground motion[END_REF], this can readily be accomplished because the proposed method prop erly takes into account the frequency content of the base motion through its elastic response spectrum. In addition, simple transfer methods, which are standard practice in struc tural dynamics, can be used to compute not only the maxi mum ground acceleration but the whole surface response spectrum (Igusa and Der Kiureghian, 1985).

A potential pitfall in determining upper bounds is if analyses are performed for individual strong-motion parame ters and then hazard analysis is carried out allowing each parameter to reach its maximum maximus value. The total energy liberated by a fault rupture will always be limited since the seismic energy released per unit of the source volume is approximately constant; extreme ground motions are there fore likely to be created not so much by generation of very high amounts of energy but rather by the focusing of the energy into particular period ranges. Source effects such as rupture directivity will generally affect longer-period motions; site effects will often cause focusing of the energy into a narrow frequency band. Therefore, a vector approach may be needed, following the proposal of Bazzurro and Cor nell (2002), in which PSHA is performed for two ground motion parameters simultaneously, including a model for the covariance of the residuals of each parameter to calculate the joint probability distribution. As has been suggested by [START_REF] Restrepo-Velez | An exploration of the nature of the scatter in ground-motion prediction equations and the implications for seismic hazard assessment[END_REF], this concept could be extended to the definition of upper-bound ground-motions since the latter correspond to the maximum residuals.

Another point that needs to be kept in mind is that although the definition of upper bounds on ground-motion amplitudes will always correspond to the largest absolute val ues, for source effects the most extreme ground motions, in a general sense, will probably correspond to the shortest possi ble duration of shaking for a given magnitude. There may be situations, particularly with regard to liquefaction hazard for example, where long duration of shaking will define the worst-case scenario, but for structural response (except for strongly degrading structures) the worst-case scenario is likely to be related to the amplitude of the ground motion even though it has shorter duration. Short duration of shaking, whereby the major part of the Arias intensity is accumulated in a short interval, can be driven by bilateral rupture and high rupture velocities (Bommer and Martfnez-Pereira, 1999) and by forward directivity effects (Somerville et a!., 1997).

A final point worthy of mention is terminology, espe cially since engineering seismology is a field in which many ill-defined terms-such as attenuation relations and stress drop (Atkinson and Beresnev, 1997)-have become part of common usage. The generic term "upper bounds" used in the title of this paper refers to the final effect of all of the factors that either increase or limit ground motions, whereas it could be helpful to use terminology that distinguishes between these influences. For the maximum amplitudes of ground motion that can be generated in rock due to the limiting val ues of source parameters and wave propagation effects, "bounding ground motion'' has been proposed as suitable, whereas the surface amplitudes imposed by the limiting strength of near-surface strata could be referred to as "upper limits" (C. Stepp, personal communication, 2003).

CONCLUDING REMARKS

This paper has summarized the case for defining upper bounds on earthquake ground motions. A great deal of work remains to be done and this article is not intended to provide definitive answers; rather we have offered some thoughts on ways the problem can be addressed and presented a frame work for how the solutions might be expressed. The number of engineering projects for which extreme motions are likely to be of relevance is limited, primarily by the fact that the importance-in terms of cost or consequences of failure or loss of function-of few projects is sufficient to warrant the expense of designing against the absolute worst-case ground motions. For those projects whose criticality and design life result in design for ground motions with annual frequencies of exceedance smaller than 10-6 , such as the Yucca Mountain waste repository, the issue of upper limits on ground motions requires urgent attention. Although research into this topic may be driven initially by the needs of seismic hazard assess ments for critical facilities, the work can be expected to pro vide very useful insight into issues of more general applicability, defining the limits of the hazard space, identify ing factors that increase or limit hazard levels, and improving, in the process, the capability to model ground motions and site response. E�

ACKNOWLEDGEMENTS

This paper, and ongoing work that it establishes the basis for, has arisen directly from the PEGASOS project in Switzer land, in which most of the authors participated as members of expert panels on ground-motion models (EG2) and site response (EG3). The Technical Facilitator Integrator of both groups was the second author of the paper, assisted by Dr. Patrick Smit for EG2 and by Dr. Martin Koller for EG3, bringing to the project the invaluable experience of the Yucca Mountain study. The authors are very grateful to all the orga nizers of the PEGASOS project who provided such a stimu lating atmosphere to explore these challenging and exciting ideas: Philip Birkhauser, Jim Farrington, Rene Graf, Andreas Holker, Philippe Roth, and Christian Sprecher. We are also grateful to the resource experts, Dr. Arben Pitarka and Dr.

Enrico Priolo, and their colleagues at URS Corporation and OGS Trieste, respectively, for providing numerical simula tions of ground motions that helped guide our first attempts to define upper bounds on ground motions in bedrock.

The authors wish to express gratitude to Luis Fernando Restrepo-Velez for his assistance with Figure 1; to Philippe Roth of Proseis, T. Y. Chang, and Philippe Lussou for input to Figure 2; and to David M. Boore for input to Figure 3. The manuscript has been improved at different stages by con structive criticisms and helpful suggestions from the follow ing individuals, to whom the authors are indebted: Nicholas N. Ambraseys, John G. Anderson, David M. Boore, John Douglas, Thomas C. Hanks, Susan E. Hough, Martin Koller, Art McGarr, Luis Fernando Restrepo-Velez, Carl J. Stepp, and Gabriel R. Toro.

1 .

 1 Civil and Environmental Engineering, Imperial College London 2. Pacific Gas and Electricity Company, San Francisco 3. Geodynamique et Structure, France 4. LCPC, Division MSRGI, and Solid Mechanics Laboratory, Ecole Polytecbnique, France 5. NORSAR, Norway 6.

Figure 1 .

 1 Seismic hazard curves derived using the ground-motion prediction equation of Ambraseys eta/. (1 996) truncated at different levels of scatter. The curves are for a site at 25 km from the boundary of a hypothetical seismic source zone with a maximum magnitude of 7.5 and the bvalue in the recurrence relationship is of 0.7[START_REF] Restrepo-Velez | An exploration of the nature of the scatter in ground-motion prediction equations and the implications for seismic hazard assessment[END_REF]. The upper graph is for an A value of 2, the lower graph for an A value of 3.5.

Figure 2 .

 2 Normal probability plots prepared from the strong-motion data sets used to derive the PGA prediction equations of Berge-Thierry eta!. (2003), Chang eta!. (2001 ), and Lussou eta/. (2001 ). The plots compare the distribution of the residuals with the normal distribution and are gener ally used to check the lognormal assumption. These plots illustrate the fact that it is common to have data points at least three standard deviations above the logarithmic mean.

Figure 3 .

 3 Relations between maximum recorded ground-motion amplitudes in terms of the horizontal components of PGA (uppel) and PGV (lowe!) and the proposed upper limits, as a function of time. In each plot, the solid line shows the maximum recorded motion, regardless of site classification of the accel erograph station. The dashed line shows the largest proposed upper bound for ground-motion amplitudes on rock (where the authors have not specified the

TABLE 1

 1 Proposals for Limiting Values of PGA at Soil Sites

		PGA
	Study	(g) Soil Type
	Ambraseys (1 970)	0.1 5	Very soft marine deposits (PI = 1 0)
		0.30	Inorganic clays of low and medium
			plasticity (PI = 50)
		0.50	Deposits of high plasticity
	Ambraseys (1 974)	0.1 5	Normal ly consolidated clays
		0.35	Highly plastic clays
		0.60	Saturated sandy clays and medium
			dense sands
	Mohammadioun	0.50	Near-source al luvial site
	and Peeker (1 984)
	Dowrick (1 987)	0.36	High p lasticity normal ly consolidated
			clays
		0.61	Medium dense sands and saturated
			sandy clays
		1 .89	Overconsolidated clays