Shallow subsurface structure estimated from dense aftershock records and microtremor observations in Furukawa district, Miyagi, Japan

Key words: array observation, H/V spectral ratio, microtremor survey, phase velocity, S-wave velocity structure, shallow subsurface layers, single-site observation, transfer function.

Severe structural damage due to ground motion occurred in some limited areas of the Tohoku and Kanto regions, north-east Japan, during the Great Eastern Japan Earthquake of 2011, off the Pacific coast of Tohoku. The Furukawa district in Miyagi prefecture was one of the most significantly damaged areas (Goto and Morikawa, 2012). The frequency content of ground motion records from the K-NET MYG006 and JMA Furukawa stations, which are located in the Furukawa district, is similar to the frequency content in JMA Kobe and JR Takatori records from the 1995 Kobe earthquake. This helps to explain why ground motion in the Furukawa district was so structurally damaging (Goto and Morikawa, 2012; Hayashi et al., 2013).

Figure 1 shows the location of the Furukawa district, and the area of significant damage investigated by Goto and Morikawa (2012). The structural damage recorded by the Osaki city office, and liquefactions, were observed locally to have occurred within an area of ~1.0 × 0.5 km², in the centre of the downtown area. In order to investigate the reason for the damage being concentrated in that limited area, Goto et al. (2012) focused on the ground motion differences across the Furukawa district, and established a temporary network of seismometers in the downtown area, named the Furukawa Seismometer Network (FuSeN). FuSeN consists of more than 30 accelerometers with a spatial interval of ~100 m, and is one of the densest seismometer networks in the world. The observed peak ground acceleration (PGA) and velocity (PGV) indicate that ground motion was greatly amplified in the significantly damaged area (Goto et al., 2012). Thus, we had a major question: why does this spatial difference of ground motion amplification appear in the Furukawa district?

Fig. 1.  Location of the Furukawa district, and the area damaged during the Great Eastern Japan Earthquake of 2011. Red symbols indicate the location of severe structural damage recorded by the Osaki city office.

The relationship between shallow subsurface geological structure and site amplification of ground motion has been widely discussed in the field of earthquake engineering. As one example, seismic microzonation (e.g. Borcherdt, 1970; Lachet et al., 1996) aims to evaluate the seismic risk at a particular site on the basis of the available dataset for geological structures; e.g. soil profile, logging data, geomorphological data, surface geology, geophysical exploration results, etc. Although the dataset is required to be densely distributed for this objective, the available soil profiles and logs may be relatively sparse due to their cost. Microtremor observation is one of the most efficient techniques in geophysical exploration, and has been widely adopted to add to the available dataset (e.g. Lermo and Chávez-García, 1994; Bour et al., 1998; Tuladhar et al., 2004; Panou et al., 2005).

Inatani et al. (2013) estimated the shallow ground structure in the Furukawa district on the basis of the ground motion records obtained by FuSeN, and soil profiles. Thicker surface soil was estimated at the sites with the larger amplifications. Inatani et al. (2013) used relative transfer functions obtained from the ground motion records, but these constrained only the relative differences between velocity structures. In addition, the resolution was limited to the scale of spatial intervals between the sensors in FuSeN.

In the project reported here, single-site observations and array observations of microtremors were conducted to revise the subsurface structure in the Furukawa district. We carried out array observations of microtremors to identify the S-wave profile, and extremely dense single-site observations of microtremors in order to obtain the densely distributed dataset. The two sets of results were merged to model the shallow (~100 m) subsurface structure. We then revised the subsurface structure in the Furukawa district, and discuss the contribution of the geological conditions to the site amplification of ground motion.

Array observations of microtremors can be used to estimate phase velocities of surface waves. The pioneering research for the estimation of phase velocities was performed by Aki (1957) and Capon (1969): the former technique is called the spatial auto-correlation (SPAC) method and the latter is the frequency-wave number (F-K) method.

In order to estimate the shallow ground structure in the area, array observations of microtremors were carried out at three sites, labelled MMK, SXZ, and F4S. The locations of these sites are shown in Figure 2. MMK is in a park named Mikkamachi, where a strong motion observation station of the Japan Meteorological Agency (JMA) is located; this site is within the area severely damaged during the 2011 Tohoku earthquake, as shown in Figure 1. The arrays consist of two equilateral triangles of different sizes set up using seven sensors: one array has side lengths of 5 and 15 m, and the other 15 and 45 m. SXZ is in a parking lot of the Memorial Museum of Sakuzo Yoshino, located in the north-eastern area of Furukama where little damage occurred during the earthquake. The array configurations at SXZ are the same as at MMK. F4S is in the athletic field of Furukawa Fourth Elementary School, located north of the downtown area, and here also little damage occurred during the earthquake. This array is configured with 5 and 15 m side lengths, using seven sensors.

Fig. 2.  Target area for the single-site and array observations of microtremor. The triangles indicate the location of Furukawa Seismometer Network (FuSeN) stations (Goto et al., 2012). The circles and stars show the sites for strong motion observation operated by the government and array observation sites for this work, respectively.

The observations were performed over three days from 14 to 16 December 2014. It snowed on the 16 December, thus, we were unable to collect the data for the largest array (45 m side length) at F4S. We made simultaneous observations with seven sensors at each site. Each sensor was a three-component moving-coil-type velocity sensor with a natural period T₀ = 0.5 s. The recording system includes amplifier, digitiser, and global positioning system (GPS) clock in one chassis. The signals were amplified 256 times using an amplifier with a low pass filter with a cut-off at 30 Hz, and then recorded in digital form with 0.005 s sampling interval and 24 bit quantisation. The sensors were synchronised using the GPS timing, with errors less than 1 ms. This system can also record circuit noise and reduce the self-noise by post processing. Thus, the system displays a dynamic range of more than 140 dB with a very low noise floor (Araki et al., 2011).

Each dataset was acquired for 30 min. We obtained two datasets for MMK and SXZ, and one dataset for F4S. The observed records were contaminated by artificial noise such as the vibration from vehicles, trains, persons walking, and so on. We selected good-quality portions of vertical component records, 20.48 s in length, from each dataset. Corrections were applied for instrumental characteristics, which were determined from the step responses of the sensor pendulums. The step responses were recorded at least 10 times at each site, and we searched the optimal natural periods and damping factors for theoretical impulse response functions to fit the observed responses. Power and cross spectra between two sites in an array were calculated using the fast Fourier transform (FFT) technique, and the SPAC coefficients were obtained (Aki, 1957). Optimal values of the phase velocities of Rayleigh waves were determined to fit the observed values of SPAC coefficients by means of a grid search across all of the discretised frequencies. We assume here that the vertical component of microtremor motion consists mainly of Rayleigh-wave fundamental modes. Some researchers have pointed out that microtremors may be contaminated by higher mode Rayleigh waves, which affects the accuracy of the estimated velocity profile (e.g. Tokimatsu et al., 1992; Bonnefoy-Claudet et al., 2006). In the case where S-wave velocity monotonically increases with depth, the contribution from higher modes is limited. Therefore, we constrain the velocity profile so that velocities monotonically increase with depth.

The H/V spectral ratios were also calculated, from three-component data, following the same procedure to calculate the power spectra for horizontal and vertical components. In general, microtremors consists of Love and Rayleigh waves and their fundamental and higher modes. The H/V spectral ratio would be controlled by all the component waves (e.g. Arai and Tokimatsu, 2004, 2005; Bonnefoy-Claudet et al., 2006). A detailed discussion would be required to link the H/V ratio to the velocity profile. On the other hand, Arai and Tokimatsu (2000) suggested that the frequency variation, and the peak periods of the H/V spectral ratio, might be explained by the ellipticity of the fundamental Rayleigh-wave mode. We here assume the peak periods can be directly compared to the ellipticity of the fundamental Rayleigh-wave mode component of the signal.

We obtain velocity profiles by forward modelling at F4S, SXZ and MMK to fit the estimated phase velocity in a frequency range of 2.0–10.0 Hz. 1D stratified media with six layers are assumed. S-wave and P-wave velocities, and density, excepting the S-wave velocity of the first layer, are constrained to be common to all the sites. S-wave and P-wave velocities, and a density profile, to 20 m depth are available at the K-NET MYG006 station, and we refer to this profile for initial modelling of the shallower part (first to third layers) of the section. As discussed later, the estimated parameter profiles at the microtremor sites are not exactly equal to the profile at the K-NET MYG006 station. We focus on better fitting between the estimated and calculated phase velocities. At the same time, peak periods of H/V spectral ratios are compared to the calculated ellipticity for the fundamental mode of the Rayleigh wave. Phase velocity and ellipticity for the Rayleigh-wave fundamental mode are calculated to solve an eigenvalue problem for the displacement-stress vector (see section 7.2 in Aki and Richards, 2002).

The phase velocities and H/V spectral ratios are shown in Figures 3 and 4, respectively. Only structures shallower than ~100 m can be estimated from the phase velocities, because the frequency range of the estimated phase velocities is limited. However, the ranges of the H/V observations extend lower frequencies, around 0.2 Hz, and the observations imply that the deeper structures that correspond to the peaks around 0.2 Hz are common to all sites. Table 1 shows the estimated velocity structures at each site. Structure deeper than 500 m depends on the peaks of the H/V spectral ratio around 0.2 Hz, and the same model is considered appropriate for all three sites.

Fig. 3.  Estimated phase velocities at each site: (a) F4S, (b) SXZ and (c) MMK. Different colours denote the different lengths of the legs of an array. The solid lines show the theoretical dispersion curves obtained from the estimated S-wave velocity structures.

Fig. 4.  H/V spectral ratio (red) and theoretical ellipticity for the Rayleigh-wave fundamental mode (green) at each site: (a) F4S, (b) SXZ and (c) MMK.

Table 1.  Estimated velocity structures from array observation of microtremor at F4S, SXZ, and MMK sites.

Kanai and Tanaka (1961) revealed that frequency contents of microtremors are well correlated with the soil types. Nakamura (1989) proposed using the Fourier spectral ratio of horizontal to vertical components, widely referred to as the H/V spectral ratio of the microtremor, and he claimed that the ratio is related to the S-wave transfer function at the ground surface. Although the claim is not exactly correct because the microtremor mainly consists of surface waves (e.g. Akamatsu, 1961), the period (frequency) that gives the peak value of the spectral ratio is well correlated with the natural period of the subsurface ground (e.g. Lachet and Bard, 1994). Tokimatsu and Miyadera (1992) found that the peak period of the H/V spectral ratio is explained by the theoretical H/V amplitude ratio of Rayleigh waves.

The traditional approaches have only focused on the peak periods. Recently, the amplitude of the H/V spectral ratio has been discussed theoretically on the basis of energy partitions for diffuse fields (Sánchez-Sesma et al., 2011). They suggested that the ratio of horizontal and vertical components of the energy is associated with the H/V spectral ratio of microtremor. Because the amplitude of the horizontal and vertical components are calculated from the material parameters of the subsurface ground, the observed amplitude of the H/V spectral ratio may be applicable to estimation of the subsurface structure (e.g. Sánchez-Sesma et al., 2011). However, the approach should address the issue that the amplitude of the H/V spectral ratio varies depending on the type of smoothing window. We, therefore, adopted only the peak periods obtained from the H/V spectral ratio of single-site microtremor observations to constrain the shallow subsurface structure.

Single-site observations of microtremors were performed over seven days, from 30 September to 2 October, and 14 December to 17 December, in 2014. The target area was 1.5 × 1.5 km² in the downtown Furukawa district, as shown in Figure 2. It almost covers the area damaged during the 2011 Tohoku earthquake. Twenty-four seismometers belonging to FuSeN (Goto et al., 2012) are included in the area. We acquired microtremor observations at 529 sites, at a spatial interval of ~50 m, which is almost the same density as Minato et al. (2014) and Goto et al. (2014).

The observation instruments that we used were high sensitivity servo-type data loggers (JU210, Hakusan Corporation) with accelerometer sensors (JA-40GA04, Japan Aviation Electronics Industry Ltd). The equipment was covered by a small box to prevent direct vibration by the wind. Microtremors were observed for 11 min at each site by default. The sampling frequency was 100 Hz. All of the candidate sites were marked roughly on a map, and then the detailed positions were selected on site.

The dataset contains two horizontal records and one vertical component record at each of 529 sites. Because the H/V spectral ratio method assumes that the ground motion is a stationary process, the observed records were analysed using the following procedures. Records are divided into segments with a duration of 20.48 s. The segments may overlap by 10.24 s. Record segments within 30 s of the beginning and the end of recording are not included in the data processing, in order to avoid contamination by artificial noise. The trend component of each segment was removed by assuming its linearity. Ten segments with the smallest root mean square (RMS) values are selected from each dataset.

Fourier amplitude spectra for the three components of the selected segments are calculated. The spectra are smoothed by a Parzen window with a 0.4 Hz window width. The H/V spectral ratio in each segment is calculated as the ratio of horizontal to vertical spectra. The horizontal value is the square root of the sum of the power in the two horizontal components. We then adopt the peak period of the averaged curve of H/V spectral ratio.

Figure 5 shows a distribution of peak periods obtained from the microtremor H/V spectral ratio calculations. The values are plotted on the map at the actual sites of the measurement. The panel also shows typical results of the microtremor H/V spectral ratio, in which black lines are spectral ratios for each segment, and the red line is their average. The spectral ratios for all 10 segments show almost the same shape and peak periods. Peak periods in the centre of the target area are relatively short, ~0.5 s, while in the northern and southern areas the peak periods are relatively long, close to 1.0 s. The southern area where long peak periods were observed is well correlated with the area of structural damage during the Tohoku earthquake, whereas damage was not reported in the northern area, where almost the same peak periods are observed. This implies that the distribution of peak period only is insufficient to directly discuss the relationship between subsurface structure and earthquake damage.

Fig. 5.  Distribution of peak periods obtained from microtremor H/V spectral ratio analysis. Four samples of H/V spectral ratios are also shown (red line: averaged curve, black line: spectral ratio for each individual segment).

The significance of the peak period differences should be quantitatively discussed on the basis of a statistical model. Spatial statistics modelling is one of the efficient tools to treat spatially distributed samples on the basis of a hierarchical Bayesian model (Banerjee et al., 2014). A conditional autoregressive (CAR) model is applied to consider spatially correlated random effects (e.g. De Oliveira, 2012).

Let T_i be the peak period at site i. The hierarchical Bayesian model for T_i is constructed as follows:

where N(,) is the normal distribution. The expectation μ_i is represented by a sum of the average μ and the spatial fluctuation Δμ_i, which is modelled by CAR as described in Equation 3. J_i is a set of neighbour samples around the site i, excepting site i itself. τ is a parameter to represent the spatial roughness.

The unknown variables, σ, μ, τ, are estimated from the observed data of the peak periods. Because σ and τ are trade-off variables, σ² is explicitly set to be 0.05 a priori, estimated from the sample variance. Prior distributions of μ and τ² are selected, as follows:

where IG(,) is inverse Gamma distribution. We numerically generate samples that follow an a posteriori distribution of μ, Δμ_i and τ from the likelihood based on the observed peak periods and prior distributions. The Metropolis-Hastings algorithm (Hastings, 1970), which is one of the efficient techniques for the Markov chain Monte Carlo method, is applied to generate the samples. We select neighbours within a circle of 100 m radius.

Figure 6 shows distributions of the 10% and 90% percentiles of exp(μ_i) for the 100 m case, which are used to discuss the significance of the spatial differences. In one example, the centre of the map of the 90% percentile shows peak periods of ~0.6 s, while the northern area of the 10% percentile map shows ~0.8 s. This indicates that the peak periods between the two areas are significantly different. We, however, cannot conclude that significant differences exist between the damaged area and the other areas in the spatial statistic model. This is because the spatial difference between peak periods is minor in the target area, and it may be buried in its variation even if it exists. Therefore, we can conclude that site classification based on the H/V spectral ratio only is difficult in the Furukawa area.

Fig. 6.  Distribution of 10% (left) and 90% (right) percentiles of exp(μ_i) for the 100 m case.

Inatani et al. (2013) estimated shallow subsurface structure on the basis of ground motion records before November 2012, observed in by the FuSeN network. They assumed a three-layered structure based on the velocity logging data at the MYG006 station, and estimated S-wave profiles at each station by comparing synthetic to observed spectral ratios. The reference site was the F17 station which is the closest site (~65 m) to MYG006 (see Figure 2). Because the velocity log was limited in 20 m depth, the basement (with S-wave velocity of 400 m/s) was set at 17 m depth at the MYG006 station, and differences between the shallower structures was assumed to contribute to the spectral ratio. However, a deeper basement with S-wave velocity of 600 m/s was estimated from microtremor array observations, and it varies spatially among three sites. Therefore, we revise the shallow subsurface structure in Furukawa district by taking into account the deeper basement.

We assume a four-layered model, which consists of three subsurface layers with S-wave velocities of 70–140 m/s (layer #1), 250 m/s (layer #2), and 400 m/s (layer #3), and a half-space basement with S-wave velocity of 600 m/s . The thickness of each layer (#1–3), and the S-wave velocity of layer #1 can vary spatially. This model is consistent with the velocity model estimated from microtremor array observations (see Table 1). Velocity logging data at MYG006 station indicates 130 m/s S-wave velocity and 17 m thickness of layer #1 and no intermediate layer #2. The thickness of layer #3 is determined via forward modelling by comparing the peak periods of microtremor H/V spectral ratios with peak periods from calculations of the theoretical ellipticity of the Rayleigh-wave fundamental mode. Figure 7 shows the comparison with the microtremor H/V ratio observed at the site closest to the MYG006 and F17 stations. The peak period of the H/V ratio (0.76 s) is closest to the peak period of ellipticity from the model (0.61 s) shown in Table 2. We here adopt this as the velocity model at the F17 station. More precise modelling at the F17 station may be required when conducting microtremor array observations.

Fig. 7.  Comparison of microtremor H/V spectral ratio and theoretical Rayleigh-wave ellipticity. The microtremor is observed at the closest site to MYG006 and F17 stations. The amplitude ratio of Rayleigh-wave is calculated from assumed model at F17 station.

Table 2.  Assumed velocity structure at MYG006 and F17 stations.

The velocity models under each station belonging to FuSeN are estimated from the observed spectral ratio with F17 records. We selected 65 event datasets for which the averaged PGA exceeded 10 cm/s² before August 2015. The target spectral ratio is calculated from the average over the spectrum of ratios to the F17 horizontal component records. A genetic algorithm is applied to search for the optimal S-wave velocity in layer #1, and the layer thicknesses of layers #1–#3, by minimising the residual between observed and synthetic spectral ratios in a frequency range of 0.5–5.0 Hz. The synthetic ratios are calculated using the Haskell-Thomson matrix method (Haskell, 1960). Density models for each station are assumed to be the same as at the F17 station (see Table 2). Damping ratios in layers #1–3 are set to be 2%. Figure 8 shows a comparison between the observed and synthetic spectral ratios at each station. Figure 9 and Table 3 show the estimated velocity structures. The new models are best constrained in the north-eastern area (F21, F26, F32, etc.) because the velocity profile estimated at array site SXZ is used to set the S-wave velocities of layers #2–#4. The models satisfy the observed spectral ratios well, especially between 2.0 and 4.0 Hz.

Fig. 8.  Comparison of observed and synthetic spectral ratios relative to the F17 station.

Fig. 9.  Comparison of the original models by Inatani et al. (2013) and the updated velocity models.

Table 3.  Estimated velocity model at F01–F36 stations of the FuSeN network.

We have thus obtained four-layered models at the 36 FuSeN station sites and at the three microtremor array observation sites. The models are now spatially interpolated using the median values of the peak periods of the single-site microtremor observations, using a statistical model. We create Delaunay triangles (Delaunay, 1934; Lee and Schachter, 1980) whose nodes are located at array observation sites and FuSeN stations, and find a triangle that includes each single-site observation site. The S-wave velocity and thickness of layer #1, and the thicknesses of layers #2 and #3 are linearly interpolated from the values at the triangle nodes. Sites not belonging to any triangles are omitted in this study, and thus 452 sites are adopted from this analysis.

However, the models are not tested to match the lowest peak period (0.4–1.0 s) of Rayleigh-wave ellipticity to the peak periods of observed H/V spectral ratios. Where observed H/V spectral ratios are available, model space is searched in the range of ± 25% and ± 50% for the optimal thicknesses of layer #2 and layer #3, respectively, and then the model between the H/V sites is linearly interpolated. RMS differences of the peak period between observed and calculated values improve to 0.092 s from 0.158 s before the correction of thicknesses of layers #2–3. Table 4 is the summary of datasets, target model parameters and frequencies, and the number of available data in modelling the shallow subsurface structure.

Table 4.  Summary of datasets for the shallow subsurface structure in Furukawa district.

Figure 10 shows the spatial distribution of S-wave velocity and depth to interfaces. The S-wave velocity of layer #1 in the southern area is ~70 m/s, while the layer thickness is less than 10 m. On the other hand, the second and third interfaces are deeper than 30 m and 70 m, respectively. This trend is also seen in the models at FuSeN stations (circle symbols), which are not updated in the procedure based on the peak periods of the H/V spectral ratio. This implies that shallow structure is spatially different in the southern area, which corresponds to the damaged area indicated in Figure 1.

Fig. 10.  Velocity models in Furukawa district. Circle symbols show the velocity models at each FuSeN station (Table 3). Square symbols show the models obtained from microtremor array observations (Table 1).

Figure 11 shows the averaged value of transfer functions calculated from the each model derived in the previous procedures. The frequency range of the averages is set to be 0.5–1.0 Hz, 1.0–2.0 Hz, and 2.0–4.0 Hz. The larger values of the average over 2.0–4.0 Hz are distributed in the southern and eastern parts. In the south, severe structural damage was not reported. The tendency was also seen in the values calculated from the previous models by Inatani et al. (2013), whereas the absolute values may be reliable under the revision. The frequency range of 2–4 Hz is higher than the effective frequency range for residential houses of ~0.67–1.0 Hz (e.g. Sakai 2009), but the frequency shift may be explained by material non-linearity due to strong motion excitation. Evidence of soil non-linearity was shown in Goto and Morikawa (2012). The peak periods of the H/V spectral ratio between JMA Furukawa and K-NET MYG006 are less than 0.5 s (2.0 Hz) for small events, but were shifted to 0.5–1.0 s during the main shock. The nonlinear site response of shallow subsurface structures is a topic for our future work. The effect on the distribution of damage in the Furukawa district will be discussed on the basis of the nonlinear dynamic response of soil and structures.

Fig. 11.  Average values of transfer functions in 0.5–1.0 Hz, 1.0–2.0 Hz, and 2.0–4.0 Hz calculated from the each velocity model. The top panel shows the average values corresponding to the previous models by Inatani et al. (2013).

We estimated shallow subsurface structure in Furukawa district, Osaki, Japan, where severe residential damage was reported during the Great Eastern Japan Earthquake of 2011. Single-site and array observations of microtremor were conducted in the area. The array observations estimated the phase velocity of Rayleigh waves at three sites, and S-wave velocity models are established. The velocity models beneath the FuSeN stations were constrained to fit the H/V spectral ratio of FuSeN records. The single-site observations estimated the peak period distribution in the surface layer on the basis of the H/V spectral ratios. We then revised the shallow (~100 m) subsurface structure on the basis of the microtremor results and ground motion records of the FuSeN stations. The overall model implies that slower S-wave velocity in the surface layer, and deeper basement are to be found around the southern and eastern areas. Damage to residential structures was concentrated in the area where the average value for the transfer functions in the frequency range of 2 to 4 Hz is large.