TWO DIMENSIONAL NONLINEAR SITE RESPONSE ANALYSES

This paper presents the results of computer studies of the seismic site response of two dimensional alluvial valleys with a variety of geometries and material properties. The alluvial material is modelled as a nonlinear hysteretic solid and results are presented to illustrate the effect of material nonlinearity on surface ground response. Comparative studies with one dimensional analyses are presented and conclusions drawn as to ground conditions that are appropriate to one dimensional site analyses.


INTRODUCTION
Site specific ground response studies form an important part of the design process for important structures.The last twenty years has seen a large growth in knowledge and applications in the field of seismic ground response analyses.Most studies performed are one dimensional equivalent linear analyses in the frequency domain [9].While analyses are thought to follow the general trend of recorded surface accelerations on a variety of ground conditions, an improvement may be made by using nonlinear time domain solutions.
Thus following the introduction of these equivalent linear solutions, efforts were made to develop truly nonlinear one dimensional methods of site response analysis [ 4,10].These methods model more closely the stress-strain relationship of insitu soils and hence offer an improved method of calculating ground surface response spectra especially for high magnitude events and near source sites.
Two dimensional equivalent linear site response analyses in the frequency domain have been developed [7] and utilised mainly in the nuclear industry in the United States.Nonlinear two dimensional methods were developed by Joyner [3].These methods and modified computer programs form the basis of the work performed for this study.Other multidimensional nonlinear methods have been developed [8] and have been applied mainly in the area of the seismic response of earth darns.
This study was carried out to investigate the nature of two dimensional nonlinear ground response of alluvial basins using a variety of aspect ratios (length to height of valley, L/H).Important conclusions may be drawn from the results about the magnitude and frequency content of the surface waveforms.Comparative one dimensional studies are also presented and conclusions drawn about the appropriate ground conditions for acceptable one dimensional modelling.'Civil Engineeiing Department, University of Auckland.
The growth in this area of research has followed the realisation by earthquake engineers that the subsurface ground conditions play a very important role in determining the characteristics of surface response and hence the seismic loads on structural elements.Classic examples of this effect are the 1985 Mexico and 1989 Loma Prieta earthquakes and the associated local variation in recorded peak ground surface acceleration.

METHOD OF ANALYSIS
The method of analysis and computer programs were developed by Joyner and are described in detail by Joyner [3] and will be presented here in outline.Modifications to the theory were made by incorporating a modified nonlinear soil model, Larkin [6], which is briefly described in this work.
The computational method uses an explicit finite difference scheme (program TENSA) applied over a grid of discrete points chosen to model the valley, as shown in Figure 1.The soil mass is divided into a number of elements, either quadrilateral or triangular, and is underlain by a semi infinite elastic medium of known (large) compression and shear wave velocities.This medium is the source of the seismic energy.The computational method utilises a transmitting boundary to allow wave energy to leave the solution space.The computations are carried out in terms of particle velocity and two independent uncoupled solutions are achieved.

(i)
The in-plane solution (ii) This analysis solves for particle velocities V, and V 3 , as shown in Figure 1, and is referred to as the PSV solution since it utilises the propagation of compression waves and in-plane shear waves.

The out of plane solution
The dependent variable in this analysis is V2, as shown in Figure 1 and is referred to as the SH solution since it utilises the propagation of shear waves causing out of plane distortions.The analyses are carried out in terms of total stress and utilise the mean stress, (Im, and the deviatoric stress (Y,j, as defined by where the repeated index denotes summation and &;; 1s the Kronecker delta.
Mean strain and deviatoric strain are similarly defined.The computations pass over the entire grid at each time step and velocity gradients are used to calculate the new values of strain.
The values of stress are computed from the rheological model, described in the section that follows, from which forces and new velocities are obtained.

Constitutive Relationship
The stress-strain model used is based on classical incremental plasticity theory as described by I wan [2].The one dimensional version of this model has been used by Joyner [4] and Taylor and Larkin (10].Joyner [3] utilised the two dimensional version of the Iwan model, and the framework of this work is used in this study with a modification which is described here.
The model uses a family of yield surfaces in stress space each surface following the yield criteria of von Mises, where k., is a characteristic of the n th yield surface and 0: 0 ,., is the origin of the surface in stress space.The total deviatoric strain e,j consists of an elastic strain eE,i plus a plastic strain component eP,j associated with the outer most yielding surface.This association of eP,j with the outer most yielding surface is a modification to the Iwan model [6).The I wan model as used by Joyner [3] has the concept of a plastic strain increment for each yield surface and sums the results for each surface that is associated with yielding.The modification results in savings of computing time of about 10-15 % .
Kinematic hardening of the Prager type is used such that each yield surface translates in space and remains attached to the stress point on reaching yield until unloading occurs. ( where c., is a constant associated with the ntl' yield surface.The normal flow rule is invoked such that the plastic strain increments are normal to the corresponding yield surfaces.The relationship between mean stress and mean strain is assumed elastic de'" where K is the bulk modulus. The model may be based on laboratory data, such as dynamic torsion tests or cyclic simple shear tests.The values of k., and c0 are then chosen to fit the laboratory data.Alternatively a hyperbolic stress-strain model (or any other suitable analytical or empirical relationship) may be used.The hyperbolic model has been widely used in soil dynamics and is described by Hardin and Drnevich [l].The work in this study was performed using a hyperbolic initial loading curve to represent the dynamic properties of the alluvial valley.The parameters required to specify the material properties are the low strain shear and compression wave velocities Vs and Yr, the undrained shear strength Su, and the mass density p.These values are representative of values commonly measured in alluvial soils.

One Dimensional Analyses
The results of two dimensional studies (TENSA) are compared with a one dimensional total stress nonlinear solution using a computer program DENSOR [5].This technique uses a one dimensional version of the Iwan stress-strain model [10] and an explicit finite difference scheme.This numerical scheme was found to be the most computationally efficient method of solution when compared with implicit and other numerical integration schemes.A transmitting boundary is incorporated [4] and thus the ID and 2D solutions use very similar but independently developed computer codes.This JD method forms a useful means of evaluating the degree of fit between one and two dimensional solutions, where differences in solution technique are minimised.

RESULTS OF ANALYSES
Meshes of varying length to height ratios were constructed with central depths of 100 m and 500 m to model shallow and deep alluvial basins respectively.Material properties were chosen to give a relatively large impedance contrast between the alluvial basin and the bedrock material.To accurately model the frequencies encountered in an earthquake the maximum element dimension with these material properties is 20 m.A typical mesh of length to height ratio (L/H) 8 and central depth 100 m is shown earlier in Figure 1.Shoulders of constant slope 1: 1 were used to simplify any geometrical considerations.The alluvial valley is modelled by the mesh, and the vertically incident seismic motion is applied to all nodes along the base of the mesh with no phase change.

Nonlinear Effects
To study the effects of nonlinearity a basin with a L/H ratio of 8 was used, with central depths of both 100 m and 500 m.This L/H ratio was chosen as one where two dimensional effects contribute only slightly to surface motions at the centre of the basin, and so the effect of nonlinerarity could be studied with confidence.
Two magnitudes of earthquake excitation were applied to the two depths of basin, one large magnitude and one small, and both SH and PSV solutions were investigated.The earthquake excitation used was that recorded at the Castaic Ole! Ridge Route Station during the San Fernando earthquake of 1971, having both horizontal and vertical components of motion and being recorded on a rock outcrop.These earthquake records were scaled to provide the two magnitudes of excitation, the larger magnitude being a factor of 8 times the smaller one with a peak acceleration of 3.09 mls 2 i.e. 0.315g.These two magnitudes provided large strain and small strain solutions, with strains in the nonlinear and the elastic ranges of the soil model respectively.
Results have been produced for both the SH solution and the PSV solution for three angles of wave front incidence, 0°, 20° and 45° to the vertical, but for reasons of space constraints only results for vertically propagating waves are presented here.
Figure 2 shows results for vertically propagating SH waves for both magnitudes of excitation and both basin depths.The traces plotted are surface velocity time -histories calculated at the centre of the alluvial basin.It is important to recognise the difference in frequency content of the motions for the two strain levels.This shows that the basin is responding to and hence filtering the incident waves in a strongly strain dependent manner.The consequence for structures is the difference in the fundamental frequencies of the surface motions.This is one of the central reasons for developing nonlinear hysteretic methods of site response analysis.The stronger seismic motion produces surface motions of lower frequency and higher amplitude which are likely to influence large flexible structures.The magnitude of the surface motion is not a linear function of the incident bedrock motion.Hysteretic strain dependent damping in the nonlinear model will result in larger energy loss in the high magnitude earthquake with associated relative reduction in the ratio of the surface to bedrock velocity amplitude compared to the small magnitude event.Similar results were found for the PSV solution for both the horizontal and vertical components of in-plane motion.Response spectra were calculated for surface motions at several points across the basin, however only results at the centre of the basin are presented here.Figure 3 shows the results for a 100 m deep basin being subjected to vertically propagating SH waves.Considerable energy is present at longer periods for the higher strain level, and there is not a single scaling relationship between the magnitude of acceleration at a single period for the two strain levels.Figure 4 is the horizontal in-plane response (PSV solution) for a basin of central depth 100 m, and in Figure 5 the vertical in-plane response for a basin of central depth 500 m is shown.These results show similar features to those discussed above for SH waves.In the case of vertical motion more energy is evident in the high frequency area of the spectrum, as is usual!y the case.

Two Dimensional Effects
For this part of the study the more intense earthquake motion was used (0.315 g).The length to height ratio (L/H) of the mesh was varied to investigate the effects of two dimensional topographies.Two depths of valley were modelled as before, for both SH and PSV solutions.Only results for vertically propagating waves are presented in this paper.
Initially for a large value of L/H the results of the 2D analysis were compared with those from the lD analysis [5].This 1D analysis is limited to SH motions and so the comparison is restricted mainly to this solution.Figure 6 shows the surface velocity time histories produced by both methods of analysis for a basin depth of 100 rn and in the 2D case at the centre of a basin with a L/H ratio of 16.The results compare extremely closely.For this alluvial basin, with a large L/H ratio, it would be a very good approximation to use a ID analysis to compute the motions at the centre of the basin.
For large values of L/H the horizontal response of the PSV solution compares well with the SH solution.This illustrates the principle of using plane strain approximations for 3D configurations.Convergence with the lD solution also occurs for the PSV solution for very wide basins, thus a 1D analysis is appropriate.
Two considerations of two dimens1011al topographies were investigated.First, at what UH ratio do the surface motions at the centre of the basin diverge from the ID solution; second, at what L/H ratio do significant variations in surface response occur across the valley.These two considerations lead to a decision of where a lD analysis is an acceptable approximation.Surface velocity time histories, acceleration response spectra.contours of strain and surface velocity fourier spectra plotted at twenty points across the basin were considered.Contours of strain for a 500 m deep basin subjected to vertical SH waves and L/H ratios of 8, 6 and 4 are shown in Figure 7.For a L/H ratio of 8 the strain profile at the centre of the basin and beyond the quarter points compares almost exactly with the lD solution.
However, for L/H of 6 the profile at the quarter points starts to diverge from the lD solution, and for L/H of 4 even the profile at the centre of the basin is different.Similar results are obtained for the PSV solution.
Figure 8 shows the variation in surface velocity fourier magnitudes with frequency and position across the basin for five values of L/H.These results are for an SH wave solution and a 500 m deep basin.A decreasing L/H ratio results in the two dimensional effects of the basin becoming more noticeable, changing the computed predominant frequency.

CONTOURS OF STRAIN
Figures 9 and 10 show acceleration response spectra at the centre of the 100 m deep basin for the SH solution and the horizontal component of the PSV solution respectively.The effect of varying the L/H ratio is marked at all natural periods, with the maximum ground acceleration (at zero period) being larger for the smaller L/H ratio in both cases.
It can be seen therefore that the two dimensional effects of the basin configuration may be an important factor in a site response analysis.It appears that for alluvial basins of L/H ratio 6 or less this is the case, and for such basins a JD analysis may produce inaccurate results.

CONCLUSIONS
The effects of nonlinearity on the seismic response of an alluvial basin are important.For a range of depths the response has been found to be strongly strain dependent.Stronger seismic motions produce larger surface response at lower frequencies.
The nonlinear hysteretic soil model results in greater energy loss for high magnitude seismic motions compared to smaller magnitude events.
Considering surface response near the centre of the valley, studies have shown that one dimensional analyses are divergent from two dimensional results for length to height ratios of approximately 6 or less.When calculating site response nearer the edges of the valley even greater care should be taken when using a one dimensional analysis.
The method presented here for calculating two dimensional site response analyses (TENSA) has been developed to the point of practical application to large projects.The depths of alluvium likely to be encountered will require modest computation times.
Many other applications have been investigated with respect to this solution technique.These include non-vertically incident seismic waves, non-symmetrical alluvial basins, surface wave generation, soil-structure interaction and several case studies.
Fig.I Typical basin mesh detail.
Fig. 4 Fig. 5The vertical response to vertical PSV waves for a 500 m deep basin.

Fig. 10
Fig. 10Horizontal acceleration response spectra for vertical SH waves for the centre points in a 100 m deep basin for two L/H ratios.