Machine learning ground motion models for a critical site with long-term earthquake records

Ground motion models derived using data from widespread stations and earthquakes induce significant aleatory variability when applied to predict ground motions for a single site. This issue can be even worse when the sites of interest show complex site response characteristics and do not have enough similar sites in the ground motion model development dataset. This work provides a machine-learning-based ground motion modeling framework for sites with numerous earthquake recordings. The Onahama Port Array, a liquefiable site in Japan that has recorded more than 1200 earthquakes in the last 25 years, was selected as a case study. The gradient boosting method was employed to predict response spectra with periods from 0.01 to 10 s for earthquakes with Japan Meteorological Agency (JMA) magnitude ranging from 2.4 to 9.0. Different predictor combinations were tested, and two gradient boosting models were recommended: the basic model that requires earthquake magnitude and epicentral distance as inputs, and the optimal model that requires three additional inputs, focal depth, azimuth, and rainfall. The basic and optimal gradient boosting models have the root mean square error of 0.004 and 0.001 g, average r2 of 5-fold cross validation of 0.915 and 0.983, respectively. The results shed light on seismic hazard assessment for critical sites with long-term earthquake recordings.


Introduction
Ground motion models (GMMs) are widely used to predict ground motion intensities across a large region during different earthquakes.Due to the scarcity of earthquake records, GMMs usually assume that the distribution of ground motions over time at a given site is the same as their spatial distribution over all sites for the same magnitude, distance, and site condition [1], which is commonly referred to as the ergodic assumption [2].GMMs usually output a median and a logarithmic standard deviation (generally referred to as sigma) when they are used to predict site-specific ground motions.Sigma is considered aleatory variability associated with ground motion modeling and appears not to have decreased despite a tremendous increase in the number of available records, an advancement of processing tools, and an increased complexity of the predictive models [3].Recent studies have shown that shifting from ergodic GMM to nonergodic GMM can reduce sigma by as much as 40% [1].
Parametric regression approaches are commonly used to develop ground motion models, which require a predefined function form (e.g., [4]).Landwehr et al. [1] developed a nonergodic GMM for California by varying regression coefficients for each source and site coordinate in the dataset.Nonparametric machine learning techniques, which do not assume a fixed model structure, have been applied in ground motion modeling and show comparable or better performance than conventional regression-based ground motion models.For instance, Thomas et al. [5,6] explored the predictability of IOP Publishing doi:10.1088/1755-1315/1334/1/012056 2 the neuro-fuzzy method and support vector method on ground motion predictions using the Pacific Earthquake Engineering Research Center (PEER) NGA database.Derras et al. [7] compared the ability of various site-condition proxies to reduce the aleatory variability of GMMs using the artificial neural network (ANN) method.Withers et al. [8] used ANN to build a ground motion model from a synthetic database of ground motions extracted from the Southern California CyberShake study.Kubo et al. [9] proposed a hybrid approach that jointly uses machine learning and conventional ground motion models to solve the underestimation problem for strong motions due to data bias.
This study developed nonergodic ground motion models using a single-site ground motion dataset consisting of more than one thousand recordings from earthquakes located in different source regions.A nonparametric machine learning technique was employed to build the nonergodic GMMs and easily integrate predictors that are not considered in common ergodic GMMs.The study site, earthquake data, and ground motion data processing were described in Section 2. The modeling strategy and machine learning technique were explained in Section 3, and the results were documented in Section 4. The main findings were summarized in the last section.

Soil profile and geotechnical data
The Onahama Port Array is located in Onahama Port in Iwaki City, Fukushima Prefecture, Japan (GPS coordinates: 36.945°N, 140.912•E).Onahama Port is one of nine strategic commodity ports selected by Japan's Ministry of Land, Infrastructure, Transport and Tourism (MLIT).The Onahama Port is built on shallow reclaimed land overlaid on alluvial sediments [10] (Figure 1).Onahama Port has a longterm seismic monitoring array managed by the Port and Airport Research Institute (PARI) of Japan since 1995 [11].The Onahama Port Array consists of two triaxial accelerometers at the surface and 11 m deep.The site lithology consists of (1) a reclaimed loose surface layer down to 1.3 m depth, (2) a dense sand layer between 1.3 and 7 m in depth, and (3) a silt layer between 7 and 11 m in depth.With the groundwater level at 1.2 m, the loose saturated sandy layer (1.3 to 3.0 m depth) was liquefied during the 2011 Tohoku earthquake [12].The average shear wave velocity of the top 11 m soil ( S11) at the site is 493 m/s.The detailed profile of soil type, SPT N-values, and shear wave velocities are shown in Figure 2. The results of multiple multichannel analysis of surface waves (MASW) [13] and 3-D microtremor horizontal-to-vertical spectral ratios (MHVRs) [14] at the Onahama Port Array suggest a substantial lateral heterogeneity may exhibit at the subsurface and affect the local site response.These complex site conditions can pose challenges for conventional ground motion models in predicting ground motions at the Onahama Port.

Earthquake ground motion data
This study collected all seismic recordings at the Onahama Port Array between 1995 and 2020 (data available from the PARI website).In total, we obtained 1238 pairs of ground motion records whose spatial distribution is shown in

Ground motion data processing
For the ground motion intensity parameters, we compute the 5% damped pseudo-spectral acceleration (Sa) for periods from 0.01 to 10 s following the NGA project [4].The Sa of both downhole and surface ground motions are shown in Figure 5. Here, we only aim to predict the Sa of the surface ground motions.
We select predictors based on considerations of both data availability and physical correlation.The seven candidate predictors are listed in Table 1, along with their group in terms of source, path, and site effects on ground motion modeling.The earthquake magnitude and focal depth are directly obtained from the PARI website, along with the epicentral coordinates and occurrence time of each earthquake.
The epicentral distance and azimuth are calculated based on the coordinates of earthquakes and the Onahama Port Array.The rainfall data was grouped into site response variables as it may influence groundwater activities and then soil properties.We collected hourly rainfall data recorded in Onahama City by the Japan Meteorological Agency.We computed three types of accumulated rainfall intensity to indicate the possible temporal delay of rainfall effects on seismic soil properties.They are rainfalls that accumulate within 1, 3, and 24 hours before the earthquake.rain_day was found to have higher accuracy in predicting Sa than rain_3hour and rain_hour using the methodology described in Section 3. repi * Rainfall accumulated within the 1 hour before the earthquake (mm) rain_hour * Rainfall accumulated within the 3 hours before the earthquake (mm) rain_3hour * Rainfall accumulated within the 24 hours before the earthquake (mm) rain_day *

Modeling strategy
In this study, we aim to predict the Sa for periods ranging from 0.01 to 10 s for a specific site (i.e., the Onahama Port Array).In other words, we aim to predict every curve in Figure 5b using machinelearning methods.There are different strategies to deal with these multiple outputs regression problems.The first strategy is to develop many individual models for different periods.Kim et al. [15] used this strategy to predict period-dependent amplifications based on the KiK-net dataset for 15 selected periods between 0.01 and 7 s.The second strategy is to set all selected spectral accelerations as simultaneous outputs of the model.Derras et al. [16] developed an artificial neural network model outputting 64 variables, including 62 spectral accelerations for 62 periods (ranging from 0.01 to 4 s).This method can consider the interaction of spectral acceleration between the selected periods but requires large computation efforts.We treat this multiple output issue using a different strategy.We put the period as one predictor variable and put the corresponding spectral acceleration as the unique response variable.
In other words, we sample many points (here 1000 points between 0.01 and 10 s with equal intervals of 0.01 s) from each curve in Figure 5b and put the period values to the new predictor column and the corresponding spectral accelerations to the output variable column.Using this strategy, we are able to accurately predict the shape of the spectral acceleration curve and local variations among the whole period range.The effect of the period on the model performance is explained in Section 4.
Figure 5.The 5%-damped pseudo-spectral acceleration (Sa) of downhole and surface ground motions for periods from 0.01 to 10 s.

Gradient boosting method
This work used the gradient boosting method to develop the site-specific ground motion models.A decision tree model can partition the input space into a set of intervals and their Cartesian products and make predictions over these intervals.A single decision tree tends to have an issue of over-fitting and is sensitive to the input data.The gradient boosting method, one of the decision tree ensemble techniques, fits consecutive trees where each new tree reduces the net loss of the prior trees.More details of the decision tree and gradient boosting method can be found in [17].The GB-based site-specific ground motion models can be expressed as: where (b; ) is the b-th decision tree using the explainable variables b and hyperparameters ,  is the weight of the b-th decision tree, B is the total number of decision trees,   is observed ground motion intensity of the i-th earthquake and   is the model residual for the i-th earthquake.We employed the h2o.ai machine learning package to develop the gradient boosting models.The dataset is randomly split into the training and test datasets, consisting of 80% and 20% of the whole dataset.In the training, five-fold cross-validation was implemented, and the hyperparameter tuning for the gradient boosting model was conducted by the grid search function.We used two commonly used regression metrics to evaluate the performance of different GB-GMMs.One is the root mean square error (RMSE), and the other is r 2 .The effects of different combinations of predictors on the gradient-boosting model performance are in Section 4.

Results
Table 2 shows the performance metrics of four GB-GMMs using different predictor combinations.The r 2 averaged from the five cross-validation datasets, r 2 cross-validation, increases when more predictors are fed into the model.Meanwhile, rmse decreases from 0.0035 g to 0.0013 g for the training dataset, and from 0.0039 g to 0.0017 g for the test data, indicating the machine learning models are not overfitted and have a small level of errors.As Model 1 uses predictors widely used in conventional ground motion models and Model 4 adds three new predictors, we further compared the benefits of adding new predictors for IOP Publishing doi:10.1088/1755-1315/1334/1/0120566 improving the accuracy of site-specific ground motion models.Figure 6 compares the measured and predicted Sa in the test dataset using the two ground motion models.The response spectra of all earthquake events estimated using the basic and optimal model are shown in Figure 7a and Figure 7b, respectively.Both Figure 6 and Figure 7 suggest that adding new predictors reduces the errors for ground motion modeling, especially for strong earthquakes.Figure 8 shows the relationship between the Model 4 residuals (i.e.,   in Equation 1) and the six predictors.The results suggest that the machine-learning-based ground motion model has no obvious bias, as no obvious trend between model residuals and predictors was found.Lastly, we checked the variable importance of the six predictors in the gradient-boosting ground motion model.The variable importance is calculated by the variance reduction of each variable to the sum variance reduction of all variables [15].The results indicate that the most influential variable in predicting single-site (i.e., Onahama Port Array site) spectral acceleration prediction is the earthquake magnitude (emag), followed by the period for spectral acceleration computation (period), epicentral distance (repi), earthquake focal depth (evdp), azimuth angle between site and epicenter (azimuth) and accumulated rainfall intensity within 24 hours before the earthquake (rain_day) (Figure 9).

Conclusion
Ground motion modeling is a vital recipe for seismic hazard assessments.This work provided an example of using machine learning techniques for ground motion modeling for a long-terminstrumented site, specifically the Onahama Port Array site in Japan, that has 25-year earthquake recording data.The gradient boosting method was employed to predict 5%-damped response spectra with periods from 0.01 to 10 seconds.Different combinations of predictors are tested and compared.Two gradient boosting models are recommended: the basic model that uses earthquake magnitude and epicentral distance as inputs and the optimal model that uses three additional predictors, namely, focal depth, azimuth, and rainfall.The basic and optimal gradient boosting models have r 2 averaged from the 5-fold cross-validation of 0.915 and 0.983 and root mean square error of 0.004 and 0.001 g, respectively.The results shed light on long-term seismic hazard assessment for critical sites with a certain amount of earthquake recordings.

Figure 1 .Figure 2 .
Figure 1.Geology map of Onahama port area and liquefaction observations of the 2011 Tohoku earthquake (modified from Yamaguchi et al. [10])

Figure 3 .
The JMA magnitude of recorded earthquakes ranges from 2.4 to 9.0.The epicentral distance ranges from 4.3 to 1252.2 km (see Figure4a).The minimum, maximum, and median surface PGA are 0.6, 1793.4,and 4.4 gal, respectively (see Table7.2), and the relation between the PGA and earthquake magnitude is shown in Figure 4b.The 2011 Mw 9.0 Tohoku earthquake induced the strongest ground motion and the only liquefaction observation at the Onahama Port Array during the instrumentation period.Onahama Port Array was inundated by the great tsunami induced by the 2011 Tohoku earthquake, and a new apparatus was installed in 2013 at the same place (personnel communication with Yosuke Nagasaka of PARI).Hence, the station only records nine aftershocks within 5 hours of the Tohoku main shock and has an instrumentation gap between 2011 and 2013.

Figure 3 .
Figure 3. Location map of the Onahama Port Array and its observed earthquake events.The location of the Onahama Port Array is indicated using a red triangle.Epicenters of earthquake events analyzed in this study are shown in large circles.The size of the circle indicates the magnitude of each event, and the color shows the depth, with red being shallow and blue being deep.The base map is the United States Geological Survey (USGS) topographic map.

Figure 4 .
Figure 4. Magnitude, epicentral distance, and PGA distribution of surface ground motions recorded at the Onahama Port Array between 1995 and 2020.

Figure 6 .
Figure 6.Comparison between measured Sa and precited Sa in the test dataset: (a) Model 1 and (b) Model 4.

7 Figure 8 .
Figure 8. Relationship between Model 4 residuals and different predictors.The red dots denote the moving average of residuals.The vertical black bars on the red dots denote the standard deviation.

Table 1 .
Summary of predictors used in the nonergodic ground motion modeling.

Table 2 .
Model performance of the gradient boosting models using different predictor combinations.