Proceedings of the International Symposium on Engineering Lessons Learned fromthe 2011 Great East Japan Earthquake, March 1-4, 2012, Tokyo, Japan

INVERSION OF EQUIVALENT LINEAR SOIL
PARAMETERS DURING THE 2011 TOHOKU
EARTHQUAKE, JAPANFlorent DE MARTIN1, Hiroshi KAWASE2 and Fabian BONILLA3
1Researcher, BRGM, Division of Risks, Seismic Hazard Team, Orléans, Francef.demartin@brgm.fr
2Professor, Disaster Prevention Research Institute, Safety Control of Urban Space, Disaster
Management for Safe and Secure Society, Kyoto University, Kyoto, Japan
kawase@zeisei.dpri.kyoto-u.ac.jp
3Researcher, IFSTTAR, Paris, Francefabian.bonilla@ifsttar.fr
ABSTRACT: This paper presents the inversion of linear and equivalent linear soil
profile (i.e., S-wave velocity and damping factor) at the KiK-net station MYGH10. Thisstation has been selected among other KiK-net stations thanks to a proxy stress-straintechnique showing that nonlinear behavior of the soil structure under this station may bepossible during strong ground shaking. The inversions, performed by genetic algorithm
coupled with Thomson-Haskell propagator method, show that the entire soil column
might have behaved nonlinearly during the Great East Japan earthquake.
Key Words: Great East Japan earthquake, borehole data inversion, equivalent linear soilcolumnINTRODUCTION
On the afternoon of March 11th, 2011, an undersea megathrust earthquake of magnitude Mw 9
occurred off the Pacific coast of Tohoku, Japan, with an epicenter approximately 70 kilometers east of
the Oshika Peninsula of Tohoku district. Thanks to various Japanese networks, excellent quality data
have been recorded. Particularly, many strong-motion data have been recorded with long durations andhigh peak ground accelerations and velocities (PGAs and PGVs, respectively) in Miyagi, Ibaraki,Tochigi, and Chiba Prefectures. Such long durations and high PGAs are probably characteristics ofrecords lengthened and amplified by site effects (i.e., modification of the waves coming from thesource due to the impedance contrast between the bedrock and the soft sediment basins). Owing to thelarge PGVs observed, the strain level sustained by the soft sediment is large and the response of thesediments is no longer linear. Soil nonlinearity is, however, a very complex phenomena and researchon this topic is essential for the future seismic design input. Thanks to the data recorded during the2011 Tohoku earthquake, different types of nonlinear soil behavior have been observed ranging fromclassic high-frequency de-amplification to liquefaction.
In this paper, we tried to determine the one-dimensional “equivalent soil profile” (i.e., a soil330

profile having their S-wave and damping factors inverted to match the observed nonlinear response interms of shear modulus reduction and damping ratio curves) during the main shock of the Great East
Japan earthquake at some KiK-net stations. First, among a large number of KiK-net stations, weidentify using a stress-strain proxy technique (Idriss, 2011), which KiK-net’s soil structures are
susceptible to behave nonlinearly. Then, among the identified sites, we check the linear soil structureusing weak motion and invert the equivalent linear soil structure during the main shock by using a
genetic algorithm inversion technique (De Martin et al. 2010).

KIK-NET SITES IDENTIFICATION
In order to select KiK-net sites that may behave nonlinearly during strong ground motions, we use thestress-strain proxy technique introduced by Idriss, 2011. In this technique, the stress-strain relationshipis roughly approximated by the PGA-PGV/Vs30 relation where PGA is proportional to the stress andPGV/Vs30 is a proxy for the strain. This technique is well suited to indicate which sites may behave
nonlinearly.
Fig. 1 shows the proxy stress-strain relation for different values of Vs30 at KiK-net sites using
earthquake database from September 2000 to April 2011. The proxy technique shows that low Vs30
sites may behave nonlinearly whereas high Vs30 site may not. Among these sites, we select the site
MYGH10 for which its proxy stress-strain relation is shown in Fig. 2.

Fig. 1 Proxy stress-strain relation for different value of Vs30. “Shear Strain” is approximated by
PGV/Vs30. The Vs30 binning is the same as Idriss (2011).

331

Fig. 2 Proxy stress-strain relation for MYGH10.

Fig. 3 Initial velocity profile of MYGH10 (left)
provided by KiK-net.

MYGH10 VELOCITY PROFILE AND DATABASE
The initial velocity profile of MYGH10 provided by KiK-net website is shown in Fig. 3. The
shear-wave velocity Vs range from 110 m/s at the free surface to 770 m/s at the downhole (205 m).
The initial soil column used for inversions at MYGH10 is presented in Table 1. Ground motions used
for inversions are listed in Table 2.

Table 1 Soil column at MYGH10. S-wave factors are used to determine the search space of genetic
algorithm inversions of weak motions data.
N°
1
2
3
4
5Thickness
(m)
12
31
80
91Vs(m/s)
110
250
390
590
770332

Density(kg/m3)
1750
1800
1850
19502030

S-wave factors
for search space
0.5 – 1.50.5 – 1.5
0.5 – 1.5
0.5 – 1.5
0.5 – 1.5

Table 2 Ground motions used for inversions
N°
EQ1
EQ2
EQ3
EQ4
EQ5
EQ6

Date
YY.MM.DD/HH.MM
11.03.11/14.46
09.09.22/20.40
09.06.01/00.33
08.12.04/17.29
08.10.30/00.48
06.03.29/19.28

Lon.
(°)
142.8610
141.6645
141.6313
141.4370
141.7323
141.5482

Lat.
(°)
38.1035
37.5968
37.777837.5140
38.046237.2018

Mag.
Mw
9.0
4.7
4.6
4.4
5.1
4.8

PGANS
(cm/s2)8712727143732PGAEW(cm/s2)853
2329153634

PGAUD(cm/s2)
622
111351611OBJECTIVE FUNCTION AND S-WAVE WINDOW SELECTION
As shown by De Martin (2011), the use of a frequency-domain objective function for inverse problems
at borehole stations can lead to biased results. In practice, the length of the cosine-shaped windows
used to isolate the S-wave portion is fixed as a percentage of the length of the S-wave portion.
Common values range (as a rule of thumb) between 15% and 50% and one can show that for such arange, the spectral ratios may strongly change: the location of the resonant peaks along the frequency
axis can be slightly affected and the amplitude of the peaks can be greatly affected. As a result, theinversion of the S-wave velocity and damping factors can be distorted because of an arbitrary
application of the cosine-shaped window.
Therefore, the inversions performed in this paper use a time-domain objective function which isthe summation of the square residual between observed and simulated acceleration time histories on
the S-wave window. By using this technique, the objective function is computed on raw time history
data (or filtered time history data). In our case, we use time history band-pass filtered in the frequency
range [0.1 – 5.0] Hz.
The S-wave portion used as objective function is determined by examination of the particle
motion on acceleration time history in the [transverse – up down] and [radial – up down] planes. An
example of particle motion for EQ2 is shown in Fig. 4.

GENETIC ALGORITHM INVERSIONSInversions’ information
For each earthquake listed in Table 2, we perform eight independent genetic algorithm inversions inorder to compute a mean and standard deviation of the inverted parameters. Each genetic algorithm
inversion runs overs 11,000 populations. A Monte-Carlo search is performed on the first 1000
populations composed by 2048 individuals (one individual represents the velocity and damping soil
structure from free surface to downhole). Then, the genetic algorithm inversion runs over 10,000
populations composed by 1024 individuals. Consequently, for an earthquake, a total of 12,288,000 soilcolumns are compared with the observation. Such an inversion runs for about 27 h 41 min on an AMD
Opteron @ 2.3 GHz. For a single inversion, elitism is activated so that the 10 best soil columns arealways present from one population to another. Selection is performed by tournament and theprobability of reproduction is 85%. The probability of mutation is 0.1%.

333Linear soil structure
The linear soil structure is obtained by inversions EQ2 to EQ6. The search space of the S-wave
velocity is, for each layer of MYGH10, [0.5 – 1.5] times the S-wave KiK-net initial velocity. Theinterval [0.5 – 1.5] is decomposed in 16 equally spaced values. The damping factor of each layer issearched in the interval [0 – 5] % decomposed in 16 equally spaced values as well. As a result, a totalof 1.1E+12 solutions are possible. Only 1.1E-3% of the solutions is computed during a genetic
algorithm inversion. The convergence curves for eight independent inversions are shown in Fig. 5. At
the end of the inversion, one unique optimum is found.
Fig. 6 shows that for EQ2 the inversion improves the time history simulation compared to theinitial KiK-net velocity profile. A mean and standard deviation over all the best soil profiles of each
inversion computed for the EQ2 to EQ6 are shown in Fig. 7. We see that except for the first two layers,
the standard deviation includes the KiK-net velocity profile.

Fig. 4 Acceleration particle motion in the [transverse – up down] plane at downhole for EQ2 between
25 and 32 seconds (left to right) and in the frequency range from [0.1-1.0] to [0.1-10.0] Hz (top tobottom). S-wave arrival between 27 and 28 s, within the red box, shows a clear horizontal motion in
the transverse direction.

334

Fig. 5 Convergence curves of eight independent inversions for EQ2.

Fig. 6 Comparison between transverse time histories for EQ2 for MYGH10 in the frequency band [0.1– 5.0] Hz: Observations (black), KiK-net velocity profile simulation (red) and inversed velocity profile
simulation (blue).

Equivalent linear soil structure
The mean soil column found from the inversions of EQ2 to EQ6 is used as starting point for inverting
the main shock (EQ1). The search space for the S-wave velocity is [0.5 – 1.0] times the mean soilcolumn of Fig. 7. The damping is search in the interval [0.0 – 50.0] %. Over 8 independent inversions,a common single optimum soil profile is found.
The observed S-wave window used for the computation of the objective function is shown in Fig.
8 together with the observed and the time history computed via the inversed velocity profile. We cansee that the S-wave window is well reproduced in term of frequency content, but could be improved in335

terms of amplitude. The inverted soil profile is shown in Fig. 9 and the shear modulus reduction withthe strain is shown in Fig. 10. According to the inversion, all soil layers sustain a shear-modulus
reduction but only the first three layers have an increase of damping.

Fig. 7 Mean and standard deviation computed over all best soil profile of all inversions for EQ2 toEQ6. Left panel: S-wave velocity profile from 0 to -205 m depth. Middle panel: Zoom from 0 to -5 m
depth. Right panel: Damping factors from 0 to -205 m depth.

Fig. 8 Comparison between observed and simulated transverse time histories for the main shock(EQ1).

336

Fig. 9 Mean and standard deviation computed by inversion for EQ1. Left panel: S-wave velocity
profile from 0 to -205 m depth. Middle panel: Zoom from 0 to -5 m depth. Right panel: Damping
factors from 0 to -205 m depth.

Fig. 10 Shear modulus degradation in layer 1 (L1) to layer 5 (L5). G0 is computed using the mean
inverted velocity profile shown in Fig. 7.

CONCLUSIONSThis study shows preliminary results for the inversion of the linear and equivalent linear soil profilesof the station MYGH10. This station has been selected among other KiK-net stations thanks to a proxy337

stress-strain technique showing that nonlinear behavior of the soil structure under this station may bepossible.
After the inversions of the linear and equivalent linear soil parameters (i.e., S-wave velocity anddamping factors), the match with the observed time histories is fair, but it could be improved byincreasing the number of layer of the initial soil profile to give more freedom to the inversions.

ACKNOWLEDGMENTS
This work was supported by the ONAMAZU project (JST/ANR Flash Japan). The authors arethankful to KiK-net for proving earthquake data and velocity profile on the Internet. F. DE MARTIN
thanks A. DUCELLIER for her precious help and advices.

REFERENCESDe Martin, F., Kawase, H., and Modaressi-Farahmand Razavi, A. (2010). “Nonlinear Soil Response ofa Borehole Station Based on One-Dimensional Inversion During the 2005 West off Fukuoka
Prefecture Earthquake”. Bull. Seism. Soc. Am. Vol 10, No 1, pp 151-171
De Martin, F. (2011), “Comparison of frequency and time-domain objective functions for boreholestation’s inverse problems”, proc. 5th international conference on earthquake and geotechnicalengineering, Santigo, Chile, 2011.
Idriss, I. M. (2011). “Use of Vs30 to represent local site conditions”, 4th IASPEI / IAEE International
Symposium Effects of Surface Geology on Strong Ground Motion, Santa Barbara, California August26, 2011338