We are using two purely empirical approaches to predict the surface acceleration time histories and the associated Fourier spectrum for both the foreshock and the mainshock of the Kumamoto 2016 earthquake. Both methods use statistically defined correction factors that are applied to the site-specific linear transfer function to account for the non-linear soil behavior.
The linear transfer function (TFllin) can be obtained from weak ground motions recordings or numerical simulations. It can be borehole or outcrop transfer function. In this paper, we use the weak motions recorded at the reference rock site SEVO and at the target sedimentary site KUMA (The location of the stations and epicenters of seismic events are illustrated in Fig. 1).
FIGURE 1
The non-linear transfer function (TFNL) is corrected for the non-linear soil behavior as follows:
\({TF}_{NL}\left(f\right)={TF}_{lin}\left(f\right).{RSR}_{NL-L}\left(f\right)\) (Eq. 1)
\({TF}_{NL}\left(f\right)={TF}_{lin}\left(f.\sqrt{fsp}\right)\) (Eq. 2)
In horizontally stratified absorptive earth with vertically traveling plane waves, transmitted waves are translates of minimum phase functions (Sherwood and Trorey, 1965). We also choose for this transfer function a minimum phase assumption to associate a realistic phase to the modulus of the non-linear transfer function TFNL and obtain a complex transfer function. This complex transfer function is then applied to the Fourier transform of the acceleration time histories recorded at SEVO (derived from the velocity recordings), and an inverse Fourier transform provides the predicted acceleration time histories at KUMA. A time delay to consider the difference in seismic waves travel times between SEVO and KUMA will be considered.
Method 1 : Correction for soil non-linear behavior using RSRNL−L :
The approach is based on the work of Derras et al. (2020). in which, the authors have predicted a linear to non-linear correction curve, called RSRNL−L, to correct both the amplitude changes and the frequency shift associated to the non-linear soil behavior, with respect to the modulus of the linear transfer function derived from weak motion recordings at both KUMA and SEVO.
This correction function is estimated on the basis of a statistical analysis of a selection of KiK-net recordings, which were used to train an artificial neural network (ANN) in order to predict the frequency-dependent correction function RSRNL−L as a function of one ground motion intensity measure (GMIM) and site-condition proxies (SCPs). The ANN approach was inspired by investigations into the human brain structure, which consists of interconnected neurons. An ANN is made up of input, hidden and output layers, and the connections between neurons of two consecutive layers are characterized by linear combinations of inputs: the corresponding synaptic weights (W) and biases (b) are the learnable parameters of a ANN. The W values characterize the influence of input (GMIM, SCPs) on the output (RSRNL−L). The Quasi-Newton Back Propagation technique also called “BFGS” (Robitaille et al.,1996) has been used in this work for the training phase. To avoid the overfitting, the regularization method presented in (Derras et al., 2012) is used in this study. More details upon the structure and the design of ANN can be found in Derras et al. (2020).
The performances of the ANN results are measured by the standard deviation of residuals between observations and model predictions, compared to the standard deviation of the original advanced data set through the variance reduction coefficient Rc (Derras et al., 2017) as defined in Eq. 3
$${R}_{c}=\left[1-\frac{{\left[\sqrt{\frac{1}{M}\sum _{1}^{M}{\left({log}_{10}\left({RSR}_{NL-L,obs}\right)-{log}_{10}\left({RSR}_{NL-L,pred}\right)\right)}^{2}}\right]}^{2}}{{\left[\sqrt{\frac{1}{M}\sum _{1}^{M}{\left({log}_{10}\left({RSR}_{NL-L,obs}\right)-{log}_{10}\left(mean\left({RSR}_{NL-L,pred}\right)\right)\right)}^{2}}\right]}^{2}}\right].100 \left(\%\right)$$
eq. 3
Where RSRNL−L,obs represents the "observed" RSRNL−L as derived in the step 3 advanced data set. RSRNL−L,pred is the neural prediction of the RSRNL−L. This function is predicted either in the normalized frequency domain with a frequency parameter called fNL defined in (Régnier et al., 2016) for the sensitivity study, or in the absolute frequency domain for the final model. M is the size of the data set.
The final neural models thus consist of a series of three layers. The first layer represents N inputs, (one GMIM and {N-1} SCPs). The second, hidden layer has the same number N of neurons. The last layer represents the values of RSRNL−L for 49 absolute frequency bins. The selected architecture is therefore of the N-N-97 or N-N-49 type.
The model chosen in this article to predict RSRNL−L is based 6 parameters. 5 site parameters: Vs30 (mean shear Vs over the first 30 m), B30 (gradient of the velocity profile over the first 30 m), VSmin (the minimum shear wave velocity), f0 (the fundamental resonance frequency, picked on H/V ratio), and A0HV (the corresponding peak amplitude), and on a loading parameter characterizing the intensity of ground shaking at KUMA (GMIM) site, here PGA. Whereas in Derras et al., 2020, we have used two SCPs (Vs30 and f0HV) and PGV/Vs30 to represent the GMIM (03 neurons at the hidden layer).
Method 2 : Correction for soil non-linear behavior using fsp
Castro-Cruz et al. (2020) used a subset of the KiK-net database to define the fsp parameter. This subset is composed of 688 stations. Among them, 650 sites are characterized with Vs and Vp profiles, soil description, and information on the stations (location and information of recording devices. For each recording, the acceleration time histories are provided, with the event origin time, the epicenter location, the depth of the hypocenter, and the magnitude of the event determined by the Japanese Meteorological Agency (JMA). Recordings were selected according to two criteria: signals from earthquakes with a magnitude higher than 3 and with an epicentral distance lower than 500 km. These two criteria lead to 75 232 pairs of 3-components signals (surface and downhole) from 5535 events with magnitudes between 3 and 9, recorded between November 1997 to December 2017. For all of the selected seismograms the weak ground motions corresponding to assumed linear behavior are selected based on their maximal peak accelerations at downhole from 10− 4 m/s2 to 6. 10− 3 m/s2. The linear Borehole Transfer Function (\(BT{F}_{Lin}\)) is the arithmetic average of all of the weak ground motions borehole transfer functions. It is compared to each \(BT{F}_{NL}\) from all stronger ground motions recorded at the same site, allowing us to define a parameter that characterizes the observed frequency shift. The logarithmic frequency shift is the gap in logarithmic scale between BTFLin and BTFNL. In linear scale, it is a coefficient that changes the frequency scale. The algorithm to find this logarithmic shift minimizes the misfit between BTFlin and BTFNL as defined in Eq. (2). Note that the misfit is weighted by the logarithmic sampling.
$$misfit=\sum _{i}^{}|BT{F}_{lin}\left(\stackrel{-}{f}.Ls\right)-BT{F}_{NL}(\stackrel{-}{f}\left)\right|{\Delta }x$$
$${\Delta }x={\text{log}}_{10}\left(\frac{{f}_{i+1}}{{f}_{i}}\right)$$
$$\stackrel{-}{f}=0.5({f}_{i+1}+{f}_{i})$$
eq. 4
In Eq. 4, Ls is the frequency scaling factor applied to BTFlin. The misfit is defined as a discrete approximation of the area between the shifted BTFlin and BSR, considering a logarithmic scale as the length of the base (Δx). The computation is done over a frequency window going from 0.3 Hz to 30 Hz.
Finally, we define a frequency shift parameter, so called fsp, as the square of the Ls, which produces the minimum value of misfit (fsp = Ls2 when misfit is minimized).
Application To Esg6 Blind Tests
Between the blind prediction performed last year and the results of the predictions presented in this article, we modified only the linear transfer function TFlin and we calculated the TFNL for the two horizontal components separately instead of performing the average. First, the Fourier spectra were smoothed (using a Konno-Ohmachi smoothing) with a coefficient of 60 instead of 40 to minimize the reduction of the amplitude.
Definition of the TFlin
TFlin is the average of the weak motion ratio between the Fourier spectrum of the velocity recorded at KUMA and SEVO. We obtain a transfer function for the three recordings components.
The accelerometric data correspond to all aftershocks listed in Table 1. The signals were pre-processed following the recommendations of Boore and Bommer (2005). From the acceleration time histories at KUMA, we integrated the signal to obtain the velocity time histories. Before performing the integration, we apply pre-processing that includes removing the mean, the tapering, the addition of zeroes before and after the signal, and a 2-poles butterworth filtering between 0.1 and 30 Hz. The Fourier spectra were calculated and smoothed using a Konno-Ohmachi smoothing (Konno & Ohmachi, 1998) with a parameter of 60. We choose a parameter of 60 to minimize the decrease of the amplitude due to the smoothing. We calculated the Fourier spectra ratio between KUMA and SEVO for all seismic events and the average and standard deviation considering a lognormal distribution called TFlin. We compared the arithmetic with the logarithmic averages on the transfer functions. The arithmetic average provides lager amplitudes. However, the geometric average, that fulfills the hypothesis of data being lognormally distributed, is more adapted to seismological observations and for this reason was considered in this study. The averaged horizontal components to vertical component spectral ratio is also calculated at KUMA site and so-called HVLin. The results are illustrated in Fig. 2. The transfer functions from both components provide similar main frequencies peaks around 0.3, 1.3, 2.0 ,7.4, 16 and 26 Hz. The average H/V spectral ratio performed with the recordings at KUMA sites confirms that the fundamental resonance frequency of the site is at 0.3 Hz. Note that the variability of the linear transfer function between weak events is important. Depending on the location of the weak motion epicenters, one of the main assumptions of the site-to-reference spectral ratio, that the site-to-reference distance is much smaller than the epicentral distance, so that the path traveled by the seismic wave beneath the referent site and at the bedrock beneath the study site is equivalent, is not respected. The large variability of the weak events leads to an average transfer function that can have lower peak amplitudes. For instance, at the predominant frequency (1.3 Hz), the average amplitude is around 8.2 while the 16 percentiles (-1\(\sigma\)) is at 5.2 and the 84 percentiles (+ 1\(\sigma\)) is at 16.6 for the East-West component.
Table 1
Characteristics of the recorded earthquakes at SEVO and KUMA
Eq N° | Mw | Depth (km) | Lat. | Long. | SEVO | KUMA | \({\Delta }t \left(s\right)\) |
Depi (km) | PGA (m/s2) | Depi (km) | PGA (m/s2) |
48 | 3,3 | 15,5 | 32,716 | 130,805 | 21 | 0,03 | 13 | 0,06 | 0,98 |
20 | 4 | 12,3 | 32,797 | 130,813 | 18 | 0,06 | 12 | 0,38 | 0,84 |
32 | 4,8 | 11,9 | 32,787 | 130,774 | 15 | 0,11 | 8 | 0,51 | 1,00 |
71 | 4,4 | 14,6 | 32,758 | 130,778 | 16 | 0,18 | 9 | 0,74 | 0,87 |
69 | 4,8 | 8,9 | 32,962 | 131,079 | 45 | 0,05 | 42 | 0,21 | 0,61 |
09 | 4,5 | 10,3 | 32,687 | 130,776 | 21 | 0,07 | 13 | 0,25 | 1,25 |
22 | 4,4 | 10,6 | 32,678 | 130,721 | 18 | 0,12 | 11 | 0,18 | 1,11 |
02 | 4,2 | 10,2 | 32,870 | 130,873 | 24 | 0,05 | 20 | 0,12 | 0,56 |
39 | 4 | 11,0 | 32,785 | 130,832 | 20 | 0,08 | 13 | 0,23 | 0,79 |
47 | 4,6 | 11,2 | 33,000 | 131,134 | 52 | 0,04 | 49 | 0,13 | 0,50 |
83 | 4,9 | 10,8 | 32,993 | 131,122 | 50 | 0,05 | 47 | 0,14 | 0,66 |
28 | 3,9 | 16,4 | 32,831 | 130,814 | 18 | 0,08 | 13 | 0,18 | 0,49 |
43 | 6,9 | 11,4 | 32,742 | 130,809 | 19 | 1,07 | 12 | 4,15 | 1,04 |
47 | 7,3 | 12,5 | 32,755 | 130,763 | 15 | 1,96 | 7 | 5,42 | 0,71 |
FIGURE 2
Required site (SCPs) and input motion (GMIMs) parameters
The values of the site parameters provided in Table 2, were derived from the preferred soil model proposed by the organization of the ESG6 benchmark (Matsushima et al., 2022) and illustrated in Fig. 3. The fundamental resonance frequency f0 (Hz) and associated amplitude A0 were obtained from the empirical site response curves and confirmed by the Horizontal to vertical spectral ratio at KUMA site.
Table 2
KUMA Site-Condition Parameters (SCPs) used for the prediction of RSRNL−L.
Vs30 (m/s) | B30 | f0 (Hz) | Vsmin(m/s) | AOHV |
160 | 0.3 | 0.3 | 95 | 4.8 |
Table 2
The PGA of the foreshock and mainshock were predicted using the linear regression between the log10 of PGAKUMA/PGASEVO modulated by the ratio between epicentral distance at KUMA and SEVO (DepiKUMA/DepiSEVO). We calculated the PGA at SEVO of the foreshock and mainshock. For the mainshock recorded at SEVO, we choose to select the PGA without considering the very high amplitude observed at 23.6s because this peak appears to not be related to the vibrations due to the earthquake. Indeed Tsuno et al., (2017) suggest that the station went out of power and they measured the PGA on the signal before this spurious peak amplitude. The average PGA over the two horizontal components at SEVO are thus 1.07 and 1.96 m/s2 for the foreshock and mainshock respectively. Considering these values and the predicted site response using the weak motions recorded at both stations, the predicted PGA at KUMA is 2.36 m/s2 for the foreshock and 4.85 m/s2 for the main shock. However, after a first prediction of the Fourier spectra at Kuma, we found a PGA for the foreshock around 4.12 m/s2, we decided to update the calculation using this value to correct for the non-linear soil behavior.
Comparing these blinded predicted values with the actual recordings at KUMA site, we found that the predicted PGAs are very close to the recorded ones even underestimated for the mainshock. For the foreshock the recorded PGA is 4.15 m/s2 and 5.42 m/s2 for the mainshock which represent a discrepancy of less than 1% for the foreshock and 10% for the mainshock. It means that the nonlinearities have been slightly underestimated during the mainshock.
Correction from method 1
To predict the correction factor for the nonlinearity, we used the site parameters defined in Table 1 and the PGA predicted at KUMA. Such a set of parameters is not the one that performs the best in terms of variance reduction, but it is a good compromise between ease of use (knowing the velocity profile) and model performance: the corresponding standard deviation reduction (Rc) is 18%.
The use of the additional parameter fNL, that is a frequency from which we observe a de-amplification in the site response, would have increased Rc to 25%. However, in this case, considering that no other large events were recorded at KUMA and SEVO sites, it would have been necessary to infer a value for fNL from f0HV, and such a derivation is associated with a large uncertainty.
FIGURE 3
Correction from method 2
Following the procedure described in Castro-Cruz et al. (2020), we obtain the frequency shift for each TF with respect TFlin. The values of fsp for each horizontal component are summarized in the Table 3. The fsp is then plotted in function of the intensity of the ground motions here the PGA at SEVO has been chosen (see Fig. 4). Using a non-linear correlation (hyperbolic model) between the fsp and the PGA at SEVO, we predict the fsp for the two strongest events. In Fig. 4, the gray markers represent the fsp and PGA of the aftershocks while the two red markers represent the two target events. The fsp prediction is based on the aftershock results only and is illustrated by the solid blue line for the mean and the dashed blue lines for the associated uncertainties. Note that the predicted fsp is not constrained above the maximal PGA of the aftershocks. It differs depending on the component: for the North-South the predicted fsp is equal to one meaning that no non-linear soil behavior is predicted. While in the East-West component, although the uncertainties are very high, the average fsp predicted for the foreshock is 0.95 and 0.71 for the mainshock. Comparison with the measured fsp of the target events indicate that the method underpredicts the shift of the transfer function.
Table 3
Values of frequency shift parameter (fsp) for all events relative to the linear transfer function (TFlin)
Eq N° | Fsp NS | Fsp (EW) |
48 | 1 | 0.97 |
20 | 1.08 | 1.1 |
32 | 1 | 1 |
71 | 0.97 | 0.92 |
69 | 0.95 | 1 |
09 | 1.1 | 1 |
22 | 1.1 | 1 |
02 | 0.87 | 1 |
39 | 1.04 | 1 |
47 | 1.01 | 1.04 |
83 | 1.07 | 1.05 |
28 | 0.88 | 1.03 |
43 | 0.6 | 0.54 |
47 | 0.65 | 12,5 |
Table 3
FIGURE 4
Returning to time domain
Using the corrected transfer functions from non-linear soil behavior, the Fourier spectra of the two events can be predicted at Kuma site by multiplying the Fourier response spectra of the SEVO recordings by the predicted transfer functions. The next step is to predict the acceleration time histories. The arrival time and phase of the acceleration time histories need to be defined.
To predict the arrival time at Kuma site, we propose to predict the delay of P-waves arrival times between KUMA and SEVO for the foreshock and mainshock using the provided aftershocks. In Fig. 1, the delay of P-wave arrival (manually peaked) is illustrated with the grey color scale. Note that the events n° 47, 69 and 83 listed in the Table 1 are not located in the Fig. 1 because they are too far from the sites.
The delays depend on the location of the epicenters with respect to the location of KUMA and SEVO stations. The difference of epicentral distances (and consequently the delay) decreases when the epicenters are further from the KUMA/SEVO alignment. The closest event to the foreshock is the event 48 and for the mainshock the event 71. We therefore choose to apply a delay of 0.98s for the foreshock and 0.87s for the mainshock that corresponds to the peaked delay of the event 48 and 71 respectively.
Then, we consider the phase of the recordings at SEVO and a minimum phase assumption for the transfer function to consider the effect of seismic wave propagation vertically in a soil column. We followed the procedure described in (Pei and Lin, 2006) that consist in taking the even part of the cepstrum (i.e. real part of the inverse Fourier transform of the logarithm of the transfer function) and to reconstruct a causal sequence by multiply it by a function with values equal to 0 for the imaginary part and 2 for the real part. In time domain, the complex transfer function is the real part of the inverse Fourier transform of the exponent of the causal sequence.