Linear and Nonlinear Behaviors of Soft Soil Layers Using Lotung Downhole Array in Taiwan

We analyze the acceleration records of a vertical array to study the linear and nonlinear behavior of the soft soil layers at LSST site, Lotung, Taiwan. This array includes five triaxial accelerometers deployed at depths of 0, 6, 11, 17 and 47 m. During a 6-year operation from 1985 to 1990, 30 earthquakes (4.0 s; ML s;6.5) triggered this array. The maximum PGA value is 258 gals recorded at the surface station. Spectral analyses show that the strong motion (PGA>150 gals) causes the peaks of ratios to shift to lower predominant frequencies. The averaged spectral ratios of 15 well-recorded weak motion data (PGA <100 gals) are selected as a reference; the shift of the maximum predominant frequency can reach 20%. Compared with the weak motions, the strong motions also decrease the amplification factor. The maximum reduction of the amplification can reach 50%. The results of waveform simulation show that the linear model based on the Haskell method can well predict the weak motions at various depths. However, this linear model does not work for the strong motion data. Here a nonlinear numerical scheme, such as DESRA-2, is required and can significantly im­ prove the simulation results although the PGA value at the surface station is still underestimated. Overall, the nonlinear numerical calculation is fea­ sible to predict the strong motions for a horizontal layered structure. (

fication of the seismic waves originates from the strong contrast between the physical proper ties of the rocks and the sediments.To evaluate this amplification, the seismic response of the soil is treated as a linear behavior under low levels of strain.But for larger stress-strain levels, the results of laboratory testing of soil samples showed a nonlinear relation that represents the nonlinear character of the soil response.Some authors have been trying to find the observa tional evidences of the nonlinearity from seismological data and estimate how much it influ ences strong ground motions (Chin and Aki 1991;Beresnev et al. 1995a;Beresnev et al. 1998a, b;Su et al. 1998;Cultrera et al. 1999).In those studies, the nonlinear effect causes the reduc tion of waveform amplitude in the time domain and the shifting of predominant frequencies and peak reduction in the frequency domain.They are due to the nonlinear response of the material that causes change in the elastic properties of the medium dependent on the waveform amplitudes.Aguirre and Irikura (1997) studied the nonlinearity, liquefaction, and velocity variation of soft soil layers in Port Island, Kobe, during the 1995 Hyogo-ken Nanbu earthquake.The S-wave velocity structure before and after the mainshock was found to be different.Specifically, the S-wave velocity of the second layer (5 to 16 m depths) after the mainshock was 20% lower than before.The liquefied state remains at least 3 hr after the mainshock but no more than 24 hr.The strong influence of nonlinearity during the mainshock yielded a signifi cant reduction (25 % ) of the horizontal ground motions at surface.
Several techniques have been used to detect the nonlinear effect.One is spectral ratio evaluation of the observed data between surface and bedrock during strong and weak ground motions (Ordas and Faccioli 1994;Beresnev et al. 1995a;Hartzell 1998;Su et al. 1998).An other option for estimating the soil layer effect is to use the recordings from a vertical array of seismometers (Wen et al. 1994(Wen et al. , 1995;;Aguirre and Irikura 1995;Beresnev et al. 1995b;Satoh et al. 1995;Borja et al. 1999Borja et al. , 2000)).The reduction and I or shift of the peaks during strong motion are indications of the nonlinearity.Another technique used to evaluate nonlinearity is based on the comparison of the observed ground motions during strong motion with those simulated by a linear method.The difference with the observed data can be interpreted as nonlinearity.Two commonly used linear methods are the ID Haskell method and the empiri cal Green's function method (e. g., Aki and Irikura 1991;Aguirre et al. 1994).Numerical approaches to predict the nonlinear response of soil can be classified as either an equivalent secant approach (e.g., the SHAKE program by Schnabel et al. 1972) or a direct nonlinear approach (e.g. , the DESRA-2 program by Lee and Finn 1982, the CHARSOIL program by Streeter et al. 1974 and the SPECTRA program by Borja and Wu 1994).Yu et al. (1993) used a direct nonlinear approach, DESRA-2, to examine the differences between linear and nonlin ear soil response with various levels of base excitations.They showed that in strong excitation the soil nonlinearity causes deamplification and also a shift in peak frequencies to lower val ues for an unsaturated shallow soil deposit of 20 m thickness.Ni et al. (1997) used a modified version of DESRA-2 constitutive model for saturated soil to study the nonlinear seismic re sponse including liquefaction of medium dense soil deposits of various thicknesses.By using the stress-dependent model with impulse base excitation, the nonlinear behavior of various soil deposit shows lower deamplification and higher first-mode (resonant) frequency than that of the stress-independent soil properties model.
In this article, we use the recordings at the Lotung vertical array to study the linear and nonline ar behavior of soft soil layers.First, the spe ctr al ratio is calculated by dividing the spectrum at surface relative to that at the depths of 6, 11, 17 a nd 4 7 m for the strong motion and weak motion events.Based on the shear-wave velocity structure by uphole and crosshole shooting methods, we simulate the 20 May 1986 ear thquake sequence incl uding the foreshock, after shock and mainshock events for the linear and nonlinear models.The nonl ine ar numerical scheme, DESRA-2, is al so used to predict the ground accelerations of the strong motion events.

SITE AND DATA
The dow nhole arr ay (represented by the star symbol in Fig. 1) is in the Lanyang alluvial de posits at the LSST (Lotung large-Scale Seismic Te st) site in the southwest quadrant of the SMART!array, Lotung in northeastern Taiwan.The sediments are mostly composed of interlayered silty sand and silty cl ay bed with gravel (Anderson and Tang 1989).This site is geotechnically classified as "deep cohesionle ss soil site" .A geophysical measurement pro gram was conducted to de termine P-wave and S-wave velocitie s (HCK 1986;Anderson and Tang 1989).The shallow S-wave velocities are estimated by crosshole and uphole shooting methods and tabulated in Table II of We n's paper (Wen 1994) depth did not catch any events, therefore, only the recordings at five depths are discussed here.
Digital data were recorded as 12-bit words at the rate of 200 samples per second.During a 6year operation from 1985 to 1990, 30 earthquakes ( 4.0:::; ML :::; 6.5) triggered this array.Figure 3 shows the epicentral locations of these earthquakes.The relative ground motion parameters, PVA and PHA values at the surface station FAl-5 for these 30 earthquakes are tabulated in Table 1.The maximum PGA (peak ground acceleration) value is 258 gals recorded at the surface station FAI-5 (Event no.4).Here events with peak ground acceleration at the surface larger than 150 gals are attributed to the "strong motion" and less than 100 gals are considered as the "weak motion".The ground acceleration between the two classes is the transition zone where the nonlinear effect begins.
On spectral ratio analysis, a 5 or 10-sec window, which primarily brackets the large am plitude S-wave parts of seismograms, is used with a 5% cosine taper.We choose the 10-sec window for strong motion and the 5-sec window for weak motion.Fourier amplitudes are calculated for all accelerograms used.In order to avoid the pseudo-peak, twenty consecutive smoothing using a 3-point running Hanning average is applied to the raw spectra and then spectral ratios are calculated.We chose 19 events including 4 strong motion events (marked with "S" in Table 1) and 15 weak motion events (marked with "W" in Table 1) based on record quality and completeness.Since all the waveforms used are near-field recordings, the signal noise ratios are all above a factor of 10 that errors in the spectra introduced by the noise may be insignificant.

TECHNIQUES
Site effect can be considered as empirical transfer functions of the surficial layers.Several methods have been used to study the problem.Two techniques commonly used are the stan dard spectral ratio and the numerical simulations.Both are considered here and their results compared in the evaluation of the site effect.The standard spectral ratio, S T , can be calculated by dividing the horizontal Fourier spectrum of the ground motions recorded on an alluvium site, SHS, by that recorded on a nearby rock site or at the downhole , S H8 • s = S HS T S . (1)

HB
The latter station is considered to be a reference station.The spectral ratios are taken at surface relative to the different depths of 6, 11, 1 7 and 4 7 m.Basically, the ratios have some variations due to the difference of the incidence angles and the heterogeneify of the medium.The transfer function could be considered as an average of ratios for several events at each depth.
For the numerical simulations, they are divided into two parts.One is the ID Haskell method (Haskell 1953(Haskell , 1960(Haskell , 1962)), based on the linear model, for the weak motions.Huang (1994) revised the technique to simulate the ground motions at the different depths using that at surface.The best advantage is that the upgoing and downgoing waves need not to be distin guished and separated from the input motion.It can be used to directly simulate the whole waveform at downhole.The other is the DESRA-2 numerical scheme, based on the nonlinear model, for the strong motions.

Haskell Method (SH wave)
For a given laterally homogeneous layered medium, if the thickness and the shear veloc ity of the medium of layer m are dm and /J m , the ground motions at the bottom and the top of layer m can be correlated by a matrix (Haskell 1960 where [A] is a 2 X 2 matrix given by taking the product of all the propagating matrix between the free surface and the n1b interface, i.e.
( 4) On the free surface, the shear stress vanishes and (3) can be further simplified as NOTE : ( 1) EQ type : S and W represent strong and weak motion event.
where A1 1 is the element of matrix [A] at row 1 and column 1.For a given subsurface structure, el ement Au is a very stabl e value even for a zero damping material .Therefore, ( 5) is a stabl e formul a for calculating the downhole displ acement using the surface displacement.

DESRA-2 Numerical Scheme
The program DESRA-2 developed by Lee and Finn (1982), that is a direct nonlinear computational method, used to calculate the nonlinear soil response.In the DESRA-2, the initial loading phase is defined by the hyperbolic str ess-strain relationship in which y is the strain ampl itude, r is the corresponding shear stress , Gmo and rm o are the initial maximum shear modulus and shear stress.Thus as the strain increases, the stress tend ing to restore the medium also increases but is asymptotic to a maximum level.For ver y large value of r mo the stress-strain rel ati onship is approxim atel y linear, while for smaller values of rm0 the stress-strain relationship becomes more severely nonlinear as the strain amplitude y becomes lar ger.In the subsequent unl oading and reloading phases, the stress-strain rel ation ship is given by the Massing stress-strain curve (Masing 1926;Lee and Finn 1982) that is described by ( 7 ) where Gmr and rm1 are the maximum shear modulus and shear stress during the unloading and reloading, which are changing continually due to the development of pore-water pressure and strain hardening dur ing unloading and reloading.In this study, we ignored the effects of pore water pressure and strain hardening so that the maximum shear modulus Gmt and shear stress rm1 dur ing the unloading and reloading ar e equivalent to the initial maximum shear modulus G mo and shear stress r mo.Two parameters, r r and y , , ar e the shear stress and strain corre sponding to the turning point of the hysteresis loop.As implemented in the program, a layered soil model, in which properties depend only on the depth, is excited by the horizontal shear wave propagating vertically.
According to Seed and Idriss (1970), the quantities Gm o and Tmo were given val ues that are representative of a medium dense sand with relative density of about 65%.In units of kN/m2, they can be expressed as and ( 9) where a: is the mean principal effective stress, a: is the vertical effective stress, K0 is the coefficient of lateral stress at rest, and <fl is the effective angle of shearing resistance.Values for <f>' and K0 were selected to be 35° and 0.43, which are typical values for medium dense sand.The mean effective stress can be computed as (10) The parameter, rmo• is related to the strength envelope of the soil and can be easily derived from a Mohr circle failure envelope as given by equation ( 9).

SPECTRAL ANALYSIS
The results of spectral ratio analysis are described first.We calculate spectral ratios at surface relative to the depths of 6, 11, 17 and 47 m. Figure 4 shows the RMS results during the May 20, 1986.earthquakesequence (No. 6, 7 and 8 tabulated in Table 1) .RMS is the root mean square of two horizontal component results.The thin, thick and dot-thin lines individually represent the spectral ratio of the foreshock (No. 6), mainshock (No. 7) and aftershock events (No. 8) of this strong motion sequence.Compared with foreshock (4.5 Hz, 3.0 Hz, 2.2 Hz and 1.3 Hz) and aftershock (4.2 Hz, 2.9 Hz, 2.2 Hz and 1.2 Hz), mainshock event (2.8 Hz, 1.7 Hz, 1.2 Hz and 0.85 Hz) apparently has lower fundamental predominant frequency and smaller amplification factor.Four frequencies in the parentheses represent the fundamental predomi nant frequency at surface relative to the depths of 6, 11, 17 and 47 m for that event.Many researchers attribute the phenomenon to the nonlinear behavior of soil after a large earthquake.It is similar to the theoretical calculation predicted by Yu et al. (1993).Besides, for the after shock event, in 12 minutes after mainshock, the soil property is recovering but its predominant frequency at the first layer is still lower than that at the foreshock event.It probably can be explained that the return time, only 12 minutes, was too short to go back to the initial state of the soil and then the aftershock followed closely.Furthermore, to investigate the changes of site responses between the strong motion and weak motion, we chose 15 weak motion events and the average spectral ratios of them are as the reference.Figure 5 shows the comparison between the spectral ratio of the mainshock (thin line) and the average spectral ratio of the 15 weak motion events (thick line) at surface relative to the depths of 6, 11, 17 and 47 m.The light thin lines represent one positive and negative standard deviation on the average weak-   motion spectral ratio.For the weak motion events, the fundamental frequency at 0/6, 0/1 1, 0/17 and 0 m/47 m individually appears at 4.9 Hz, 3 Hz, 2.2 Hz and 1.2 Hz.When the fre que ncy band is less than lOHz, the strong motion event has lower fundamental fre quency than the weak motion events at the vari ous de pths .The largest shifting between them is about 20%.
The amplific ation factor al so de creases for the strong moti on event and the largest difference betwee n strong and weak events is about 50%.Whe n the frequency range is between 10 Hz and 20 Hz, the ampl ification factor of the strong motion becomes higher than the weak motions.Besides, these two lines intersect at about lOHz.The results are similar to Wen et al. (1994Wen et al. ( , 1995)).6 where the black thick line is also the average weak-motion spectral ratio and the other four lines represent the results of different strong motion events.It is appar ent that all of these strong motion events have lower fundamental frequency than the weak motion events at the various depths.

NUMERICAL SIMULATIONS
In numerical simulation, it can be divided into two parts.One is the lD Haskell method based on the linear model.Another one is for the nonlinear model, that is, the DESRA-2

DHB (RMS)
-AVG.OF WEAK EVENTS -1986,0 1,16 13: 04 -1986,07,30 11 :31 -1986,05,20 05:25 -1986, 11,14  FREQUENCY (Hz) 20 numerical scheme.We will discuss them as the following.The Lotung downhole array is located at the Lanyang alluvial plain.The shallow soil structure around the site can be represented by a horizontal layered structure (Wen and Yeh 1984).In order to study site amplification, we use the Haskell method (Haskell 1953(Haskell , 1960(Haskell , 1962;;Huang and Chiu 1996) to simulate the ground motion of the horizontal and layered structure at the Lotung downhole site.The ground motion at the surface site is chosen as the input to predict those at the depths of 6, 11, 17 and 47m.The great advantage of this approach is that the upgoing (or direct S) wave need not be separated from the downgoing (or reflected)

Linear Model
wave.By means of this technique, the whole waveform corresponding to each layer can be calculated directly.According to Wen (1994), the used velocity structure and density are listed in Table 2.No attenuation term (l/Q , =0) is added to the simulation pro cess.This is mainly because the propagating distance is short (47m) and the bandpass used has removes high fre quency components which would be most affected by Q 5 In this part, the recordings of the 20 May 19S6 earthquake sequences are analyzed.Ac cording to the back-azimuths of the earthquakes, two horizontal components are used to con struct the SH wave.The velocity structure is divided into 11 layers and is considered to be at the half-space under the 47m depth.It is assumed that the SH wave is vertically incident from the half-space.The ground motion at surface site is used to predict those at downholes.Figure 7 shows the synthetic results (thick lines) compared with the observed data (thin lines) for (a) the 29 March 19S6 earthquake and (b) the S April 19S6 earthquake.The synthetic waveforms at the depths of 6m, 1 lm, 17m and 47m match the observed in amplitude and in phase.The simplified model provides a good simulation for the weak ground motion recordings.However, for the aftershock event (No .S), the synthetic waveforms match the observed one neither in amplitude nor in phase (Fig. Sa).This indicates that the shallow soil behavior may be changing.The fundamental resonance frequency of one layer isf=v.J4H, where v s is the S wave velocity and His the layer thickness.According to this formula, we can roughly estimate the predomi nant frequency for the horizontal layered model.Therefore, a modified S-wave velocity struc ture is constructed.That is the S wave velocity at the first layer is adjusted from 120 m/sec to 95 m/sec.Figure Sb shows the synthetic waveforms for the modified velocity structures.Both fit better in amplitude and in phase.But for the mainshock event (N o. 7), the original model (Table 2) does not produce a simulation of the strong motion event that matches the data well (Fig. 9).Especially at the depth of 47m, the synthetic result is bad.However, the result can not be improved by means of adjusting velocity or incident angles.It means that the linear model based on Haskell method does not work well for the strong ground accelerations (PGA > 150 gals).This phenomenon is probably caused by the nonlinear effect of the soft soil after strong  -40 '---'---'---'-----' ---'---'---'-----' ---'----'---'--__J  � -20�� ��'--�--'-� �.L. .-�--'-��-'--�--'-� �-'---'---'   gro und shaking.In order to solve thi s problem, the nonlinear numerical method (DESRA-2) is used next to predict the stro ng ground motions.

Nonlinear Model
To examine the DESRA-2 numerical technique , we did some tests to re produce the re sults of Yu et al. (1993).According to Wen (1994), the velocity structures and some input parameters fo r the nonline ar mode l are listed in Table 3 where Gmax = pv;, rmax = 0.0005Gmax• and b value is relative to 5% damping ratio (Yu et al . 1993;Ni et al. 1997).First, the recording at the 47 m depth is selected to be the input motion to predict the gro und motion at surface and other depths.Figure lOa shows the synthetic results (thick line s) compared with the observed data (thin lines) for the 20 May 1986 mainshock event (Table 1, No. 7).The amplitude at surface is underestimated about a facto r of 0.5, but the y still have good fits in phase.It indicates that some errors appear in the velocity structure or the no nline ar behavior probably happens in the more shallow structure after strong motion.Therefore, we change the recording at the 17 m depth as the input motion.Figure 1 Ob shows the simulate d re sults co mpare d with the observed data.Apparently, there are better agreements in both am plitude and phase .We al so ge t the stress-strain relationship and stress variation wi th time in this case.Figure 11 shows the stress-strain re lationship and stre ss variatio n with ti me at the 10 m depth for the mainshock event.Before 10 seconds, the strain value is small and the stress strain relationship is approximately line ar ("A" part in Fig. 1 la).As the strain becomes higher, (2) 1 Pa"' 1 N/m 2 519 the stress-strain relationship follows the Masing curve ("B1" and "B2" part).The "C" point is where the maximum stress and the maximum strain happen at about 10 seconds.After that, the strain decreases again, the stress-strain relationship also turns to be a linear trend ("D" part).
The level of strain value is between -0.09%-0.06%.According to the above discussions, the DESRA-2 NONLINEAR MODEL  soil has produced evident strain changes during the mainshock event.Therefore, the nonlinear numerical scheme, DESRA-2, is required and can significantly improve the simulation results although the PGA value at the surface station is still underestimated.Overall, it is feasible to predict the strong motions for a horizontal layered structure.

CONCLUSIONS
In order to study the linear and nonlinear behavior of the soft soil layers, the ground accelerations recorded by the DHB vertical array at LSST site in Lotung are analyzed in the paper.There are 30 earthquakes detected by this array during a 6-year operation.The maxi mum PGA value is 258 gals recorded at the surface station.Spectral ratios at surface relative to the depths of 6, 11, 17 and 47 m are calculated.Spectral analyses show that the strong motion (PGA > 150 gals) causes the peaks of ratios to shift to lower predominant frequencies.Many researchers attribute the phenomenon to the nonlinear behavior of soil after a large earthquake.When the averaged spectral ratios of 15 well-recorded weak motion data (PGA < 100 gals) are selected as a reference, the shift of the maximum predominant frequency can reach 20%.In comparison with the weak motions, the strong motions also decrease the ampli fication factor.The maximum reduction of the amplification can reach 50%.It indicated that the strong motion events have induced the nonlinear effect of the soft layer at LSST array.
The shallow structure around the LSST site can be represented by a horizontal layered structure.The SH wave is vertically incident.The results of waveform simulation show that the linear model based on the Haskell method can well predict the weak motions at various depths.However, this linear model does not work well for the strong motion recordings.In order to solve the problem, the nonlinear numerical scheme (DESRA-2) is used and can im prove the simulation results.Basically, the synthetic result is consistent with the observed data in amplitude and phase although the PGA value at surface is still underestimated.Overall, the nonlinear simulation is fe asible to predict the strong motions for a horizontal layered structure.
Fig. 1.A map show s the location of the Lotung downhole array (star symbol) that is inside the LSST (Lotung large Scale Seismic Te st) site in the Lanyang alluvial plain.The site is in the southwest quadrant of the SMARTl ar ray, Lotung, northeastern Taiwan.

Fig
Fig. 2. A schematic cross-section of the site, the deployment of instruments at the depths of 0, 6, 11, 17 and 4 7 m and the velocity profiles of the shear wave.

Fig. 5 .
Fig. 5. Comparisons between spec tral ratio of the 20 May 1986 mainshock (thin line) and the average spectral ratio of the 15 weak-mo tion events (thick line) at surface station relative to the depths of 6, 11, 17 and 47 m.The light thin lines represent one positive and negative standard deviation on the average weak-mo tion spectral ratio.
Yu et al. (1993) also got the similar phenomenon using th _ e theoretical simulation.It indicates that the strong motion event induces the nonlinear behavior of the soft soil at LSST array.Besides, the similar results also appear on the other three strong motion events (No. 4, 12 and 17) shown in Fig.

Fig. 6 .
Fig. 6.Comparisons between four spectral ratios using strong motion events (No. 4, 7, 12 and 17; represented by dif ferent kinds of lines) and the average spectral ratio of the 15 weak-motion events (black thick line) at surface station relative to the depths of6, 11, 17 and47 m.
Fig. 7a.The synthetic ground acceleration (t hick lines), using Haskell method, compare d with the observed data (thi n lines) at the Lotung downhole site in the 29 March 1986 earthquake (No. 5) .The recording at surface is selected as the input motion and the SH wave is vertically incident from the half-space.HASKELL SYNTHETIC WAVE
Fig. Ba.The synthetic ground acceleratio n (thick line s) , using Haskell method with the velocity model in Table2, compared with the observed data (thin lines) at the Lotung downhole site in the 20 May 1986 aftershock event (No. 8).The re cording at surface is selected as the input motion an d the SH wave is vertically incident from the half-space .
Fig. Bb.The synthetic ground ac celeratio n (thick lines), using Haskell method with the mod ified velocity model where the velocity of the first layer is changed from 120 to 95 m/sec, compared with the observed data (thin lines) at the Lotung downhole site in the 20 May 1986 afte rsho ck event (No. 8).The re cording at surface is selected as the input motion and the SH wave is vertically incide nt from the half-space .

Fig. 9 .
Fig. 9.The synthetic ground acceleration (thick lines), using Haskell method, compared with the observed data (thin line s) at the Lotung do wnhole site in the 20 May 1986 mainshock event (No .7) .The recording at surface is selected as the input motion and the SH wave is vertically incident from the half-space.
Fig. JOa.The synthetic ground acceleration (thick lines), using DESRA-2 nu merical schedule, compared with the observed data (thin lines) at the Lotung downhole site in the 20 May 1986 mainshock event (No. 7).The recording at the 47 m depth is selected as the input motion.
Fig. I Ob.The synthetic ground acceleration (thick lines), using DESRA-2 nu merical schedule, compared with the observed data (thin lines) at the Lotung downhole site in the 20 May 1986 mainshock event (No. 7).The recording at the 17 m depth is selected as the input motion.

Fig. 11 .
Fig. 11.The stress-strain relationship (a) and stress variation with time (b) at the 10 rn depth of the Lotung downhole site in the 20 May 1986 rnainshock event (No. 7).

Table 1 .
The source parameters and PGA values of 30 earthquakes recorded at LSST array site.

Table 2 .
Velocity structure used in Haskell method at the Lotung downhole array.

Table 3 .
Velocity structure used in DESRA-2 scheme at the Lotung downhole array.