Inversion of Love waves in earthquake ground motion records for two-dimensional S-wave velocity model of deep sedimentary layers

We propose a new waveform inversion method to estimate the 2D S-wave velocity structure of deep sedimentary layers using broadband Love waves. As a preprocessing operation in our inversion scheme, we decompose earthquake observation records into velocity waveforms for periods of 1 s each. Then, we include in the inversion only those periods for which the assumption of 2D propagation holds, which we propose to determine through a principal component analysis. A linearized iterative inversion analysis for the selected Love wave segments filtered for periods of 1 s each allows a detailed estimation of the boundary shapes of interfaces over the seismic bedrock with an S-wave velocity of approximately 3 km/s. We demonstrate the effectiveness of the technique with applications to observed seismograms in the Kanto Plain, Japan. The differences between the estimated and existing velocity structure models are remarkable at the basin edges. Our results show remarkable differences from previous existing structural models, particularly near the basin edges while being in good agreement with the surface geology. Since a subsurface structure at a basin edge strongly affects the earthquake ground motions in a basin with the generation of surface waves, our method can provide a detailed model of a complex S-wave velocity structure at an edge part for strong ground motion prediction.


Introduction
Predicting earthquake ground motion using a theoretical approach such as a FD simulation requires an appropriate model of the subsurface S-wave and P-wave velocities, density, and attenuation factors. Particularly, the S-wave velocities of deep sedimentary layers over seismic bedrock with an S-wave velocity of approximately 3 km/s must be accurately modeled when estimating sitespecific ground motions at periods of less than approximately 10 s in a large basin. Accordingly, several 3D velocity models have been constructed for large basins in Japan (e.g., Sato et al. 1999;Yamanaka and Yamada 2006;HERP 2009HERP , 2012HERP , 2017Yoshimoto and Takemura 2014). Of these models, hereinafter, we refer to the models by HERP (2009), HERP (2012, and HERP (2017) as the HERP 2009, HERP 2012, and HERP 2017 models, respectively. The HERP 2009 and HERP 2012 models cover all of Japan. The latter model was constructed from procedures suggested by Koketsu et al. (2012) as an updated version of the HERP 2009 model, with improvements made with additional data such as the results of earthquake ground motion simulations and geophysical surveys. The HERP 2017 model was constructed only for the Kanto Plain by combining S-wave velocity structures deeper than a Open Access *Correspondence: kasamake@kajima.com 1 Kajima Corporation (Former Tokyo Institute of Technology), 2-19-1 Tobitakyu, Chofu-shi, Tokyo, Japan Full list of author information is available at the end of the article layer with an S-wave velocity of 0.5 km/s and shallower than that simultaneously, which had ever been modeled separately. These velocity models have been validated by comparing simulated earthquake ground motions with observed motions for moderate-size seismic events. For example, the HERP 2017 model was examined for its appropriateness on the basis of estimations of goodnessof-fit analyses between synthetic and observed waveforms for several earthquakes.
An S-wave velocity structure at the edge of a sedimentary basin is well known to affect the ground motion characteristics in an entire basin, because it generates/ induces surface waves when affected by incident body waves. For example, Kawase (1996) showed that a large amplitude area was caused by the constructive interference of direct S-waves and basin-induced Rayleigh waves generated at a basin edge during the 1995 Hyogo-ken Nanbu earthquake through strong ground motion simulations. Most of the existing 3D models were mainly established by interpolating accumulations of 1D profiles with several velocity discontinuities from geophysical explorations. The 2D profiles detected by seismic reflection and refraction surveys are useful for constructing a 3D model, but the amount of profile data is limited. Therefore, the existing 3D models may lack detailed features in the edge parts because of difficulties in modeling their complicated shapes.
As one of the alternatives for modeling deep sedimentary layers with a detailed basin-edge structure, inversion analyses using seismograms from moderate-size events have been performed for 3D basin profiles (e.g., Aoi 2002;Iwaki and Iwata 2011). However, these 3D inversions rely on iterative forward calculations for partial derivatives using a 3D FD scheme with numerous computational costs. Thus, a long computational time is required to conduct many inversion analyses using various conditions and initial models. Therefore, a waveform inversion method for retrieving a 2D profile of the deep sedimentary layers from earthquake records can be practically effective because of its low computational costs. As examples of previous studies on 2D waveform inversions, Aoi et al. (1995) proposed an inversion method with an assumed plane incident wave that employs a 2D boundary-element method. Ji et al. (2000) set a 1D structural model outside a basin and performed an inversion analysis to estimate a 2D basin structure by perturbing the depths of several control points in numerical experiments. Most of the studies on 2D waveform inversions have examined the performances of their methods through numerical experiments with synthetic ground motions. Few case studies have applied 2D inversion techniques to observation records (e.g., Zhao et al. 2004;Hikima and Koketsu 2010). One of the reasons for the difficulty in applying them to actual data is the need to approximate the actual 3D wave propagation to the 2D one. For example, Aoi et al. (1997) pointed out that a simple vertically impinged plane wave might not sufficiently represent a real complex incident wave and cause errors in the 2D inversion results. Hikima and Koketsu (2010) assumed that body and surface wave propagations can be fully reproduced by 2D modeling. However, they could not adequately verify the assumption of the 2D approximations of the body and surface wave propagations. The violation of the 2D assumption might distort an estimated 2D model from a true velocity structure.
In this study, we propose a new waveform inversion method to estimate the 2D S-wave velocity structure of deep sedimentary layers using Love waves from a moderate-size earthquake. As a preprocessing operation in our method, we first analyze the observation records for the validation of the 2D Love wave propagation assumption to select the corresponding segments to be used in the inversion. We demonstrate the applicability of our method to actual data recorded in the Kanto Plain, Japan, from comparisons of estimated 2D structures with existing models.

Earthquake records
Here, we describe the characteristics of long-period earthquake ground motion records acquired during a moderate-size M J 6.4 (JMA magnitude) earthquake in the Kanto Plain, Japan, that we will subject later to waveform inversion. Figure 1 shows the locations of the earthquake epicenter estimated by the JMA and seismic observatories, whose records are used in this study. The top depths of the seismic bedrock with an S-wave velocity of 3.2 km/s from the HERP 2009 model are also illustrated in the figure by red contours. The event is a strike-slip crustal earthquake with a focal depth of 14 km and has a fault mechanism with a compression axis in the northwest-southeast direction (Fukuyama et al. 1998). Densely distributed strong-motion stations have been installed as parts of strong-motion networks such as K-NET, KiKnet (NIED 2019), and MeSO-net (Hirata et al. 2009) in the study area. The K-NET stations have accelerometers installed on the surface, while two accelerometers, one on the surface and one at depths greater than 100 m, are used at the KiK-net stations. The accelerometers at the MeSO-net stations were installed at a depth of 20 m. We used the surface records at the K-NET and KiK-net stations and the borehole records at the MeSO-net stations as the surface motions.
The velocity Fourier amplitude spectra at the SIT003, TKY007, and E.KNJM stations are shown in Fig. 2 as examples of the observed ground motion records. These stations are located in areas with thick sedimentary layers (e.g., Shima 1977), which can amplify long-period ground motions with spectral peaks at a period of approximately 8 s. The sediments in the central part of the Kanto Plain have a predominant period of 8 s of motion, as determined in previous studies with a focus on long-period surface wave propagation (e.g., Tanaka et al. 1980;Yamanaka et al. 1989;Kinoshita et al. 1992). Figure 3 shows the velocity waveforms at eight stations in the northwest-southeast direction, which is almost equal to the transverse direction, along Line-1 from TKY003 to E.MZUM. The vertical axis of each panel denotes the distance from TKY003 at the edge of the basin. This distance is hereinafter referred to as the basin-edge distance. The velocity waveforms were calculated from the acceleration records by integration with the fast Fourier transform. A Gaussian bandpass filter represented by Eq. (1) was also applied to the records in the integration at central periods, T m , of 1, 3, 5, and 9 s for the filtered velocities. This filter shape is expressed as where T j represents a discrete period calculated by the Fourier transform. γ is a coefficient for determining the filter width, which is set to 50 with reference to Uetake and Kudo (2001), who investigated the characteristics of long-period Love wave propagations in the Kanto area.
(1) Fig. 1 Locations of the epicenter and strong ground motion observatories used in this study. They are indicated by the focal sphere and squares, respectively. Red contours denote the top depths of the seismic bedrock modeled by HERP (2009). The 2D S-wave velocity structures of the deep sedimentary layers were estimated along Line-1 to Line-3, represented by the thick black lines Figure 3 shows that large amplitude surface waves propagate to the northeast direction, which is almost equivalent to the radial direction, after the arrival of S-waves. For example, the record at E.IMIM at a period of 3 s has the surface wave arriving at approximately 70 s with a much larger amplitude than that of the S-wave. We can also see surface waves with amplitudes larger than those of the S-waves in other periods.

Particle motions
The horizontal orbits of the velocity motions at the eight stations along Line-1 are shown in Fig. 4a. These orbits are calculated for the segments, including the surface waves at the times indicated by the crosses from 0 s in Fig. 3. We use the motions at periods from 1 to 10 s, because the surface wave motions whose generations and propagations are influenced by the deep sedimentary layers over the seismic bedrock are predominant, especially in this period range. The particle motions for periods with a 1 s interval are arranged from the top to the bottom in Fig. 4a. The horizontal axes in the figures denote the basin-edge distances of each station. The principal axis directions, evaluated by principal component analysis (Montalbetti and Kanasewich 1970), are superimposed as red lines. The motions are polarized in the transverse direction, especially at long periods, such as a period of 10 s, which means that the transversely oriented components of the ground motions are mainly Love waves propagating in the radial direction.

Preprocessing
We select the segments used for the waveform inversion by inspecting the polarization features of the horizontal particle motions described in the previous section. The amplitude of the motion in the radial direction becomes large as the basin-edge distance increases, such as in the horizontal motions at SIT003 and E.MZUM at a period Observed velocity waveforms in the transverse direction. These waveforms are bandpass filtered for periods of 1, 3, 5, and 9 s. Triangles and circles indicate the arrival times of the P-waves and S-waves, respectively. The vertical axes denote the basin-edge distances. Segments at the times indicated with crosses from 0 s were used in the calculations for horizontal particle motions in Fig. 4a of 7 s in Fig. 4a. Accordingly, we observe deviations of the principal directions from the transverse direction and increases in the ellipticities of the horizontal particle motions. We also see that the principal axis directions of the particle motions tend to differ largely from the transverse direction for short periods. The short-period ground motions may be greatly influenced by scattering waves due to 3D irregularities. These biases of the principal axis directions of the motions clearly suggest a difficulty in assuming the 2D approximation of the Love wave propagations, especially for short periods. Therefore, bandpass filtered segments of Love wave motions must be selected for our inversion to estimate a 2D S-wave velocity structure. Hereinafter, the assumption of the 2D approximation in the Love wave propagations is referred to as two-dimensionality. From the above analyses of the observed Love waves in each period, we can suppose that the two-dimensionality is satisfied in the propagation of the long-period Love wave motion at stations with long basin-edge distances. The bandpass filtered segments to be used in the 2D inversion are selected from the polarization features (e.g.,

Fig. 4
Results of principal component analyses for records along Line-1. a Horizontal particle motions and principal axis directions, indicated by the black and red lines, respectively. The horizontal axis indicates a basin-edge distance. b E values for the observatories at individual periods. The waveforms at stations located within the distances indicated by the arrows are used in the 2D inversion analysis. c Relationship between periods of ground motions and basin-edge distances for confirmation of the two-dimensionality in 2D inversion Vidale 1986). As described above, the propagations of the short-period Love waves contain non-two-dimensional waves due to 3D irregularities. We thus define an E value to examine the two-dimensionality quantitatively as follows: where θ m i and ∅ m i represent the principal axis direction and the ellipticity of the motions calculated by the principal component analyses using the segments in the transverse and radial directions at the ith station and the mth period, respectively. Here, θ m i is measured from − 90° to 90°. We suppose that the two-dimensionality of the Love wave propagation is sufficiently satisfied at a period of 10 s, because this period is longer than the predominant period of the Love wave amplification expected from the S-wave velocity structure of the deep sedimentary layers, as explained earlier. Thus, the principal axis direction θ i as their product. The reference period can be set to the longest period considering the ratios of the Fourier amplitude spectra between the pretrigger and principal motions of the observation records. We set the reference period to 10 or 9 s in this study. Figure 4b shows the E values of the Love wave segments at the eight stations along Line-1. The horizontal axes in the figure indicate the basin-edge distances. The E values over periods from 1 to 9 s are arranged from the top to the bottom. The E values of the stations for periods of 1 s each are plotted with circles in Fig. 4b. The E values have a clear dependency on the periods. Namely, we see that the long-period segments have small E values at stations far from TKY003. For example, the principal axis direction for a period of 8 s deviates from the transverse direction by approximately 45° at distances greater than 50 km (as shown in Fig. 4a) with an E value of approximately 0.2. Furthermore, the biases of the principal axis directions of the motions for a period of 5 s are also large at distances greater than 15 km (as shown in Fig. 4a) with E values greater than 0.2. We, therefore, considered that small E values less than approximately 0.2 indicate satisfaction of the two-dimensionality assumption for the 2D waveform inversion. The maximum basin-edge distances of the stations with E values less than 0.2 for individual periods are shown by the arrows in Figs. 4a-c. In Fig. 4c, the horizontal and vertical axes denote the basin-edge distance and the period of motion, respectively. We also illustrate (gray-shaded area in the figure) a range of basin-edge distances and periods, where the two-dimensionality can be satisfied. We use the bandpass filtered segments for periods and distances inside this range for the inversion analysis.
The most important point of the above preprocessing is that the number of segments used in the 2D waveform inversion differs depending on the periods of the motions. Table 1 shows the number of stations whose E values are less than 0.2 and the maximum basinedge distances used in later inversion analyses for each period. The shorter the periods of the Love wave motions are, the smaller the number of segments used in the inversion. The short-period Love waves at stations around only the edge part can be used in the 2D inversion, while the long-period Love waves at stations at distances far from the edge can be used.

Model parameterization
In the waveform inversion, we assume a homogeneous layered model with irregular interfaces. The boundary shapes of the layers in the 2D model are deduced in the inversion of the bandpass filtered segments of the Love wave motions selected in the preprocessing. Physical parameters such as the S-wave and P-wave velocities, density, and attenuation factors are given in advance.
We employed the boundary shape parameterization suggested by Aoi et al. (1997). The difference in discretized boundary shapes between the present and initial 2D models, ζ(x) − ζ 0 (x), is described by a linear combination of the coefficients p n,k and basis functions c k as in Eq. (3): where x denotes the horizontal distance from the left end of the model. The suffixes n and k represent the nth boundary interface and the kth basis function, respectively. c k is common to all the interfaces and is represented as follows: The model to be estimated has a length of 2L in the range of − L ≤ x ≤ L. The model in the horizontal direction is equally divided into K + 1 pieces with K + 2 node points that are numbered from 0 to K + 1. Δ indicates the horizontal interval of the basis functions, which is equal to 2L/K + 1. x k denotes the horizontal distance of the kth basis function contributing to the depths at distances from x k-1 to x k+1 . p n,k denotes the unknown parameters in the inversion used to minimize the residuals between the observed and synthetic Love waves, as described in the next subsection.

Inversion Method
We employed a linearized iterative inversion algorithm using the Gauss-Newton method to obtain an optimized solution. The linearized system to be solved can be expressed as where δp, A, G, and e are the correction matrix, the Jacobian matrix, the smoothing matrix, and the residual matrix, respectively. The partial derivatives of the Jacobian matrix A can be numerically computed with model perturbations, dp, using the backward FD scheme. dp is set to 2.5 times the vertical FD grid spacing in the 2.5D forward calculation, which is described in the next subsection.
The residual matrix e is defined from the residuals between the observed and synthetic segments. We evaluate the residuals with the L2 norm as follows: To make the inversion robust, we include the smoothing matrix G, for which the differences between p n,k and p n,k+1 at the nth boundary interface do not change largely, as follows: where L denotes the number of layer interfaces to be deduced. β is a coefficient for determining the strength g = β 2 L · (K + 1) p n,k − p n,k+1 , of the smoothing, which is set to 30 in consideration of the stability of solutions from the results of several inversion analyses using different coefficients in this study. The number of rows of the matrix G is given by multiplying L with K + 1. We solve Eq. (7) using the singular value decomposition method (Lawson and Hanson 1974) to obtain δp for the optimal p n,k values.

Forward FD Modeling
We use a FD staggered-grid formulation of the 2.5D elastic equation of the SH motion (Liner 1991) with the fourth-order approximation in space and the secondorder approximation in time to simulate the Love wave behaviors in the 2D model. The temporal evolutions of the velocity and stress are solved using the explicit FD scheme. The A1 absorbing boundary condition (Clayton and Engquist 1977) and the buffer region (Cerjan et al. 1985) are introduced at the right, left, and bottom edges of our numerical model to reduce the amplitude of nonphysical reflections. The top surface is implemented using the improved vacuum formulation suggested by Zeng et al. (2012). Attenuation is introduced by employing the method of Graves (1996) with a reference frequency of 0.4 Hz (Kasamatsu and Kato 2020). A point load is given at horizontal distances greater than 20 km from the basin edge at the top of the bedrock as an external force to generate the Love waves. One of the major difficulties in waveform inversion is the treatment of effects due to a fault rupture process. We, therefore, not only give the point load far from the basin part but also convolute the observed and synthetic waveforms to solve this problem (e.g., Ji et al. 2000;Amrouche and Yamanaka 2015). The convoluted waveform s m i (t) can be expressed as where T m i is a transfer function at the mth period defined from the Fourier spectra of the synthetic waveforms at the ith station, C m i , and a reference station, C m r , in Eq. (11). An asterisk indicates a complex conjugate of the Fourier spectrum. T m i is multiplied by the observed Fourier spectrum of the reference station, O m r , to generate the convoluted waveform, s m i (t) , using the inverse Fourier transform, F −1 , as shown in Eq. (10).

Inversion result
We applied our method to the records of the Mt. Fuji region earthquake in 2011 (M J 6.4) to retrieve the 2D Fig. 1. Here, we explain the inversion results for Line-1 in detail.

S-wave velocity profiles of the deep sedimentary layers for the three survey lines shown in
We prepared an initial model by HERP (2009) with seven homogenous layers and S-wave velocities of 0.5 to 4.6 km/s. The sediments over the seismic bedrock consist of three layers with S-wave velocities of 0.5, 0.9, and 1.5 km/s. The optimal bottom depth shapes of these three layers were detected from the inversion of the filtered segments of the Love wave motions selected in the previous section (see . Since the HERP 2009 model is an old version of the HERP 2012 model, we used the HERP 2009 model as the initial model in the inversion analysis. We then compared our 2D models with the HERP 2012 model, which includes improvements with additional data, such as the results of earthquake ground motion simulations and geophysical surveys obtained after the release of the HERP 2009 model. We estimated a velocity structure at horizontal distances from 0 to 57 km using the 15 basis functions. The horizontal interval of the basis functions, Δ, represented in Eqs. (4) and (5) was approximately 4 km. A total of 45 p n,k values for the three boundaries of the sediments expressed by Eq. (3) were optimized. We used TKY003 as the reference station for the convolution processing. The grid spacings of the 2.5D FD forward modeling were 50 m in the horizontal direction and 20 m in the vertical direction. The vertical grid spacing was 50 m at depths greater than 5 km to reduce the computational costs. Figure 5 shows the inverted and observed velocity waveforms over periods of 1, 3, 5, and 9 s. The amplitudes of the synthetic motions for a period of 9 s are slightly underestimated at E.MZUM, E.IMIM, and SIT009. However, the overall characteristics over all periods are well reproduced at every station. Even the observed and inverted motions for relatively short periods, such as 1 and 3 s, show high similarities at SIT012, E.KSRM, and E.OKDM. The synthetic waveforms in the initial 2D model (HERP 2009) were calculated by a forward simulation, as shown in Fig. 5. These synthetic waveforms for a period of 9 s contain Love waves with arrival times slightly later than those observed at distances greater than 47 km. Furthermore, the amplitudes of the synthetic motions over periods of 3 and 5 s in the initial model are overestimated at E.OKDM and E.KSRM. Figure 6 shows the estimated 2D structure along Line-1. The physical parameters of each layer are listed in Table 2. The third layer can be regarded as an outcrop at distances from 0 to 5 km. The bottom depth of the first layer becomes deep at distances greater than 5 km. Then, it reaches 0.2 km at distances of approximately 10 km. The bottom depth of the first layer becomes shallow at distances from 12 to 25 km, and the second layer exists near the surface at approximately 17 km.
The initial model (HERP 2009) and the HERP 2012 model are also shown in Fig. 6. At distances greater than 30 km, all of our estimated boundary depths of the sedimentary layers are shallower than those in the HERP 2012 model. The differences between the inverted and HERP 2012 models are remarkable in the shallow parts of the basin-edge area. The estimated bottom depth of the second layer becomes shallow at distances less than 10 km and deep at distances greater than 10 km. As a result, the estimated bottom interface of the second layer becomes steeper than that in the HERP 2012 model at distances from 5 to 15 km. The bottom depth of the first layer becomes slightly shallow around SIT009 in the HERP 2012 model, but the top of the second layer does not exist near the surface, unlike our results.

Fig. 5
Comparison between observed and calculated velocity waveforms along Line-1. Black, red, and blue dotted lines represent the observed, inverted, and synthetic waveforms, respectively. The synthetic waveforms were generated from the initial model by HERP (2009) Fig. 6 Comparison of the inverted 2D velocity structure with previous models. Red, gray dotted, and black lines represent the estimated structure and structural models by HERP (2009HERP ( , 2012 along Line-1, respectively. The horizontal axis indicates the basin-edge distance. Triangles denote the observatory locations Table 2 Physical parameters of layers 1-4 in Figs. 6, 7, 8a As above, our obtained 2D profile shows more irregular structural features around the basin-edge area than does the HERP 2012 model. This outcome suggests that we could model the detailed velocity structure, especially in the edge part, due to the introduction of the short-period Love waves in the inversion.

Comparison with the HERP 2017 Model
For the Kanto Plain, the HERP 2017 model was constructed as a further updated version of the HERP 2012 model. The HERP 2017 model contains layers with S-wave velocities less than 0.5 km/s, which are not included in the HERP 2009 and HERP 2012 models. Our inversion results, therefore, cannot be directly compared with those of the HERP 2017 model. However, we can still find a similarity between them around the basin-edge part along Line-1. Figure 7 compares their 2D structural models at distances from 0 to 22 km as the edge part. The top depths of the bedrock in the HERP 2017 model and our model show high similarity. On the other hand, our estimated upper depths of the second and third layers with S-wave velocities of 0.9 km/s and 1.5 km/s are different from those of the HERP 2017 model.

Comparison with surface geological classifications
The surface geological classifications (Geological Survey of Japan 2015) around the basin-edge area along Line-1 are indicated in Fig. 7. According to a previous study on the relationship between S-wave velocities and geological ages (Yamamizu et al. 1977), those layers with S-wave velocities of 1.2 to 1.6 and 0.4 to 0.9 km/s belong to the Neogene and Quaternary ages, respectively.
In the estimated structure, the third layer with an S-wave velocity of 1.5 km/s can be regarded as an outcrop at distances from 0 to 5 km because of the thin thicknesses of the overlying layers. The surface geology in this area belongs mostly to the Neogene age. At distances from 5 to 13 km, the estimated bottom depth of the first layer with an S-wave velocity of 0.5 km/s gradually becomes deep away from station SIT012. The geological age of the surface formation in this area is categorized as Quaternary with lower terrace sediments. In the estimated structure, a second layer with an S-wave velocity of 0.9 km/s exists near the surface locally at distances of approximately 17 km. Although the surface geology around this part also belongs to the Quaternary age, the Geological Survey of Japan (2015) reported the existence of middle terrace sediments that are older than the lower terrace sediments in this part. Therefore, our estimated S-wave velocity structure near the surface is in good agreement with the surface geology.

Results of lines-2 and 3
Line-2 Figure 8a illustrates an inverted 2D velocity structure along Line-2. The horizontal length along Line-2 is approximately 58 km from the basin edge. We used 28 stations along this line in the inversion analysis. Segments of the Love waves used in the inversion were selected for  Table 1. We excluded the segments for a period of 6 s in the inversion, because the Fourier amplitude spectrum of the motion at each station has a trough during this period. The horizontal axis in Fig. 8a indicates the horizontal distance from the OK.NKYM reference station, where seismic bedrock exists near the surface. Each estimated boundary depth is shallower than that in the initial model (HERP 2009) at almost all distances. These boundary depths are also shallower than those in the HERP 2012 model at distances greater than 25 km. We found remarkable differences between the estimated and initial structures in the basin-edge part. The estimated bottom depth of the first layer is shallower than that in the initial model at approximately 10 km, and then it gradually deepens away from OK.NKYM. The shape of the bottom interface of the first layer in the edge part is similar to that in the HERP 2012 model.

Line-3
Figure 8b shows an inverted 2D velocity structure along Line-3. The horizontal length along Line-3 is approximately 34 km. We used seven records covering periods from 1 to 9 s except for 6 s in the inversion analysis, as shown in Table 1. The segments for a period of 6 s were excluded for the same reason as for the Line-2 inversion, as explained in the previous subsection. All the boundary depths of the sedimentary layers are significantly shallower than those in the initial model (HERP 2009) at distances less than 30 km, and they are in good agreement with the HERP 2012 model. A second layer with an S-wave velocity of 0.9 km/s exists near the surface at distances from 15 to 20 km. The S-wave velocities of the top layers in the deep sediments around this area are reportedly 0.7 to 0.8 km/s (Yamanaka and Yamada 2006). Although we do not consider lateral variation in the S-wave velocities, the existence of the second layer with an S-wave velocity of 0.9 km/s near the surface means that our estimated structure has features similar to those of the surface S-wave variation of the model by Yamanaka and Yamada (2006). The estimated bottom depths of layers 1 to 3 at distances greater than 30 km in our model are almost the same as those of the HERP 2009 and HERP 2012 models.

Effect of introducing short-period motions
We indicated the importance of the preprocessing operation for confirming the two-dimensionality of the Love wave propagations before the 2D waveform inversion with its applications to the earthquake records in the  Fig. 6 previous sections. As a result of introducing the preprocessing step, we were able to obtain detailed 2D profiles, especially at the edge parts, with the introduction of short-period Love wave motions as well as long-period motions. In this section, a numerical experiment using synthetic waveforms generated by a 3D computation is performed to examine the effectiveness of introducing short-period motions in the waveform inversion.
The synthetic records were generated using the 3D model by HERP (2012) for the numerical test. The ground motions of the Mt. Fuji region earthquake in 2011 (M J 6.4) were simulated at the seven stations along Line-3 by employing a 3D FD scheme (e.g., Graves 1996;Pitarka 1999). We regard them as synthetic earthquake ground motion records. These seven stations are located at the same locations as the actual observatories. A double-couple point force with a source mechanism by F-net (Fukuyama et al. 1998) was given as an external force at the hypocenter location. Fault parameters such as the seismic moment are listed in Fig. 1. The HERP 2012 model used in the 3D calculation was discretized with a horizontal grid spacing of 100 m and a vertical grid spacing of 10 m.
We performed two inversion analyses to model sedimentary layers from 0 to 34 km along Line-3 using the synthetic records covering periods from 1 to 9 s and from 7 to 9 s, respectively. Then, we compared the 2D profiles retrieved from the two inversions with the above true profiles (HERP 2012). In the former inversion, the segments for periods from 1 to 9 s selected from the E values were inverted in the same way as the inversion carried out using actual observation records in the previous section. In the latter inversion, only the segments covering periods of 7 to 9 s were used at distances from 0 to 34 km without consideration of the E values. Other conditions in the two inversions were the same as those in the inversion described in the previous section.
The 2D profile from the inversion using the segments corresponding to periods from 1 to 9 s is illustrated in Fig. 9 with the true structure (HERP 2012). The true model is well reconstructed by the inversion. The 2D velocity structure estimated using the segments corresponding to periods from 7 to 9 s is also shown in Fig. 9. Obvious differences are observed between the estimated and true structures at the basin-edge part. The bottom depth of the first layer estimated using only the longperiod segments is shallower than that of the true structure at distances from 0 to 15 km. Furthermore, the upper depths of the second and third layers in the estimated structure using the segments corresponding to periods from 7 to 9 s also become shallow at distances of approximately 5 km. Thus, inversion using only long-period Love waves cannot offer an accurate 2D profile at the basin edge. The preprocessing carried out with the E values to include the short-period segments in the 2D inversion enhances the reliability of the inverted 2D profile.

Conclusions and discussion
We have proposed a waveform inversion method to estimate the 2D S-wave velocity structure of deep sedimentary layers using broadband Love waves. As a preprocessing operation in our inversion scheme, we decompose earthquake observation records into velocity waveforms for periods of 1 s each to confirm the assumption of 2D propagations of the Love waves with a principal component analysis of the horizontal ground motions. Our linearized iterative inversion analysis allows an accurate estimation of boundary shapes from the top depth of seismic bedrock to the bottom depth of a layer with an S-wave velocity of approximately 0.5 km/s from the segments of the Love wave motions for periods of 1 s each, which is selected in the preprocessing step.
We have demonstrated the effectiveness of our technique with applications to observed seismograms in the Kanto Plain, Japan. Remarkable differences were observed between our inverted model and the existing ones (HERP 2009(HERP , 2012, especially in the basin-edge area. Our method can offer detailed 2D velocity profiles at the edges due to the introduction of short-period segments selected from the E values for the inversions. A velocity structure at a basin edge strongly affects the earthquake ground motions in a basin because of the generation of surface waves (e.g., Vidale and Helmberger 1988;Kawase 1996;Hartzell et al. 2016). 2D profiles including the details of the basin-edge structure along many lines from our proposed inversions can be used to Fig. 9 Comparison of 2D velocity structures in a numerical test. Gray line represents the true structure (HERP 2012). Black solid and dotted lines indicate inverted structures from bandpass filtered segments for periods from 1 to 9 s and from 7 to 9 s, respectively improve the existing 3D models for accurate earthquake ground motion evaluations.
In actual work for the prediction of strong ground motion from a large earthquake, many inversion analyses using various conditions are required to estimate the S-wave velocity structures. Since a high-performance computer for these inversions is not always available in practice, we have to reduce the calculation cost in the forward calculation as much as possible. In this study, our 2D inversion employing the 2.5D FD calculations took approximately one day to obtain an optimal solution after performing parallel computing using the openMP algorithm with a conventional PC equipped with four Intel Xeon processors (E7-8891 v4, 2.8 GHz, 10 cores). The calculation time required for one iteration was approximately two hours when using 40 cores in total, and the number of iterations until the solution converged was approximately ten. On the other hand, the 3D FD calculation performed to generate the synthetic records for the numerical test in the previous section took several days with the same conventional PC. This outcome indicates several-100-fold differences in the calculation costs between the 2D and 3D inversion analyses. Although the computational costs depend on the performance of the PC used in the calculation, an optimal solution is still difficult to obtain from a 3D waveform inversion at a reasonable computational cost with a conventional PC.
The above differences in the computational cost can be critical when conducting many inversions of records at dense seismic stations in a large basin. Recently, many earthquake observatories were installed in large basins, especially in seismically active areas. Our 2D inversion technique can be effectively applied to earthquake records along many lines in these areas with reasonable computational costs.