Before we start explaining the GB, we describe the decision tree algorithm, which is the main concept of GB and RF. The decision tree consists of nodes, where a tree is grown on the training dataset. The tree contains three types of nodes: root node, internal node, and leaf node, where the root and internal nodes play a role in splitting the data samples, and the leaf node makes the final decision for the prediction value.

We presented an example tree using the independent variables of the K-NET, as shown in Figure 4 to explain the internal structure. First, the root node splits 15,253 independent data samples into two internal nodes by asking if the N-value ≤ 45.5. If the condition is true, the internal node condition (i.e., N-value ≤ 14.75) works to further divide the allocated data samples into a leaf node (R ₁ ) and another internal node. If the root node condition is false, the data samples are further divided by the internal node (i.e., depth ≤ 4.5 m) into a leaf node (R ₄ ) and another internal node. Following the if-else rules, the model finally creates a tree that consists of a root node, four internal nodes, and six leaf nodes (R ₁ , R ₂ , . . ., R ₆ ).

One may wonder how the decision tree model creates the splitting criterion of the node. The model grows a tree by splitting the data samples into two groups by finding the threshold that minimizes the mean of squared errors (MSE), which is calculated as MSE = 1 n ∑ j = 1 J ∑ i = 1 n V S i m e a − V S R j e s t 2 (1) where V S i m e a is the measured V _S associated with i t h data sample of the K-NET (i.e., i = 1, 2, . . ., n; n = 15,253), and V S R j e s t is the estimated V _S determined by a specific leaf node, R _j , where j is the leaf node index (j = 1, 2, . . ., J).

The estimated V _S ( V S e s t ), which is an output of the trained tree model (T), can be expressed by the following equation: V S e s t = T x i = ∑ j = 1 J c j I x i ∈ R j (2) where x i represents the independent variables of the i t h data sample of the K-NET, j and J are the specific leaf node number and the total number of leaf nodes (i.e., six leaf nodes for the example tree in Figure 4), respectively, c j is the predicted dependent variable decided by the specific region R j , and I is an indicator function that takes a value of 0 or 1 (i.e., I = 1 if x ∈ R j and 0 otherwise).

However, a single decision tree model is prone to overfitting on a training dataset, resulting in a high variance in new data samples (Geurts et al., 2009; Czajkowski and Kretowski, 2019). The GB algorithm, proposed by Friedman (2001); Friedman (2002), is an ensemble of weak models (i.e., decision trees) and provides robust model performance over the overfitting problem. GB grows many decision trees and connects them in order like links in a chain, where each new tree is grown to modify a mistake made by a previous tree. An example of a GB architecture is presented in Figure 5.

The trees in the GB estimate the residuals between V S m e a and V S e s t instead of the V _S itself and are trained to minimize the residuals. The specific steps for training the GB model are as follows. Step 1: The GB has a constant value ( F 0 x ), which is V S m e a ¯ from the training dataset. Step 2: From b = 1 to B , where B is the last tree index, the GB repeats the following steps (3–5) for successive trees ( T 1 , T 2 , … , T B ). Step 3: An individual tree T b x calculates the residual data for each data sample as: r i b = V S i m e a − F b − 1 x i (3) where r i b is a residual associated with data sample i in the dataset (i.e., i = 1 , 2 , … , 15,253 ) and tree index b , and F b − 1 x i is the prediction value of the previous GB model ( V S i , b − 1 e s t ). Step 4: After the residual dataset ( r 1 , b , r 2 , b , … , r 15253 , b ) for T b x is developed, T b x is trained on the dataset, x i , r i b i = 1 15253 , instead of x i , V S i m e a i = 1 15253 . The leaf node ( R j b ) is determined during training, where j is the leaf node index of the tree, T b x . The mean residual value predicted in R j b is r j b . The r j b is subsequently reduced by a learning rate ( v ), which is a constant value to reduce the contribution of each tree. Therefore, the tree ( T b x can be described as follows: T b x = v ∑ j = 1 J b r j b I x ∈ R j b (4) where J b is the number of leaf nodes in T b x . The equation returns the v * r j b according to independent data ( x ). Step 5: After a tree is built, it is added to the previous tree. The updated GB model ( F b x ) can be described as follows: F b x = F b − 1 x + T b x = F 0 x + T 1 x + … + T b x (5)

After the GB finishes developing the last tree ( T B x ), we finally obtain the F B x , which is the complete GB model ready to use for predicting V _S ( V S e s t ). In this study, we set the number of trees ( b ) to 100 and the learning rate ( v ) to 0.1, as suggested by Pedregosa et al. (2011).

The RF algorithm, proposed by Breiman (2001), is a bootstrap aggregation (bagging) ensemble algorithm that grows many decision trees using a random subset of the data. Unlike GB, RF trains many weak trees in a parallel manner, where the trees are not affected by each other while being trained. Each tree in the GB predicts the residual value, but the tree in the RF directly returns V S e s t . An example of an RF architecture is presented in Figure 5. From b = 1 to B , each tree ( T 1 , … , T B ) is grown on bootstrap samples ( D 1 , … , D B ), where D b is the randomly sampled subset data from the training dataset.

While T b is trained on D b , the number of variables ( p * in T b is randomly chosen from the total number of independent variables ( p to minimize the MSE. The node is further split into two child nodes after choosing the best split points among p * number of variables.

After RF completes training all trees, it makes a final decision of V S e s t at a new data sample x by averaging the multiple results of T b , which can be described as V S e s t = 1 B ∑ b = 1 B T b x (6)

In this study, we set the number of trees ( b ) to 100 and p * = p (i.e., all variables are considered), as suggested by Pedregosa et al. (2011).

The ANN model comprises a collection of nodes grouped in layers, where each node in a layer is connected to the nodes in the next layer. The ANN model includes three types of layers: input layer, hidden layer, and output layer. Figure 6 presents an ANN model containing the two hidden layers used in this study. The number of input variables is nine for K-NET and seven for KiK-net, as described in Table 1. However, we applied the binary encoding method to categorical variables. The total number of variables was increased to 18 for the K-NET and 22 for the KiK-net to train ML models (i.e., GB, RF, and ANN). A detailed explanation of this is provided in the subsequent section. Therefore, the number of input nodes is 18 for K-NET and 22 for KiK-net. We set the number of nodes in the hidden layers to 200, as inspired by Kim et al. (2020). In Figure 6, the values of each node for hidden layer 1 ( h j 1 ), hidden layer 2 ( h k 2 ), and output layer ( V s e s t ) can be described as h j 1 = f 1 ∑ i = 1 18 w j i x i + b j , h k 2 = f 2 ∑ j = 1 200 w k j h j 1 + b k , V S e s t = ∑ k = 1 200 w n k h k 2 + b n (7) where x i is the input value of the i t h node in the input layer, w j i is the weight between i t h node in the input layer and j t h node in the hidden layer 1, w k j is the weight between j t h node in hidden layer 1 and k t h node in hidden layer 2, and w n k is the weight between k t h node in hidden layer 2 and n t h node in the output layer (i.e., n = 1 ). b j , b k , and b n are the biases of j t h node in hidden layer 1, k t h node in hidden layer 2, and n t h node in the output layer, respectively.

Each node in an input layer receives an independent variable (e.g., the N-value). At each node, the value ( x ) is multiplied by the corresponding weight ( w ), which is summed with other values multiplied by other weights in the same layer. A bias ( b ) is then added to the sum of all values multiplied by each weight in the layer. The bias provides a better generalization ability to the model by enhancing the fitting flexibility. Then, the activation function ( f ) is applied to the sum of x × w + b , where f gives a non-linear property to help the ANN model capture the complex data patterns. The activation function outputs a value that becomes the input value of the node for the next layer.

We applied a rectified linear unit (ReLU) to f 1 and f 2 , where the ReLU is a widely used non-linear activation function (Boob et al., 2020). Specifically, the ReLU is defined as f X = max 0 , X , where f X = 0 if X < 0 or f X = X if X ≥ 0 .

The three ML-based models developed in this study were evaluated for each test fold after training. Figure 7 presents the measured wave velocities (i.e., V P m e a and V S m e a ) versus the estimated values ( V P e s t ; V S e s t ) using the three models for the entire K-NET test data samples. Note that these data were aggregated from five test folds from the five experiments (i.e., five derived ML models). As a performance indicator for V _S prediction, we calculated the root mean squared error ( RMSE V S ) as: RMSE V S   m / s = 1 n ∑ i = 1 n V S i m e a − V S i e s t 2 (8) where V S i e s t is the estimated V S value for the i t h data sample, V S i m e a is the corresponding measurement, and n is the data sample size. The RMSE values calculated for each of the five test folds were averaged. The RMSE V P was computed in the same manner. The RMSE V S and RMSE V P are the smallest for the GB-based model and the largest for the ANN-based model.

Figure 8 presents the measured wave velocities (i.e., V P m e a and V S m e a ) versus the estimated values ( V P e s t ; V S e s t ) using the three models for the KiK-net dataset. RMSE V S and RMSE V P are the smallest for the GB-based model indicating a stronger alignment along the 1:1 line, and the largest for the RF-based model. The results for individual experiments are included in Supplementary Appendix Ⅰ of the Electronic Supplement.

The RMSE depends on the study area and data features including the number of sites and velocities distribution. Many previous studies have utilized varying ranges of V _S to make predictions for different geological regions, resulting in varied RMSEs. For example, Ataee et al. (2019) utilized uncorrected and corrected SPT-N with 88 boreholes to predict V _S . The results using uncorrected SPT-N and V _S under approximately 1,200 m/s presented that the RMSEs of the models ranged from 94.512 to 104.149 m/s. Those using corrected SPT-N and V _S under approximately 600 m/s presented RMSEs ranging from 59.423 to 67.473 m/s. Ghorbani et al. (2012) utilized corrected SPT blow counts, and effective overburden stress to predict V _S . They used 80 boreholes, where the V _S ranges from 66 to 363 m/s. The RMSE of the prediction model is 37.2 m/s. Sun et al. (2013) used tip resistance, sleeve friction, pore pressure, and overburden effective stress to establish the correlation with V _S . They utilized 17 sites, where the measured V _S is under approximately 400 m/s. The RMSEs of the correlation forms are from 30.42 to 38.57 m/s. Furthermore, Dumke and Berndt (2019) used 38 types of variables (e.g., depth below seafloor, surface heat flow, and distance to nearest spreading ridge) to predict V _P . They used 333 sites containing V _P above 4,000 m/s, where the velocity range is not mentioned. The RMSEs vary approximately between 400 and 500 m/s depending on the considered variables. The RMSE presented in this paper may be reasonable, given that the prediction models were made and tested for the larger number of sites distributed throughout Japan, which includes various study areas and a wider range of velocities than other studies. However, discrepancies have been observed, especially for the KiK-net: V _P dataset, implying that more region-specific depth-related variables may be needed to infer the velocity profiles better.

We further investigated the relationship between the measured and estimated velocities by employing the Regression Error Characteristic (REC) curve (Bi and Bennett, 2003). The REC curve depicts the relationship between the specified deviation tolerance on the x-axis, which is the error tolerance, and the y-axis for the proportion of data with prediction deviations smaller than the corresponding deviation. The resulting curve provides an estimation of the cumulative distribution function of the error. Furthermore, the REC curve quantifies the performance of the model by computing the area under the curve (AUC). A higher AUC value indicates better model performance. Figure 9 illustrates the REC curves for each model. The curves were individually computed for the five test folds from the five experiments and were then averaged to make a single curve. The AUC was subsequently calculated based on the single curve, representing the general performance of each model on the dataset. The results for K-NET and KiK-net, both at V _P and V _S , reveal that all AUCs are above 0.702 for the specified deviations ranging from 0 to 1.0. Notably, the GB-based model has the highest AUC across all cases, indicating its relatively strong predictive performance within the deviation range.

Figure 10 presents maps of RMSEs for the GB-based model for all the sites considered in this study, which were aggregated from five test folds from the five experiments (i.e., five derived ML models). Overall, the models for V _P and V _S of the K-NET sites (Figures 10A, B) indicate that almost 80% of the sites have RMSE values within the range of (0, 500] for V _P , and almost 99% are within the same range for V _S . In contrast, for KiK-net (Figures 10C, D), approximately 46% of the sites have RMSE values within the (500, 1,000] range for V _P , and about 63% are within the (0, 500] range for V _S . It can be noticed that the RMSE values larger than 1,000 m/s for the estimated V _P values at the K-NET sites are concentrated in the area around 139 °E and 35.5 °N (Figure 10A). Furthermore, the RMSE values greater than 1,500 m/s and 1,000 m/s for the estimated V _P and V _S values, respectively, at the KiK-net sites are mainly clustered in the region around 137 °E and 35 °N (Figures 10C, D). The RMSEs for the KiK-net sites show a certain pattern along the east coast (from 140 to 142 °E and from 36 to 40 °N) (see Figure 10D). These observations imply that there could be factors that can affect the V _S and V _P values, other than those considered in this study.

Figure 11 shows examples of the wave velocity profiles predicted by the GB-based model compared with the measured profiles at the nine K-NET sites. The eight and one sample profiles were randomly selected from the V _S RMSE bands of (0, 500] m/s and (500, 1,000] m/s, respectively, from the entire test folds from five experiments. The wave velocities predicted for the HKD024, FKO015, and KNG008 sites (Figures 11A–C, respectively) are in good agreement with the measured profiles when compared to the other illustration, producing RMSE values ≤186 m/s. In contrast, there are some discrepancies between the measured and predicted profiles at certain depth ranges for the rest of the sites. At GNM014, V _P is overestimated at depths of up to 9 m (Figure 11D). At NGN025, the wave velocities are underestimated at depths greater than 13 m (Figure 11E). At EHM012, V _P is underestimated at depths from 2 to 6 m (Figure 11F). There are also some discrepancies in the V _P profiles at the YMT006, IBR009, and KGW008 sites (Figures 11G–I, respectively). In detail, the V _P is consistently overestimated across the entire depth range at the YMT006 site (Figure 11G). At the IBR009 site, V _P is underestimated up to a depth of 5 m (Figure 11H). Similarly, the KGW008 site demonstrates underestimation up to 3 m and overestimation at depths beyond 10 m (Figure 11I).

Figure 12 presents examples of the wave velocity profiles estimated by the GB-based model compared with the measured profiles at the nine KiK-net sites. The six, two, and one sample profiles were randomly selected from the V _S RMSE bands of (0, 500], (500, 1,000], and (1,500, 2000], respectively, from the entire test folds from five experiments. The estimated wave velocities for the HRSH15, KOCH11, IBRH12, and TCGH15 sites (Figures 12A–C, E, respectively) comparatively match well with the measured profiles, producing RMSE values ≤388 m/s. Some discrepancies are observed at specific depth ranges for the other sites. At IBRH07, the wave velocities are overestimated at depths from 51 m to 650 m and underestimated at depths from 651 m to 1,050 m (Figure 12D); however, the estimated wave velocities show relatively close agreement beyond 1,050 m. At SITH03, the velocities are overestimated almost throughout the depth (Figure 12F). At AICH12, the velocities are underestimated at depths greater than 49 m (Figure 12G). Overestimations are observed at the ISKH06 site for depths greater than 44 m (Figure 12H). For the YMTH08 site, velocities are overestimated for depths exceeding 18 m, while the V _P is underestimated for depths up to 6 m (Figure 12I).

Some discrepancies were observed in relation to a certain depth and profile patterns, as shown in Figure 11 and Figure 12. It is possible that the model could not well predict velocities for sites that have unusual profile patterns or for those that have not been frequently used when training or are not included in the training dataset. The model was trained to reduce the overall error for the entire sites used in training, so it might not well generalize the unseen patterns. As seen in the samples in Figure 11 and Figure 12, the model was trained to predict slower velocities near the ground surface and faster velocities at greater depths. Furthermore, model predicts velocities gradually increasing, and does not predict well the abrupt velocity changes (e.g., Figure 11E; Figure 12G). There are velocity reversals at depths greater than 44 m at ISKH06 (Figure 12H), and the V _P values are unusually faster near the ground surface at YMTH08 (Figure 12I). It turned out that the model was not able to capture these profiles.

For the systematic evaluation of discrepancies between measured and estimated wave velocities, we computed the residuals for all the continuous variables as Res V S = ln V S m e a − ln V S e s t Res V P = ln V P m e a − ln V P e s t (9) where Res V S and Res V P are the residuals for V _S and V _P , respectively. We also calculated the standard deviations of the residuals and biases (i.e., mean values of the residuals) for all trained models.

Figure 13 shows the Res V S of the GB-based model for the K-NET stations, which were aggregated from five test folds from the five experiments. The Res V S ranges from −1.798 to 1.677 and is aligned along zeros with respect to all continuous variables with a small bias (i.e., −0.066). Figure 14 shows the Res V P of the GB-based model for the K-NET stations, which were aggregated from five test folds from the five experiments. The Res V P ranges from −1.765 to 1.396 and is also aligned along with zeros with a bias value of −0.059. The standard deviation values for Res V S are 0.354, 0.358, and 0.374 for the GB-, RF-, and ANN-based models, respectively. The standard deviation values for Res V P are 0.360, 0.361, and 0.383 for the GB-, RF-, and ANN-based models, respectively. The residuals do not show any trend with the considered variables, indicating that all the variables have certain contributions to the model, or that some of the variables do not have influence on wave velocities. Moreover, the results seem reasonable when compared to those of Kwak et al. (2015), who presented the range of standard deviation of residuals for V _S prediction models made using K-NET for each soil/rock type, which was from 0.245 to 0.462. Because the Res V S and Res V P of the GB-based model for KiK-net stations show the same aspects, we describe the results in Supplementary Appendix Ⅱ of the Electronic Supplement.

We examined the contribution levels of the independent variables to the prediction accuracy of the best model, the GB-based model. The method is called variable importance (VI), which is computed as the sum of the decrease in error when a variable splits a tree node (e.g., a node split by an N-value ≤ 14.75). The variable importance for variable x (i.e., VI(x)) is calculated as follows: VI x = 1 B ∑ b = 1 B ∑ t ∈ T b , v s t = x p t ∆ i t (10) where b is the tree index ( b = 1 to B ), t is the node in a specific tree model (T _b ), and v s t is the variable used for splitting the node in which s t is the splitting criterion ( s ) at note t. p t is the proportion ( N t N ) of data samples reaching t , where N is the number of total training data samples, and N t is the number of data samples at node t . ∆ i t is the impurity reduction at node t , which can be expressed as follows: ∆ i t = i t − N t l N t i t l − N t r N t i t r (11) where i t is the MSE at node t , i t l and i t r are the MSEs at the left child node t l and right child node t r , respectively, split from node t . N t l and N t r are the numbers of data samples at t l and t r , respectively.

Figures 15A, B show the relative VIs for the K-NET independent variables for the GB-based models. The relative VI was computed by VI for each variable divided by the total VI for all variables. The VI was calculated on each test fold, and the VIs for all the five test folds were averaged. The VIs computed for binary codes were summed for the categorical variables. Three depth-dependent variables (i.e., depth, N-value, and density) have the highest VIs for both V _P and V _S models. The depth is ranked at the top for the V _P model, whereas the N-value is ranked at the top for the V _S model. Figures 15C, D present the relative VIs for the KiK-net dataset. The depth turned out to be the most critical variable for both V _P and V _S models. The effect of the site location is more significant for the KiK-net model than for the K-NET model. The slope angle and elevation have a certain influence on the models, whereas the soil/rock type and geology have the least influence. Although the influence of the geology turned out to be insignificant, the performance of the GB-based model was enhanced by including it. The RMSEs of the model were reduced from 615 m/s to 597 m/s for V _S , and from 979 m/s to 961 m/s for V _P , implying that it is also related to wave velocities at a deeper depth.

The confining pressure increases with depth, leading to an increase in the density, N-value, and wave velocities. Therefore, the depth and associated variables were determined to be most strongly correlated, as revealed by VI. The slope and elevation are related with shear stiffnesses, which eventually affect wave velocities. The site coordinates are relatively high VI, implying that they may be associated with site conditions that were not captured by other variables. The geology has the lowest VI, as it is for the ground surface. However, we included it in the model because of its certain effect in enhancing the predictive performance.