THE SEISMIC ANALYSIS OF SANDY SITES

This paper presents an approach for performing one dimensional effective stress site response analyses for sandy sites, including the evaluation of liquefaction potential. This type of analysis differs from the more common total stress response analyses in that induced pore pressures in saturated sandy soils are accounted for, including the resulting influence on soil properties. This analytical method has been refined to the point where the need for complex and expensive laboratory soil testing is no longer required, a factor which has traditionally held back developments in the effective stress area. The effective stress analysis requires the determination of five soil specific parameters. A trial and error backfitting procedure was developed to successfully determine these parameters from traditional site investigation data rather than detailed laboratory testing. This procedure was investigated using two case studies, the Edgecumbe earthquake of 1987 and the Loma Prieta event of 1989, which both exhibited significant liquefaction damage. The Edgecumbe analysis produced useful results. The predicted ground acceleration required to initiate liquefaction was 2.8 m/s2 (0.29g) which is close to the estimated value of 3 m/s2 (0.31g). This was a good result as a reasonable amount of estimated and correlated data had to be used due to a lack of specific site data. The case study of Treasure Island, in the San Francisco Bay area, also produced encouraging results with both the prediction of liquefaction and surface response spectra in good agreement with recorded data. Both case studies used liquefaction resistance curves determined empirically from SPT blow count data. While this data proved acceptable it was discovered that care must be taken in the use of such overseas derived empirical data, particularly if no corroborating site specific information is available.


INTRODUCTION
The potential for site effects during earthquake events has been recognised for many years.Deep flexible soil deposits tend to amplify source motions and significantly alter the frequency content of the bedrock motions.The extent to which such soft sites exacerbate damage is also dependent on the distance of the source in addition to soil characteristics [22].Numerous cases indicating the influence of site effects have been documented over the years, with widely differing extents of surface damage in areas of close proximity.The earthquakes of Armenia ( 1988), Loma Prieta (1989), Philippines ( 1990) and the extreme case of Mexico City ( 1985) are recent examples of such behaviour.
deep flexible soil sites and bedrock.From this observation, it is clearly important to account for the sub-surface soil conditions when considering the potential earthquake hazard for a site.
The destructive potential of another site effect, liquefaction, has also been recognised for many years, but particularly since 1964 when the Niigata and Alaskan earthquakes graphically illustrated the problem.Since then extensive research has been undertaken to a point at which the fundamental mechanisms involved are now well understood.Liquefaction has been observed to varying degrees in most large earthquake events, with the relatively widespread liquefaction observed during the Loma Prieta earthquake maintaining engineering interest.In the New Zealand context, widespread liquefaction was observed in the lnangahua earthquake (1968), and recently in the Edgecumbe (1987) earthquake and to a lesser extent Westport (1991).There is also significant potential for such behaviour in New Zealand.The most significant and studied areas of this country are Wellington and Lower Hutt.Van Dissen et al. [29] presents a comprehensive microzone study of this region illustrating the influence of soil conditions, and suggests that differences of up to four MM felt-intensity units may be experienced between 1

Civil Engineering Depanment, University of Auckland
This paper presents a computational method for determining the influence of site effects for soil deposits subjected to earthquake loading using a one dimensional effective stress analysis.This models both the propagation of earthquake motions through the soil deposit, and the induced pore pressure response in any saturated sandy soils that may be present.
The pore pressure model in particular requires a number of soil specific parameters to be determined.For most engineering applications, the involved, complex and expensive experimental testing to determine these parameters makes most analyses uneconomic.Generally, only traditional site investigation data is available.Methods of determining the required parameters from such data were investigated and evaluated against experimental and site data from past events.

SITE RESPONSE MODEL
The most common site response analyses used by engineers are one dimensional, as in this case.Non-linear one-dimensional methods have undergone significant development, initially in total stress form (7,10,27], and have subsequently been extended to effective stress solutions [6,13).
The one-dimensional analysis assumes that a soil deposit is subjected to vertically propagating shear waves from the source, and that this propagation through the soil may be described by the one-dimensional wave equation.The solution of this requires knowledge of the constitutive (stress/strain) relationship of the soil.This has been found to be non linear in nature, strain dependent and loading rate independent.This non linear relationship may be modelled using an array of elastoplastic elements which make up a piecewise linear approximation to the soil properties (26,27].Such a system has advantages over other models, such as the hyperbolic function, in that it may be matchf to the actual soil response rather than fitting the soil respon to a predetermined shape.
Effective stress soil properties of saturated sandy soils are accounted for by modelling induced dynamic pore pressures in coupled form with the dynamic analysis.The pore pressure model was originally developed from experimental work by Martin et.al [16], and used numerically in the work of Larkin [11].A coupled form of analysis means as dynamic pore pressures rise, the effective stress of the sand skeleton decreases and the properties of the sand alter as a result.The change in these properties is incorporated directly into the analysis.
The modelling of pore pressures also allows for the prediction of the onset of liquefaction, defined as the point at which induced pore pressure equals the effective confining stress of the sand.This results in a state of zero effective stress and hence zero shear modulus.In reality, dilatancy tends to reduce the pore pressure from this maximum limit, but post liquefaction behaviour is not modelled directly due to the lack of an acceptable numerical model at the present time.For this study a nominal post liquefaction shear modulus was assigned to account for this.
One dimensional models take no account of vertical (end) boundary influences.Two dimensional effects have been investigated by Larkin and Marsh [9] who concluded such influences only become significant for basins with width to depth ratios of less than approximately 6. Thus for wide or relatively shallow basins, the one-dimensional assumption is satisfactory.
As a summary, a schematic representation of the numerical solution procedure for a shear wave propagating through a soil layer is shown in Figure 1.
The solution is performed in the time domain using a numerical finite difference scheme, incorporating pore pressure and constitutive relationship models, by the computer program NESSA [I 3], which was the basis of the work presented below.

EFFECTIVE STRESS PARAMETERS
The excess pore water pressure model is based on the experimental work of Martin et al. [16], who derived a governing expression for the incremental generation of excess pore pressure in sands subjected to cyclic loading.This model requires knowledge of two other quantities, the incremental volumetric strain per load cycle and the tangent rebound modulus.
The volumetric strain expression describes the incremental volumetric strain (compaction) of sand subjected to cyclic loading.The relationship used is one provided by Byrne [ 1], who also provides approximate correlations to SPT blow counts for the two associated experimental constants.
The expression for the tangent rebound modulus describes the elastic rebound of a sand soil skeleton subjected to a reduction in effective stress.This modulus has three experimental constants which are a function of the individual sand grains and hence difficult to correlate to in situ penetration tests.These three constants are undeterminable by "standard" investigation means, and thus require difficult laboratory sample testing.
To determine the rebound modulus constants without the need for sample testing, a trail and error approach was developed involving the use of a numerical simple cyclic shear test simulation [11], which also incorporates the effective stress model.This simulation program is able to produce liquefaction resistance curves for sands at different confining pressures using known effective stress parameters.The unknown rebound modulus parameters may be determined by matching the numerical liquefaction curve to a known site liquefaction curve by altering the three parameters in a trial and error fashion.The two parameters for the volumetric strain component of the pore pressure generation model are assigned from SPT blow count correlations.
The site liquefaction curve may be determined empirically from SPT blow count data using the method of Seed and Idriss [19) or experimentally from the cyclic testing of samples.

CASE STUDIES
Detailed case study investigations using the effective stress method described above are reported by Marks [15], and outlined below.

EDGECUMBEEARTHQUAKE
On March 2 1987, the town of Edgecumbe and the surrounding Rangitaiki Flood plains experienced a magnitude 6.3 (ML) earthquake event.Relatively widespread ground damage and extensive liquefaction resulted.This provided an opportunity to investigate the liquefaction prediction capabilities of this method in a New Zealand situation.As a reasonable amount of site investigation had been perfonned, the available data also allowed investigation of the validity of commonly used CPT/SPT correlation factors, developed primarily from soils in the United States, in a New Zealand setting.
Edgecumbe lies on the Rangitaiki plains, which are situated on the east coast of the North Island in the Bay of Plenty.This area is at the intersection of two fault zones which have fanned the Whakatane Graben in which the plains lie.At the centre of the graben the greywacke bedrock is some 2 kilometres below the surface.The graben is filled with hard compacted alluvial materials, overlain by extensive ignimbrite sheets and recent near surface deposits of volcanically derived alluvium and flood plain material.The area is geologically active, and is thus vulnerable to future significant earthquakes.
A representative soil profile for the Rangitaiki plains to use in the numerical analysis was •detennined from the available data [8,23,28).A single overall soil profile for the region was considered sufficient as investigation of the soil data indicated only limited horizontal soil variations, particularly near Edgecumbe.In this profile, four distinct near surface soil layers were identified.as shown in Figure 2(a).
The data used to derive soil properties included 15 CPT and 4 SPT bores at the Bay Milk Ltd. site near Edgecumbe [28), 0 Loose sandy silt with volcanic/peat/timber deposits 3.5 -Medium to dense fine sands Loose complex stratigraphy of interbedded sands, silts, peats and volcanic deposits Dense to very dense marine sands particle grading curves from the same site, and 4 Seismic CPT logs from the Edgecumbe area [23).The CPT data was used to detennine an average penetration resistance profile since the CPT provides a more detailed account of soil stratigraphy than the SPT.The resulting CPT profile is also shown in Figure 2(b).
No liquefaction data was available for the Edgecumbe area and the empirical approach [ I 9) for detennining liquefaction curves was used.
The bulk of the international data base of liquefaction correlation data derived from past events is from SPT investigations.The CPT profile for Edgecumbe was used to detennine an equivalent SPT profile for use in the empirical generation of liquefaction curves.Since SPT data was available for the same site as the majority of the CPT logs [28), the opportunity was available to investigate of the validity of the CPT/SPT correlation factors.The four SPT and closest four CPT bores were compared, for which the separation distances were less than 10 metres for sites 1-3, and approximately 65 metres for site 4. The profiles using "best fit" correlation factors are shown in Figure 3, with the correlation factors detennined from this procedure given in Table 1.Fi,gure 3(1-4) CPT and correlated SPT profiles for Edgecumbe region particle grading curves for the Edgecumbe site [28], the empirical correlation factors from Robertson and Campanella for each soil layer were determined, and are compared to the average Edgecumbe derived values in Table 2.It can be seen from Table 2 that for the upper two surface layers the empirical and actual correlations are in excellent agreement.
For the lower two layers however the agreement is not so good.The logs of these layers indicate a rather complex stratigraphy of various soil deposits of generally volcanic origin, which may account for some of the disagreement.
The empirical correlation factors are determined from quartz-based sands, and may not be as relevant for those of pumiceous origin.Another reason may be the influence of the soft interbedded layers, creating a smearing effect and influencing the results outside their immediate area.Determination of the exact reasons why the correlations factors were not applicable to the lower layers of the Edgecumbe were not pursued in this study, but is an area well suited to additional research in the New Zealand context.
The main conclusion that may be drawn is the importance of confrrming empirical correlation factors with measured field results where at all possible.From the above results it is clear that sometimes the correlations are useful and sometimes this is not the case.As most empirical liquefaction potential prediction methods are related to standard penetration resistance values, it is very important to be aware of the possibility for such variations, particularly when using CPT data to predict liquefaction resistance.Empirical CPT liquefaction prediction relationships are also based to some extent on correlations with the SPT [17].
Using the field derived correlation factors, an average SPT profile was determined from the averaged CPT profile of Figure 2(b).This was then normalized to an overburden pressure of 100 kPa, and is shown in Figure 4.The third layer was split into two sub layers (3A and 3B) since, although similar in composition, the penetration resistance profile of the lower half of this layer was found to be higher on average than the top half.
The other soil parameters required in the analysis were the densities, coefficients of lateral earth pressure and permeabilities, all of which were estimated.The initial low strain shear modulus was determined from measured shear wave velocity data [23].
Due to the late summer timing of the earthquake, the water table was near the bottom of layer one, and hence although this layer exhibited a low blow count, it did not experience liquefaction.Due to the high blow count of layer two, the probability of this layer experiencing any significant liquefaction effects was very small.The fourth layer showed very high blow counts, and was excluded from the analysis as it had little influence on the liquefaction prediction.
The effective stress properties of layers 3A and 3B were determined by backfitting the numerically-generated liquefaction curves of cyclic shear simulation program to those determined empirically for each sand layer.The numerical and empirical liquefaction curves for layers 3A and 3B are shown in Figures 5 and 6 respectively, and show an acceptably close fit.The Edgecumbe site was then analyzed using the effective stress analysis program NESSA with the soil properties described above.The near-surface soil column analyzed consisted of layers 1 to 3B, which extended to a depth of 14 metres.The input motion at the base of this 14 metre profile was that recorded at the base of the Matahina dam, some 21 km south of the epicentral area.Although relatively distant from the area of main interest, the Matahina dam is underlain by alluvium and tertiary sediments, similar to those underlying most of the Rangitaiki plains.Thus it was considered that this was a reasonable estimate of the general form of the earthquake motion under the plains, although lower in magnitude than the near source area around Edgecumbe.
For the liquefaction analysis of the near source area, the Matahina acceleration record was scaled in accordance with the attenuation relationship derived from this event [4], which takes the form: where a= is the peak acceleration in cm/s 2 I~\~ --------~-   The near source area was assigned a MM IX intensity [14], which corresponds to a maximum surface acceleration in the range of 0.5g to 0.6g, a range also supported by back analysis of slope failures [ 5] .
Initially the input motion was scaled in magnitude until the analysis predicted the onset of liquefaction in the soil column.This onset occurred at a depth of 9 metres, within layer 3A, with a base acceleration of 2.8 m/s 2 (0.29g).This value of base acceleration may be termed the "threshold" acceleration.An acceleration of this magnitude was recorded adjacent to the Matahina dam, where there was no reported evidence of liquefaction.The closest significant liquefaction was reported 3.5 km downstream from the dam in the direction of the epicentral area.Using Dowrick's [4] attenuation curves, a base acceleration of 3.0 m/s 2 (0.31g) would be expected to have generated this liquefaction.Thus the results of the analysis were in good agreement with the observed field evidence.
Considering the conservative nature of the empirical relationships used to determined the liquefaction curves, and hence effective stress parameters, it was expected that the prediction would be more conservative and predict liquefaction at a lower excitation magnitude.Effective stress parameters derived from laboratory liquefaction test data may have provided even closer agreement.On balance the results are in close agreement with the field evidence, and give confidence in the analysis method.
Liquefaction, although generally increasing surface deformations, tends to dampen surface accelerations.Therefore those areas that did not exhibit liquefaction within the near source area would be expected to have experienced the high surface accelerations previously noted.Using a total stress analysis [12] and slightly plastic soil properties to account for these non liquefied near source areas, a base acceleration of around 4.4 rn/s 2 (0.45g) was found to produce surface accelerations in the 0.5-0.6grange.Using the effective stress analysis with this base acceleration to simulate the liquefied soils in the same near source area, the analysis predicted widespread liquefaction occurring in layer 3A, which agrees with the observed events.Figure 7 shows the liquefaction response within the soil column for the above two cases, and for comparison a lower base acceleration of 2.2 rn/s 2 (0.22g) is also included.
Liquefaction is observed where the excess pore pressure generated within the soil column intersects the initial effective overburden pressure line.

TREASURE ISLAND AND THE LOMA PRIETA EARTHQUAKE
The second case studied was Treasure Island, a small artificially-constructed island in the San Francisco bay that experienced widespread liquefaction as a result of the 1989 Loma Prieta earthquake.Due to the high seismic risk of the bay area, many strong-motion recorders were in place at the time of the earthquake, providing the largest ever recorded data set for a large magnitude earthquake.Included in this set was a surface acceleration time history on Treasure Island and an acceleration history on a rocky outcrop to the immediate south of the island.This rocky outcrop is known as Yerba Buena Island.These two data sets provided both a "hard rock" and a "soft soil" record, allowing the opportunity for a detailed investigation of site effects.This was therefore an ideal case for Pore Pressure (kPa) which to evaluate the ability of the analysis method to model both the liquefaction process and the surface seismic response.
Treasure Island was constructed in the mid 1930s.The island has a rock fill dyke enclosure and the surface consists of dredged hydraulic sand and silt.Beneath this surface layer is a natural sand bar, a relatively deep layer of soft "Bay Mud" and various other firm compacted sands and cJays.The soil profile shown in Figure 8 was that used in the analysis, and is based on site investigation data [3].
As in the other case study, many of the soil properties had to be estimated due to insufficient available data.These included: densities, permeabilities and clay compressibilities.Except for the bay mud layer [20], the dynamic constitutive pro.perties of all other soil layers were derived from empirical relationships [21,25).The 14 metre deep surface silty sand layer was the only one susceptible to liquefaction, and the effective stress parameters of this layer were again determined with the use of an empirically generated liquefaction curve.The matched liquefaction curves are shown in Figure 9.
The E-W component of the strong motion recording of Yerba Buena Island was used as the input motion for the analysis of the soil column shown in Figure 8. Figure 10 shows both the recorded and computed surface response spectra for the E-W component, along with the source motion (YBI).
The agreement of the spectra is very encouraging.The computed maximum ground acceleration of 0.185g is close to the recorded value of 0.16g.The most significant difference is the prediction of the lower period components at T ""0.33 s (3.1 Hz) and T=:o0.6 s (1.7 Hz), with underpredictions of 30% and 15 % respectively.Agreement at the mid to long period ranges is very good, with only a slight over prediction of the very long period ranges is evident.
The predominant period of the Treasure Island site is estimated to be 1.3 to 1.5 seconds, which manifested itself as a four fold spectral amplification of the source motion in the surface response spectrum at this period.
The high spectral amplifications evident in the lower period range (T ""0.3 to 0.6s) are in excess of three times the source motion and are evidence of a second predominating response period of the site.These lower period amplifications appear due to the concentration of energy of the source motion in this period region.
To achieve the closeness of fit of the predicted response spectra, the dynamic properties of the surface sand layer were modelled as a sand with a very limited amount of plasticity.Modelling the sand without this feature produced an underprediction of the surface motion.Higher plasticity soils tend to act more linearly, and hence exhibit less hysteretic damping and a lower reduction in shear modulus with increasing strain amplitude.This experience coincides with that of Dickenson et al. [3) who performed a site analysis of Treasure Island using the equivalent linear program SHAKE [ 18).They found that to achieve good agreement with the recorded spectra the sand also had to be modelled as exhibiting slight plasticity.
As no data was available on the properties of the top silty sand layer it is difficult to assess if the treasure island silty sand does exhibit some limited plasticity or if the effective stress analysis underpredicted the surface motion of the liquefied soils at this site.Because of the similar experience of Dickenson et al. [3] and the indication of some silty content (fines) in the sand, the first explanation seems the more likely.

Liquefaction Prediction
As indicated previously, the surface silty sand layer experienced significant liquefaction during this earthquake.The effective stress analysis produced a pore pressure ratio •of 100% occurring at 12.5 seconds into the motion.The recorded surface motion showed a step in the trace at around 15 seconds, widely interpreted as the onset of liquefaction.Thus the prediction is only 2.5 seconds from the observed case. Figure 11 shows the progressive generation of pore pressure over the depth of the surface sand layer, leading to widespread liquefaction at 12.5 seconds.There also appears to be a difference in triggering time of the Yerba Buena and Treasure Island motion recorders of around 1.5 seconds (allowing for wave propagation lime lo the surface and the distance between the recorders).Taking this into account.the difference in the observed and predicted times for liquefaction reduces to only one second, which is a satisfactory result.
The initial motions prior to the arrival of the first high cnngy S wave are a little over-predicted in the analysis due to the small strains induced within the model.With strains of such magnitude.the analysis remains essentially elastic and hence there is little soil damping.
The computed motion after the time of liquefaction is the result of assigning a constant nominal post liquefaction shear modulus (as stated previously) rather than a post liquefaction model computation.The value selected appears to agree well with the recorded data.

COMPARISON OF EFFECTIVE AND TOTAL STRESS ANALYSES
The majority of site response analyses performed, even on sandy sites, are total stress with no account taken of the influence of induced pore pressures.To investigate the influence of excess pore pressure on response analyses.two 25 metre deep sand deposits were subjected to the same 1952 Taft earthquake motion scaled to a peak acceleration of 1 m/s 2 (0. lg).The sand deposits were similar in nature with the exception of their pore pressure generation properties, which were varied for each deposit to give two cases.In the first case, the pore pressure parameters were set such that the sand experienced the onset of liquefaction when subjected to the source motion.For the second case, the pore pressure parameters were altered to allow significant liquefaction over a large zone in the deposit to occur.thus indicating a sand with a higher liquefaction potential.By altering the pore pressure generation properties of the sand rather than the magnitude of excitation to achieve different extents of liquefaction, the influence of excess pore pressure on the response of a deposit could be isolated.
The response of the above two cases was then compared to a non-linear total stress analysis 1121 of a similar deposit.A total stress analysis takes no account of pore pressure generation or the resulting influence on the constitutive relationship.The three response cases studied were therefore: case 1 : onset of liquefaction occurred case 2: significant liquefaction occurred case 3: total stress with no excess pore pressures Figure 15 compares the surface response spectra of the two cases of liquefaction and the total stress case.
From this Figure it is clear that the main influence of induced dynamic pore pressures on the response of a soil deposit is the general lengthening of the periods of the response peaks and the increased damping of the surface acceleration response.An explanation of this may found when comparing the induced strains of the tntal stress and onset of liquefaction cases ( case I and 3), as shown in Figure 16.Max.Shear Strain(%)

Figure 16 Induced shear strain magnitudes
The induced strains are larger in the zone of high excess pore pressures as one would expect, and peak at a depth of 7 metres where the onset of liquefaction occurs.Large strains lead to a higher soil damping, thus retarding the surface acceleration response.A reduced effective stress and higher strains also reduce the average secant shear modulus of the soil deposit, leading to a reduced soil stiffness and an increase in the predominant response periods of the site.The magnitude of the predicted strains at and after the onset of liquefaction are approximate as the model is not formulated to predict behaviour in this region.Further research into post liquefaction sand behaviour and a prediction model for these strains would be very advantageous as a large percentage of liquefaction damage results from significant soil deformations.

CONCLUSION
The aim of this work was to develop a practical method for the use of an effective stress site response analysis which did not require complex sample testing.By backfitting a numerical cyclic shear test simulation to a known liquefaction curve for any sandy deposits present, the unknown parameters may be acceptably determined.
Financial considerations may dictate that the empirical method of Seed and Idriss of generating liquefaction curves be used.However dynamic triaxial testing of samples is the preferred option of determining such curves when quality samples are available.There appears to be a trend overseas favouring increased cyclic sample testing and accepting the costs as a necessary and legitimate expense [2].
The effective stress analysis may be used to predict the potential liquefaction hazard of a site.For example, the threshold earthquake magnitude likely to lead to the onset of liquefaction or the extent of liquefaction for a given design earthquake may be determined.Using such results, the need and extent of ground improvement can be estimated.The improved site may subsequently be re-analyzed to verify it is of sufficient standard.
Both the Edgecumbe and Treasure Island case studies illustrated the ability of the analysis to predict, with reasonable accuracy, the liquefaction and site response of sandy soil sites subjected to earthquake motions.Due to the inherent variability and complexity of in situ soil properties, any theoretical analysis requires a number of assumptions and approximations to be made.This means an "exact" answer will not be obtained, but as the case studies indicate reasonable predictions, ultimately design decisions can be made using this technique.

Figure 1
Figure 1 Solution procedure for a shear wave propagating through a soil layer •

Figure 2 ((
Figure 2(a,b) Soil profile and CPT penetration profile for Edgecumbe region

Figure 4
Figure 4 Nonnalized SPT soil profile for Edgecumbe

Figure 6
Figure 6 Matched liquefaction curves for layer 3B

Figure 7 Figure 8
Figure 7 Edgecumbe liquefaction response of varying base accelerations Figure 9

2 Period (s) 3 Figure
Figure JO Recorded and computed response spectra for 90°component.

)Figure 11
Figure 11 Time history of liquefaction response for surface silty sand layer

Figure 14
Figure 13 Recorded surface acceleration time history

2. 5 •Figure 15
Figure 15 Comparison of response spectra for various levels of liquefaction

Table 1
Calculated CPT/SPT correlation factors