Simulation of soil liquefaction distribution in downtown Mashiki during 2016 Kumamoto earthquake using nonlinear site response

Several sites located between Road No.28 and Akitsu River in downtown Mashiki were liquefied during the mainshock of the 2016 Kumamoto earthquake. According to the building damage survey results, only a few buildings were damaged in areas proximate to the Akitsu River, where liquefaction occurred, however, serious building damage occurred in neighboring regions. Therefore, the effect of soil liquefaction on strong ground motions in Mashiki should be ascertained. Moreover, the distribution of visible and invisible liquefaction is required to be estimated as well. In this study, the distribution of depth of groundwater level in Mashiki was studied, which decreased from 14 to 0 m from northeast to southwest. Thereafter, the nonlinearities of the shallow layers at four borehole drilling sites were identified from the experimental data using the Ramberg–Osgood relationship. Subsequently, the dynamic nonlinear effective stress analysis of the one-dimensional soil column was performed to 592 sites in Mashiki between the seismological bedrock and ground surface to estimate the distribution of strong ground motions during the mainshock. First, the ground motions estimated by the nonlinear analysis corresponded to the ground motions observed at the Kik-net KMMH16. Second, the soil nonlinearity of shallow layers was considerably strong in the entire target area especially in the southern Mashiki, and the PGV distribution was similar to the building damage distribution after the mainshock. Furthermore, the estimated distribution of the soil liquefaction site was similar to the observed results, whereas certain invisible-liquefaction sites were estimated in the north and middle of the target area.

Total soil dilatancy C s Swelling index C c Compression index e 0 Initial void ratio of the soil material R e Shear strain radius X l Lower limit of liquefaction resistance Ni Number of loading iteration

Introduction
Two major events occurred in the 2016 Kumamoto earthquake sequence, whose epicenters were located in the vicinity of Kumamoto Prefecture of Japan (e.g., Asano and Iwata 2016;Yoshida et al. 2017). At 12:26 on April 14, 2016 (UTC), the foreshock with a moment magnitude (M w ) of 6.2 occurred at a focal depth of 11 km. Subsequently, at 16:25 on April 15, 2016 (UTC), the mainshock with an M w 7.0 occurred at a focal depth of 12 km. These events generated severe strong ground motions in the near-source regions (e.g., Doi et al. 2019;Irikura et al. 2020;Pitarka et al. 2020 Wakamatsu et al. 2017). Majority of the soil liquefactions were occurred in the near-river areas of Mashiki. After the mainshock, two microtremor observations were conducted to analyze the site liquefactions and building damages in Mashiki (Sun et al. 2020). Moreover, some researchers drilled boreholes in Mashiki and performed the triaxial stress experiments of shallow layer soils. After the mainshock, Arai (2017) conducted borehole drilling experiments and triaxial tests on shallow layers at the four borehole sites (sites A, M, O, and K, as illustrated in Fig. 1) in Mashiki. Additionally, the relationship between the shear strain (γ) and shear modulus reduction ratio (G/G max ) and that between the shear strain and damping ratio (h) of several shallow subsurface layers were obtained. The soil nonlinearities of these four sites would be applied to the dynamic nonlinear analysis in this study.
Additionally, the heavily damaged areas of buildings and the traces of soil liquefaction in Mashiki are illustrated in Fig. 1, where in the grid colors represent the distinct average damage ratio of the buildings; fuchsia and black stars denote the three strongmotion stations and the four borehole sites investigated by Arai (2017); red rectangles  NILIM and BRI (2016). Red rectangles and blue triangles denote microtremor sites observed in 2016 and 2018, respectively. Aqua and magenta circles denote liquefaction sites detected by field survey and image analysis in MLIT report (Ministry of Land, Infrastructure, Transport and Tourism 2017); purple and dark green circles denote liquefaction sites during the foreshock and mainshock determined by Wakamatsu et al. (2017). Fuchsia and black stars denote strong-motion stations and borehole drilling sites, and bold dotted line represents the region of heavily damaged buildings 1 3 and blue triangles denote the microtremor observation sites in 2016 and 2018 (Sun et al. 2020); aqua and magenta circles denote the liquefaction sites determined by the MLIT (2016) from the on-site investigation and aerial photo analysis after the mainshock. In addition, purple and dark green circles respectively denote the liquefaction sites resulting from the foreshock and mainshock identified by Wakamatsu et al. (2017). The black dotted line represents the severely damaged area by the 2016 Kumamoto earthquake sequence. The inset map at the bottom right of Fig. 1 illustrates the location of Mashiki, the foreshock and the mainshock, and the centroid-moment tensor information was obtained from the Japan Meteorological Agency (JMA 2021). A majority of the buildings between Road No. 28 and the Akitsu River suffered high damage ratios, whereas the buildings in the vicinity of Akitsu River experienced safer outcomes as compared that in the surrounding area. Therefore, the impact of the amplification and liquefaction of shallow layers in the near-river area on the strong ground motions during the mainshock needs to be determined, and damage to the wooden buildings would be estimated considering the specific building construction period. Seed and Idriss (1971) proposed a method to asses liquefaction resistance of soils. The factor of safety against liquefaction is determined by the ratio between the cyclic resistance ratio (CRR) and the earthquake-induced cyclic stress ratio (CSR). Mohamad and Dobry (1986) studied the undrained behavior of saturated sand. Total stress methods of dynamic response analysis were developed for determining earthquake induced shear stress histories in soil deposits and for deducing from these histories, with the aid of laboratory data, seismically induced pore-water pressures in saturated sands (Finn 1981). Boulanger and Idriss (2014) proposed a liquefaction analysis framework of liquefaction triggering procedures for cohesionless soils based on the previous studies (e.g. Idriss 1999;Idriss andBoulanger 2004, 2006). Idriss 1999 proposed an updated procedure for evaluation liquefaction potential referring to Seed and Idriss (1971). Next, Idriss andBoulanger (2004, 2006) displayed a semi-empirical procedure for evaluating the liquefaction potential of saturated cohesionless soils during earthquakes. As for the case studies of soil liquefaction, Zhou et al. (2020) studied the soil liquefaction of gravelly soils during the 2018 Wenchuan earthquake; Pokhrel et al. (2021) reported the liquefaction potential for the Kathmandu valley, Nepal; Quigley et al. (2013) reported the liquefaction in Christchurch, New Zealand. Zhu et al. (2015) and Bozzoni et al. (2021) also attempted the development of a logistic model to predict the probability of liquefaction based on geospatial variables and earthquake-specific parameters.
Since the 1964 Niigata earthquake, liquefaction and its associated problems have been comprehensively studied in Japan (Iwasaki 1986;Fukutake 2006). Additionally, nonlinear effective stress analysis (ESA) method has been proposed based on the solutions of coupled or uncoupled equations for granular solid and pore fluid, which were used to simulate soil liquefaction (Fukutake, 2011;Fukutake et al. 1990). The HiPER (Fukutake 1997;Fukutake and Kiriyama 2020) program developed by the Shimizu Corporation features the ESA method for soil deposits, whose effectiveness was proven in the soil liquefaction simulation during the 2011 off the Pacific coast of Tohoku earthquake (Fukutake and Jang 2012;. In these investigations, the acceleration waveform of the mainshock of the 2011 off the Pacific coast of Tohoku earthquake at the K-NET Urayasu station was precisely simulated. Moreover, the large observed settlement on ground surface were explained by their liquefaction analysis. Later, the validity of their approach was reinforced by the collaborative research of liquefaction experiments and analysis projects (LEAP), where the results of nine centrifuge experiments of a gentle sloped ground were successfully simulated by the ESA of SoilPlus during the LEAP-2017 centrifuge test (Fukutake and Kiriyama 2020).
In this study, the ESA was performed to simulate the strong ground motions of Mashiki during the mainshock, based on the identified shallow subsurface structures, the estimated seismic motions at the seismological bedrock (V s = 3265 m/s), and four sets of experimental soil nonlinearity in Mashiki. Then, distributions of the peak ground acceleration (PGA), the peak ground velocity (PGV), and the excess pore-water pressure ratio (EPWPR) would be presented. In addition, the estimated strong ground motions would be applied to simulate the damage probability (DP) of wooden houses in Mashiki by using the Yoshida DP estimation model (Yoshida et al. 2004;Nagato and Kawase 2004), and distribution of the averaged DP of each small grid of Mashiki would be illustrated. The authors showed a scientific procedure to study the site effects and DP estimation of regional buildings using the dynamic nonlinear analysis.

Research methods of three-dimensional nonlinear analysis
In general, the ESA (Fukutake 1997) comprises two models: hyperbolic and bowl models, wherein the nonlinear properties of shallow subsurface layers are required prior to application. A hyperbolic model extending in three dimensional is used for the stress-train relationship, and the strain-dilatancy relationship is modeled with a bowl model (Fukutake and Kiriyama 2020). Parameters of the hyperbolic stress-strain model are determined from the dynamic deformation tests (G/G max ~ γ, h ~ γ relationships), and parameters of the bowl model are determined from the liquefaction resistance tests.

Hyperbolic model and its parameters
One-dimensional (1D) site response analysis methods are widely used to quantify the effects of soil deposits on propagated ground motions. These methods can be classified into two major categories: frequency and time-domain analyses, among which the frequencydomain methods are widely used to estimate the site effects owing to their simplicity, flexibility, and low computational requirements (Phillips and Hashash 2009).
The Ramberg-Osgood model (R -O model) modified by Tatuoka et al. (1978) was applied to analyze the site response. Fukutake et al. (1990) presented the equations of the stress-train curve based on the Masing rule in a 1D site response, as presented in Eqs. (1) and (2), respectively: where xy and xy denote the shear stress and strain for the current loading cycle, respectively. 0 and 0 represent the current reversal points on the stress-strain curve, and the parameters and are defined in the R-O model as expressed in Eq. (3).
(1) xy = 1 3 where h max denotes the maximum damping constant of soil with a large shear strain. r (as well as 0.5 ) denotes the reference shear strain on the modulus-strain curve for G/ G max = 0.5.
In three-dimensional (3D) space as illustrated in Fig. 2 (Fukutake 1997), the shear stress versus shear strain relationships for the shear component and the axial difference component can be defined using Eqs. (4) and (5), respectively (Fukutake and Kiriyama 2020).
where G max denotes the initial shear modulus, and r denotes the reference strain obtained from the shear strength f using Eq. (6). Moreover, x , y , and z denote the axial shear strain in x-, y-, and z-axes, respectively. Therefore, the damping-shear strain relationship can be obtained using Eq. (7).
where h denotes the hysteretic damping parameter, h max represents the maximum damping ratio, and G denotes the shear modulus.
In addition, the G max , h max , and γ r are the three parameters that required for constructing the hyperbolic model in 3D site response analysis. Among these parameters, G max and h max are functions of effective stress. If G max and h max at a certain reference effective stress ( ′ mi ) assume the values of G maxi and h maxi , the G max and h max can be obtained using Eq. (8).

Fig. 2 Normal and shear stresses in 3D space
As the effective stress varies at each step of the incremental calculation, these parameters are calculated over time using Eq. (8). Moreover, the shear stress versus shear strain relationship expressed in Eqs. (4) and (5) varies with the effective stress. The values of G maxi , h max , and ri are listed in Table 1 and Fig. 3. The formulation of G maxi and ri was used when the mean effective stress ′ m was 1.0 kN/m 2 . These values can be determined from the stiffness reduction curve ( G G max ∼ relationship) or the damping increasing curve ( h ∼ relationship) obtained from the dynamic deformation tests.

Bowl model and its parameters
The definitions of the resultant shear strain Γ and cumulative shear strain G * are required to the model of soil deformation in 3D (Fukutake 1997;Fukutake and Kiriyama 2020). In particular, the bowl model proposed by Matsuoka (1989, 1993) focuses on the two parameters expressed in Eqs. (9) and (10), respectively.
where xy , yz , and zx denote the simple shear deformation, and x − y , y − z , and z − x denote the 3D axial deformation differential.
The soil dilatancy S v results from a certain soil particle repeatedly descending into the valley between the surrounding soil particles (negative dilatancy) and ascending on the surrounding soil particles (positive dilatancy). The total soil dilatancy can be evaluated by Eq. (11). where A, C, and D are parameters of the bowl model. G denotes the monotonic negative dilatancy (compressive strain) represented as a hyperbolic function with respect to G * and is irreversible. Γ denotes the cyclic positive dilatancy (also called the swelling strain) represented as an exponential function with respect to Г and is reversible. In particular, the G component represents the master curve determining the basic dilatancy during cyclic shearing, whereas the Γ component represents its oscillating component. 1/D is the asymptotic line of the hyperbolic curve, corresponding to a relative density of 100%. The mechanism of the bowl model is the movement of soil particles in a seven-dimensional strain space with xy , yz , zx , x − y , y − z , z − x , and S as axes. The unidirectional cyclic shearing model is illustrated in Fig. 4. The consolidation term is considered in the stress-strain relationship, and the effective stress is modeled under the undrained condition (constant volume). The volumetric strain S represented by Eq. (11) is the dilatancy component caused by shearing. Moreover, an additional volumetric strain resulting from the variation in the effective stress ′ m must also be considered. The total volumetric strain increment of soil d is expressed in Eq. (12). where d s denote the shear component, and d c denote the consolidation component.
Equation (13) can be used to calculate the consolidation component d c , assuming a 1D consolidation condition.
where C s is the swelling index, C c denotes the compression index, and e 0 represents the initial void ratio of the soil material. Under the undrained condition, Eq. (14) can be derived from Eq. (13).
If the mean effective stress in the initial shearing stage is ′ m0 , and if Eq. (14) is integrated under the condition � m = � m0 , then Eq. (15) can be derived as.
To suppress the occurrence of dilatancy under the small shear amplitude, a spherical region with a shear strain radius Γ = R e was considered in the strain space, wherein, no d G exited in this region. Furthermore, the R e was calculated by Eq. (17).
where ′ m0 denotes the mean effective stress in the initial shear. The positive excess porewater pressure does not increase when the amplitude of the stress ratio was equal to or less than X l (lower limit of liquefaction resistance). The physical interpretations of X l and R e are presented in Fig. 5.
The liquefaction resistance curve presents the relationship between number of cycles of repeated loading and liquefaction resistance of the lower limit value X l , as shown in Fig. 5a, wherein, X l represents the liquefaction resistance after several cycles. For simplicity, Fukutake (1997) assumed that excess pore-water pressure does not arise ( P w = 0 ) for repeated stress ratios less than the X l . Moreover, the six important parameters of the bowl model are listed in Table 2, which were determined by fitting to the liquefaction resistance curves, and their interpretation is illustrated in Fig. 6.
Physical interpretation of lower limit of liquefaction resistance X l and shear strain R e (Fukutake and Kiriyama 2020) C s is the swelling index; e 0 is the initial void ratio C c 1+e 0 C c is the compression index; e 0 is the initial void ratio X l The lower limit value of the liquefaction resistance. In the relationship between stress ratio ′ m0 and the number of cycles N c , it is represented by ′ m0 when N c is sufficiently large. It appears power-water pressure when ′ m0 > X l Fig. 6 Liquefaction resistance curves and the parameters of bowl model (Fukutake and Kiriyama 2020) 1 3

Workflow of this study
The HiPER program was included in the "SoilPlus" business software (engineeringeye, 2021), which was used to conduct the ESA of the Mashiki subsurface layers in this study. In particular, it provides a curve fitting program to determine the parameters to be used in the R-O model. In the bowl model, the SoilPlus provides a program to study the important parameters for application in the bowl model during performing the nonlinear analysis. Additionally, the Rayleigh damping of every soil column in nonlinear analysis is required. The workflow of the nonlinear analysis of the Mashiki sites is illustrated in Fig. 7.

Velocity structures in Mashiki
In total, 57 velocity structures in Mashiki were identified by the pseudo earthquake horizontal-to-vertical spectral ratio (pEHVR) at the microtremor observation sites and three strong-motion stations (Sun et al. 2020). Thereafter, the velocity structures of 535 gridded sites were obtained by linear interpolation of the 57 velocity structures, and the grid size was approximately 47 m × 47 m. The velocity structures identified from the ground surface to the seismological bedrock (V s > 3000 m/s) at six sites in Mashiki are presented in Table 3.

Seismic motions at the seismological and engineering bedrock in Mashiki
The diffused field theory (DFT) for earthquake was proposed by Kawase et al. (2011). Initially, Nagashima and Kawase (2018) estimated the seismological bedrock motions from the strong ground motions observed during the mainshock at the KiK-net KMMH16 in Mashiki (Fig. 1) using the DFT for earthquake. In their estimation the causality between the spectral amplitude and the spectral phase was violated, so unexpected vibration appeared at the beginning and ending of the estimated waveforms, which may cause irregular site response results for the ESA. Therefore, we used a cosine taper filter (Bloomfield 2000) for several seconds before the P-wave arrived and the last a few seconds of the original estimated motions. The comparisons of the three components of the seismological bedrock waves for the original and taper filtered motions are illustrated in Fig. 8, wherein the figures (a) and (b), (c) and (d), (e) and (f) represent the time and frequency domains in the EW, NS, and UD directions, respectively. The results revealed that the two motions ranging from 0.2 to 20 Hz were remarkably similar, which confirmed that tapering in this study does not impact the spectral amplitudes. These three tapered motions were applied to the nonlinear site responses described below.
In the velocity structure inversion conducted by Sun et al. (2020), the layer thickness between the seismological bedrock (V s = 3265 m/s) and the engineering bedrock (V s = 828 m/s) was assumed to be the same among all observation sites. Therefore, to reduce the calculation cost in the ESA, firstly, the ground motions at the engineering bedrock from the tapered motions were estimated by the linear analysis (LA) method (Yoshida 2001(Yoshida , 2014. Furthermore, the estimated seismic motions at the engineering bedrock were applied to the ESA on the subsurface layers above the engineering bedrock. The estimated outcrop motions at the engineering bedrock is presented in Fig. 9. 1 3

Depth of groundwater in Mashiki
Following the mainshock of the 2016 Kumamoto earthquake, Arai (2017) and Nakazawa et al. (2018) conducted borehole drilling in Mashiki to study the soil materials of the subsurface layers and published the depths of groundwater level of the borehole drilling sites, as listed in Table 4. Additionally, the site names were referred to previous researchs in Mashiki (Arai 2017). Akiba et al. (2019) studied the elevation of the groundwater level at Mashiki, as shown in Fig. 10. The elevations of ground surface in Mashiki were also obtained from the Geospatial Information Authority of Japan (GSI, Data and Resources), as illustrated in Fig. 11. The ground elevation map revealed that the ground surface was high in the NE region and low in the SW. Moreover, the ground elevation was typically low in the vicinity of river sites. The depth of groundwater level was obtained by subtracting the elevation of groundwater level from the elevation of ground surface, the results are exhibited in Fig. 12. Additionally, because the calculated water-table depth in the northeastern areas were deeper than that of KMMH16, where is the highest ground elevation area of Mashiki. Associating with the cross-section of subsurface structure from north to south of Mashiki (Akiba et al. 2019), depth of groundwater reduced from north to south, we thought the calculated depth of groundwater were not accurate in north, and some revisions should be conducted to them. Thus, the depths of groundwater level greater than 14 m (at KMMH16) were corrected to be 14 m, and depths of groundwater level less than 1 m were corrected to 1 m. The water-table depths of the five borehole drilling sites were obtained from the   Table 4. Furthermore, the water-table depth was observed as approximately 1.5 m at the near-river sites. In the northern region, water table was observed as relatively deep.

Nonlinear analysis of four borehole drilling sites
Nonlinear analyses (NA) of the soil columns in Mashiki were performed based on the obtained data using the ESA. With the SoilPlus software, the nonlinear soil properties of four borehole drilling sites were analyzed according to the experimental results reported by Arai (2017 (Arai 2017). The dynamic triaxial compression tests were conducted to determine the dynamic deformation characteristics of the undisturbed soil samples collected in the boring survey at each site. A specimen (a dense cylinder having a diameter of 50 mm and a height of 100 mm) was made by the trimming method, and then set in a triaxial cell to be saturated and isotropically consolidated. For the consolidation stress, the assumed value of the effective loading pressure at the sampling position was used. After the consolidation was completed, the P-wave velocity and S-wave velocity were measured by the bender element method (Ferreira et al. 2007), and cyclic loading was performed. The stress and strain obtained in the dynamic triaxial compression test are "axial stress" and "axial strain", and the deformation coefficient (Young's modulus) is calculated as a constant of both. Therefore, they are converted into "shear stress", "shear strain", and "shear rigidity ratio" via Poisson's ratio. Moreover, the Poisson's ratio was assumed to be 0.5 for soils at shallower than the groundwater levels and 0.3 for soils at deeper than the groundwater levels.

Nonlinear analysis of soil column at KMMH16
The NA was conducted to the subsurface layer from the ground surface to the engineering bedrock. First, with tools of the SoilPlus software, the eigenvalue analysis of the soil column at KMMH16 was used to calculate the Rayleigh damping of the shallow subsurface layers above the engineering bedrock. The results are presented in Table 5. The Rayleigh damping (Bathe 2006) of KMMH16 is illustrated in Fig. 13, wherein R and R denote the coefficient factors of stiffness and mass, respectively, and the coefficient factors of KMMH16 are listed in Table 6.
Moreover, the relationships of G/G max − γ and h − γ of shallow subsurface layers at site K, which were in the proximity of Kik-net KMMH16, were obtained from triaxial tests. Based on the SoilPlus software, the parameters of the R-O model were calculated through the curve-fitting for all types of shallow layers. The curve-fitting results based on the experimental data are illustrated in Fig. 14. Soils 1-3 represent the experimental data of borehole drilling soil at site K, whereas the parameters of the gravel layer represent the experimental results obtained from the Easy-to-understand Soil Mechanics Principles Revision Editorial Committee and Geotechnical Society (1992), which was included in the SoilPlus program.
Subsequently, the bowl model parameters of KMMH16 were obtained referring to the results of Fukutake (1997) and Fukutake and Kiriyama (2020). As compared to the borehole drilling results of site K, the third layer (Soil 3 in Fig. 14) comprised the sand, which exhibits the possibility of liquefaction during the mainshock. The relationships between the cyclic stress ratio (CSR) and the number of loading iterations (Ni), between the shear stress and shear strain, between the shear stress and effective stress of the bowl model of the third layer are illustrated in Fig. 15. Thereafter, all the parameters of subsurface layers from the ground surface to the engineering bedrock were obtained in accordance with studies of Fukutake (1997) and Fukutake and Kiriyama (2020), which are listed in Tables 7 and 8. Because the water-table depth of KMMH16 was 14.0 m, layers above the engineering bedrock were divided into two sections: above and below   (Sun et al. 2020(Sun et al. , 2021. Moreover, each subsurface layer was sliced into several pieces during the NA. The response results of each layer are displayed in Fig. 19, and the excess pore-water pressure ratio (EPWPR) of the layer during the ESA is depicted in Fig. 20. The EW, NS, and UD estimations shows the NA results were close to the observations for both the waveform and spectra comparisons at KMMH16. Figure 19 shows the response results of EW were stronger than that of NS; the soil nonlinearity was efficient in the shallow layers, especially in layers above 20 m. In addition, as observed, the liquefaction did not occur at KMMH16 because the with the maximum EPWPR was less than 100%.

Nonlinear analysis of soil column at KMMP58
KMMP58 is a strong-motion station which situated on the first floor of the Mashiki Townhall, wherein several researchers from Japan (Arai 2017;Nakazawa et al. 2018) have performed borehole drilling (site O on Fig. 1) near the station after the mainshock. In particular, Arai (2017) conducted the dynamic triaxial compression test of the sediment soil material at this site. Following a similar analysis procedure for KMMH16, the  Fig. 22. For KMMP58, soils 2 and 3 (layers 2-4) exhibited a high tendency for liquefaction, however, because the water-table depth at KMMP58 was 6.5 m, only layer 4 was under the water table at KMMP58, thus the parameters of the bowl model are illustrated in Figs. 23 and 24. All the parameters of each subsurface layer for the nonlinear analysis are listed in Tables 11 and 12. Following the NA of the subsurface layers above the engineering bedrock at KMMP58, the estimated ground motions in the EW, NS, and UD directions were compared with the strong ground motions observed during the mainshock at KMMP58, as illustrated in Figs. 25-27. Although the estimated accelerations of the three components corresponded to the observations in the EW and NS directions, the estimated velocities were weaker than the observations. In UD, the NA-estimations were close to the observations. In addition, the estimated motions were smaller than the observations in the 0.1-1.5 Hz range, while the estimations were larger than the observations for frequencies greater than 1.5 Hz. The results of the nonlinear site response for the soil column at KMMP58 are illustrated in Fig. 28. Evidently, the EW response was stronger than the NS component, and soil nonlinearity was strong for shallow layers above 20 m deep. The EPWPRs of all the layers are less than 100%, and calculated as per the maximum EPWPR accumulation procedure during the analysis, as shown in Fig. 29. Soil liquefaction is not considered to occur at KMMP58.  1 3

Four classifications of 592 sites
The NA estimations of ground motions at KMMH16, KMMP58, MSA32, and MS9-6 sites appropriately corresponded to the observations. Based on the formation history of the subsurface layers,, terrain (National Institute for Land and Infrastructure Management and Building Research Institute 2016), the water-table depth distribution of Mashiki, and the distance from each sites to the four borehole drilling sites, the nonlinear soil properties of all the 592 grid sites in Mashiki were classified into four categories, as illustrated in Fig. 30. After categorization, the red, orange, green, and blue sites utilized the nonlinear properties of site K, O, M, and A in the NA of soil columns, respectively.

PGA, PGV, EPWPR, and maximum shear strain distributions of Mashiki
The NA was applied to the 1D velocity structures of the 592 sites, And the strong ground motion of each site was estimated with the nonlinear soil properties of each category. Thereafter, the PGA and PGV distributions in Mashiki by NA were obtained, as illustrated in Fig. 31. The large PGAs were concentrated in the two northern regions, that were distinct from the PGAs observed in the NS direction. In particular, the NA PGV distributions are similar to but less than those obtained in the ELA estimated results (Sun et al. 2020). Furthermore, the NA results suggested that the large PGVs were restricted in a region corresponding to the area in Mashiki with heavily damaged buildings (Fig. 1). These results implied the PGV as an essential parameter for evaluating site risk in Mashiki. Additionally, the maximum EPWPR of each site was obtained from the NA for each site, as illustrated in Fig. 32. The sites near the Akitsu River were prone to liquefaction at both the southern area sites and certain eastern sites. However, the soil liquefaction did not occur for a majority of the northern sites, except for a few sites. The distribution map of EPWPR was consistent with the results reported by the MLIT survey and Wakamatsu et al. (2017). Furthermore, distribution map of the maximum shear strain of Mashiki is presented in Fig. 33, which intends the significant soil nonlinearity of the southern area during the dynamic ESA because the maximum shear strains of several sites were greater than 1%.

Building damage probability estimation in Mashiki
The damage probabilities (DP) of the Mashiki houses were evaluated after estimating the ground motions of 592 sites in Mashiki during the mainshock. The damage probability analysis model of wooden houses were proposed by Nagato and Kawase (2004) and Yoshida et al. (2004), Sun et al. (2021) has presented the DP estimation using the ELA-estimated ground motions, which accounted for the building construction period of each house in Mashiki and proposed an interpolation method to derive the DP based  on the DP analysis of the Yoshida model (Yoshida et al. 2004). Subsequently, the DP of each house was obtained according to the estimated ground motions from the NA site responses. The DP distribution of buildings built before 1950 (DP 950 ), 1951-1970(DP 970 ), 1971-1980 ), and after 1981 (DP 990 ) are illustrated in Fig. 34, Which  Fig. 16. Gray, orange, blue and red lines denote the observed, LA-, ELA-, and NA-estimated results, respectively reveals that majority of the heavy damage were occurred in buildings constructed before 1970. In particular, heavy building damage was primarily concentrated in the southern region of Mashiki. During every construction period, the estimated DPs of the NA were smaller than the DPs of the ELA. Thereafter, the DP of each building was calculated, and the average DP of each small region (DP cell ) in Mashiki was obtained by referring to the region size reported by Sun et al. (2021). The distribution map of the DP cell of Mashiki is illustrated in Fig. 35, which implied that the DP cell of NA is lower than the that of ELA. Moreover, large DPs were predominately located in the southern Mashiki, which is coincident with the AIJ building damage survey report. Thus, we suggest that a large DP may appear in the northwestern region during the mainshock. Moreover, the estimated ground motions were diminished in the NA, as they are affected by soil nonlinearity. According to the analysis results of MS9-6 (supplement), the soil liquefaction reduced the amplification of strong ground motions. Thus, the damage of buildings in the liquefaction areas may decreased by the soil liquefaction because the amplifications of strong ground motions were decreased in 0.5−2.0 Hz, especially in the near-river areas. In addition, new wooden houses were built in the near-rive areas which may also decreased the damage of buildings in the near-river area. Overall, the estimated DP distribution map depicts the building damage distribution in Mashiki during the mainshock.

Discussion
Four sets of soil materials for the shallow subsurface layers were applied to the NA of 592 sites, which improved the accuracy of the analysis. The NA-and ELA-estimated ground motions at KMMH16 and KMMP58 were slightly weaker than the LA estimations, and both the NA and ELA estimates closely represented the observed ground motions, as shown in Figs. 16,17,18 and Figs. 25,26 and 27. Additionally, the soil nonlinearity reduces the amplification of strong ground motions. The NA-estimated ground motions at MS9-6 more or less corresponded with the observed strong ground motions at KMMP58 (Figs. S14−S16), indicating the parameters of soil material and ESA were suitable for the analysis. Moreover, the PGV distribution obtained by NA displayed a similar trend to that of the building damage map by AIJ. Generally, areas with large PGV exhibited heavy building damage, which implies that the PGV explains the building damage ratio of wooden houses in Mashiki. As the NA-PGV was slightly smaller than the ELA-PGV, the NA-PGV distribution displayed a more appropriate pattern than the ELA-PGV distribution when compared with the building damage distribution of the AIJ. Furthermore, these findings indicate that ELA-and NA-estimations have some difference in Mashiki, especially in the large PGV areas, and the analysis method must be appropriately checked prior to the site response analysis. Overall, the liquefaction sites were located near the Akitsu River in Mashiki according to the MLIT survey and Wakamatsu et al. (2017). Distributions of soil liquefaction sites were in accordance with the distribution of observed soil liquefaction site. Moreover, several sites in northern and central Mashiki displayed a high percentage of liquefaction during the mainshock. The distribution map of soil liquefaction revealed that the site locations with an EPWPR ≥ 99% were concentrated in southern Mashiki. Moreover, NA−DP distribution was in accordance with the damage distribution of the AIJ report. The final NA-estimated DP was slightly smaller than the observed damage ratio in southern Mashiki, which indicates the anti-seismic property of buildings may be worse than we expected in the model.
Although the 1D NA displayed that the hybrid effects of high PGV and soil liquefaction may have caused the heavy building damage in Mashiki, the 3D topographic effects was not considered in this study. However, the results may not vary widely from that of the 1D analysis because we applied an extremely dense NA of soil columns. Nevertheless, the authors will continue to perform the 3D NA of subsurface layers in Mashiki.

Conclusion
This study presents a detailed procedure to estimate the building DP from the microtremor observation. In this process, the effective stress analysis was applied to the shallow subsurface layers of Mashiki, and the water-table depth of Mashiki were determined and soil materials of shallow soil layers for four categories of soil columns in Mashiki were defined based on the experimental data. Additionally, we present the PGA and PGV distributions obtained by NA and conducted simulations of the soil liquefaction sites. The building DP was evaluated considering the building construction period, and the distribution map of Mashiki was simulated.
(1) The NA-estimated ground motions were similar to the observed strong ground motions during the mainshock in both the time and frequency domains, which verified that the analysis method is suitable for reproducing the strong ground motions in Mashiki.
(2) The NA-estimated ground motions were similar to but smaller than the ELA-estimated ground motions, which implies that the nonlinear soil properties of shallow layers were stronger than the ELA used in the high shear strain range. Moreover, the distribution of maximum shear strain displayed the high shear strain were mostly located in the southern areas. The results confirmed the necessity of applying ESA for soil columns of shear subsurface layers in Mashiki.
(3) The NA−PGV distribution maps displayed a pattern similar to that of the AIJ building damage distribution map; however, the NA −PGA distribution maps did not exhibit    , (e), (f), and (g) present estimated absolute peak acceleration, estimated relative peak velocity, estimated relative peak displacement, estimated EPWPR, G/G max and damping ratio, maximum shear strain, and estimated shear stress, respectively, for every layer such close relationships, indicating that the seismic resistance ability of wooden houses in Mashiki were not as strong as we expected, and the PGV is an essential parameter explaining the building damage in Mashiki. (4) The estimated site distribution map of soil liquefaction in Mashiki was similar to that of the observed results. Furthermore, the estimated liquefaction locations were wider than the observed ones, and several sites were located in the same region of the heavy building damage reported AIJ. These results indicated that the soil liquefaction may have effects on the building damage in the southern Mashiki. (5) The NA −DP estimations were smaller in areas with heavy building damage. However, upon collating with the soil liquefaction estimation map, the estimated NA −DP map corresponded to that of the AIJ survey results of building damage. Institute for Land and Infrastructure Management (NILIM), and Building Research Institute (2016) can be found at https:// www. kenken. go. jp/ japan ese/ conte nts/ publi catio ns/ data/ 173/ index. html (in Japanese, last accessed August 2021). "DYNEQ" software was downloaded at https:// www. kiso. co. jp/ yoshi da/ Japan ese_ 02. html (last accessed August 2021). "SoilPlus" software was downloaded at https:// www. engin eering-eye. com/ SOILP LUS/ (last accessed August 2021). Several figures were generated using the QGIS3 software and GMT5 software, which can be found at https:// www. qgis. org/ en/ site/ (last accessed August 1 3 Intellectual Property We confirm that we have given due consideration to the protection of intellectual property associated with this work and that there are no impediments to publication, including the timing of publication, with respect to intellectual property. In so doing we confirm that we have followed the regulations of our institutions concerning intellectual property.

Research Ethics
We confirm that this work does not have any conflict with social ethics. All research data of this study have been approved by relevant departments.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.