The author applied the site-effect substitution method. The horizontal components were estimated based on the horizontal spectral ratio (Fig. 3). The vertical component was estimated based on the vertical spectral ratio (Fig. 4).

In Fig. 6, the acceleration waveforms estimated and submitted by the author for the blind prediction Step 2 are compared with the recorded acceleration waveforms, which were distributed by the organizer after the blind prediction. The origin of the horizontal axis corresponds to April 16, 2016, 3:03:10 (JST). It can be pointed out that there is an error in the arrival time of the main energy. When the author submitted the results, the author did not make any rigorous investigation in terms of the origin of the time history. After the recorded waveforms were distributed, the author found that, although the events No.35 and No.104 occurred at similar locations (Fig. 1 and Table 1), the arrival of the S waves were delayed for the event No.35 compared to the event No.104 both for KUMA and JMA. The author cannot find any good explanation for this tendency. Therefore, in this article, the author decided to simply shift the estimated waveforms so that the arrival time of the S wave be consistent between the estimated and recorded waveforms. Figure 7 compares the estimated acceleration waveforms after the time shift with the recorded acceleration waveforms. Figure 8 compares the estimated velocity waveforms after the time shift with the velocity waveforms calculated from the records. Figure 9 compares the estimated Fourier spectra with the Fourier spectra calculated from the records. These comparisons reveal that, except for the error in the arrival time, the estimated ground motions were fairly consistent with the observed ground motions. Especially the velocity waveforms were reproduced well including the later phases. The result indicates the effectiveness of the site-effect substitution method in estimating ground motions after an earthquake at a site where the ground motions are unknown as long as aftershock records are available at the site.

It should be noted, however, that two condition contributed to this favorable result in Step 2. One is the fact that the rupture process of the target event was simple. As can be seen in Fig. 5, the Fourier phase spectra were interchangeable between the events No.35 and No.104. This means that the Fourier phase characteristics were almost determined by the path and site effects and not by the source effects. This is the evidence for a simple rupture process of the target event No.35. Another condition was that the soil nonlinearity at the target site was negligible. This was already discussed in terms of the bottom panel of Fig. 3 and supported by the accurate reproduction of the Fourier spectra in Fig. 9. These conditions were favorable for the application of the site-effect substitution method.

### Step 3 – Simulation Of Strong Motions

The goal of Step 3 was to estimate ground motions at KUMA during the M6.5 foreshock and the M7.3 mainshock of the Kumamoto earthquake sequence. For this step, the organizer of the blind prediction suggested us to pay attention to the effects of soil nonlinearity. Therefore, the author followed the following procedure:

(1) Estimate surface ground motions in which linear soil response is assumed.

(2) Estimate 2E ground motions at GL-31m assuming linear multiple reflections.

(3) Estimate surface ground motions with effective stress analyses.

For (1), the author estimated the Fourier amplitude and phase characteristics separately. For (3), the author used a program called “FLIP” (Iai, 1991; Iai et al., 1992). This procedure implies that the effect of soil nonlinearity at the reference site JMA, if any, is neglected. This fact will be discussed later.

## Applicability of the site-effect substitution method

There were several candidate methods to estimate ground motions during the foreshock and the mainshock including ground motion simulations based on fault models. Among them, the author first examined the applicability of the site-effect substitution method, which was used in Step 2.

One of the main features of the site-effect substitution method is that the Fourier phase spectrum at the target site for the target event is approximated by the Fourier phase spectrum at the same site for another event. This operation is effective if the two events share the same path and site effects and the source effect is not a dominant factor for the Fourier phase spectra. These conditions applied to Step 2. However, it was anticipated that the latter condition does not apply to Step 3 due to the complexity of the rupture process especially for the mainshock (Asano and Iwata, 2016; Nozu and Nagasaka, 2017). Therefore, the author investigated the similarity of the Fourier phase characteristics between the target events and other available events at JMA.

In Fig. 10, the similarity of the Fourier phase characteristics is examined between the mainshock and the event No.161 at JMA. Among available events, the event No.161 produced the best results for the mainshock. This is presumably because the event No.161 (Fig. 1) occurred near the asperities of the mainshock (Asano and Iwata, 2016; Nozu and Nagasaka, 2017) and this event conveyed similar path and site effects with the mainshock. However, the large-amplitude phase for the EW and UD components at approximately 27 s was not captured by the phase-exchanged waves. This is presumably because of the complexity of the rupture process of the mainshock. Therefore, the author avoided to use the site-effect substitution method for the mainshock.

The similarity of the Fourier phase characteristics was also examined between the foreshock and other available events at JMA (Figs. 11 and 12). Among available events, the events No.64 and No.161 produced the best results for the foreshock, presumably because these events occurred near the asperities of the foreshock (Asano and Iwata, 2016) and they conveyed similar path and site effects with the foreshock. The similarity of the waveforms indicated the potential applicability of the site-effect substitution method to the foreshock. However, the author avoided to use the site-effect substitution method for the foreshock simply because the author preferred to use the same method for the foreshock and the mainshock.

## Availability of the Fourier phase spectrum at JMA

Instead of using the site-effect substitution method, the author tested the similarity of Fourier phase characteristics between KUMA and JMA for several events.

Figure 13 shows the results for the event No.161. The Fourier phase characteristics at KUMA and JMA resemble each other at least for the horizontal components, although there is a notable difference for the vertical component. The author found that the degree of similarity depends on the event. However, the author acknowledges the importance of the event No.161 because the event presumably conveys similar path and site effects with the foreshock (Fig. 12) and the mainshock (Fig. 10). If the Fourier phase characteristics are similar for the event No.161 between KUMA and JMA, it could be reasonable to assume that they are similar also for the foreshock and the mainshock. Therefore, to estimate ground motions at KUMA during the foreshock and the mainshock under linear site response, the author simply used the Fourier phase characteristics at JMA for the same events.

## Fourier amplitude spectrum

The Fourier amplitude spectra at KUMA during the foreshock and the mainshock were estimated in a similar way as in Step 2, that is, the Fourier amplitude spectra observed at JMA during the same events were multiplied by the spectral ratio shown in the top right panel of Fig. 3 for the horizontal components and the right panel of Fig. 4 for the vertical component. The author used the averaged ratios for the five events. To make sure that it is appropriate to use the averaged ratios for the five events, which did not necessarily occur near the target events, the author compared the averaged ratios for the five events with the spectral ratios for the event No.161, which presumably conveyed similar path and site effects with the target events as mentioned above. The results indicated that there is no significant difference between the spectral ratios at least for the horizontal components. Therefore, the author used the averaged ratios for the five events.

## Incident (2E) ground motions at GL -31 m assuming linear multiple reflections

Next, the incident (2E) ground motions at GL-31m were estimated assuming linear multiple reflections. The soil layers were divided as shown in Table 2 according to the PS logging results provided by the organizer of the blind prediction (Matsushima et al., 2022). For the horizontal components, the S-wave velocities were used. For the vertical components, the P-wave velocities were used. For the 2nd, 4th, 5th, 6th and 7th layers, laboratory test results were available. The corresponding samples were T-1, T-2, Tr-3, Tr-4 and Tr-5, respectively. The author used the wet densities from the laboratory tests for these layers. For the 3rd, 8th and 9th layers, the author decided to use the same wet densities as the 6th, 6th and 5th layers, respectively, considering the similarity of the soils and N values. For the 1st layer and the “bedrock”, the author assumed a density of 2.0 g/cm3. The author used the damping coefficient of 0.02 for all the layers for all the analyses, referring to the results of cyclic triaxial tests.

## Surface ground motions with effective stress analyses

Finally effective stress analyses were conducted to estimate surface ground motions under nonlinear site response. For this purpose, the author used a program called “FLIP” (Iai, 1991; Iai et al., 1992), which has been extensively used in the design of port structures in Japan. The parameters were determined based on a simplified procedure (Morita et al., 1997) from the PS logging results, N values and fine particle content provided by the organizer (Matsushima et al., 2022).

Table 3 shows the parameters used for the effective stress analyses. In this program, the shear modulus and the bulk modulus for the soil skeleton at small strain depend on some power of the effective confining pressure. They are specified for a reference effective confining pressure \({\sigma }_{m0}^{{\prime }}\). In this analysis, the effective confining pressure at the center of each layer was selected as the reference effective confining pressure. Then, the shear modulus at small strain \({G}_{m0}\) for the reference effective confining pressure was determined from the PS logging results. The bulk modulus for the soil skeleton \({K}_{m0}\) for the reference effective confining pressure was determined assuming a Poisson’s ratio of 1/3. Then they were assumed to be proportional to the square root of the effective confining pressure, which is a standard assumption in this program. For the sand and sandy silt (the 2nd, 3rd, 5th, 6th, 8th and 9th layers), \({\varphi }_{f}\) was determined from the N values based on Morita et al. (1997). For the silt and clay (the 4th and 7th layers), a shear strength equivalent to \({\varphi }_{f}\)=30° was used. For the porosity \(n\) and the maximum damping coefficient \({h}_{max}\), the values ordinarily used for FLIP analyses were used.

For the hatched layers in Table 3, the effects of excess pore pressure were considered. The parameters that control the excess pore pressure were determined from the N values and fine particle content based on Morita et al. (1997). The fine particle content was based on the laboratory test results for the 5th and 6th layers. For the 3rd, 8th and 9th layers, the author decided to use the same fine particle content as the 6th, 6th and 5th layers, respectively, considering the similarity of the soils and N values.

The input ground motions (NS and UD components in one case and EW and UD components in another case) were applied via a viscous boundary. In terms of the UD component, the result for the former case was selected to be submitted. The difference for the two cases was small.