Predicting permanent slope displacements due to multidirectional earthquake shaking from unidirectional shaking analyses

The current state-of-practice for time-domain modelling of seismic slope stability is to apply a single horizontal ground motion component in the downslope direction. However, several numerical and experimental studies have demonstrated the importance of multidirectional shaking on seismically induced permanent displacements. To understand better the effect of multidirectional shaking on seismically induced permanent slope displacements, we performed 28,100 three-dimensional finite element slope stability analyses using 230 ground motion record pairs, 48 soil profiles, and four different combinations of ground motion components. We then performed traditional regression analyses and applied machine learning techniques to develop simplified models that estimate the permanent displacements due to multidirectional shaking from the results of analyses where only one component of the ground motion is applied downslope. The simplified model has R 2 = 0.91 and uses the displacement from a unidirectional shaking analysis as its only input variable. The standard deviation decreases as the predicted unidirectional shaking displacements increase.


Introduction
Earthquake induced landslides pose a constant and major risk to many areas of the world.One of the deadliest examples is the magnitude 8.0 Haiyuan earthquake of 1920, which triggered more than 7000 landslides in north-central China that created 126 dammed lakes and killed 32,000 people [1].A more recent example is the April 25, 2015 Gorkha, Nepal earthquake sequence, which caused thousands of landslides that killed hundreds of people and blocked vital road and lifeline routes to many villages [2].
There are numerous methods to predict displacements of slopes due to earthquake shaking.These methods vary in complexity from equations based on simplified slope and ground motion parameters (e.g.Bray et al. [3], Bray and Macedo [4]), to Newmark sliding block models that incorporate the full acceleration time series but assume that the failing soil acts as a rigid body [5], to two-dimensional (2D) or three-dimensional (3D) numerical finite-element or finite-difference simulations that incorporate the full slope geometry and use an acceleration time series to represent earthquake shaking.However, common for all these methods is that the current state-of-practice uses only a single horizontal component of a selected representative ground motion applied in the downslope direction [6] to simulate the slope stability during strong shaking.This approach is founded on the basic assumption that the earthquake excitations in the vertical and cross-slope directions are negligible, which substantially reduces the complexity of slope stability analyses (e.g. by permitting analyses in 2D rather than 3D) and allows sensitivity studies with only moderate computational efforts.However, Carlton and Kaynia [7] and Kayen [8] found this assumption to be unconservative and to underestimate the seismic slope deformations in many cases.Therefore, it is important to understand the effect of multidirectional shaking on predicted slope displacements to ensure safe design of slopes and infrastructure built on or near slopes.
Much of the research related to multidirectional shaking is primarily concerned with liquefaction analyses or site response analyses for level ground conditions.Several laboratory [9][10][11], centrifuge [12,13] and numerical studies [14,15] focusing on liquefaction have shown that multidirectional shaking causes an increase in settlements and excess pore pressures compared to unidirectional shaking.However, the magnitude of the difference depends on many factors, such as the intensity of the shaking and density of the soil [16].Several researchers have also investigated the effect of multidirectional shaking on site response analyses for level ground [17,18].These studies found that using two components of a ground motion gave results that more closely matched the measured response in the field.This is due in part to the redistribution of seismic wave energy when movement can occur in multiple directions simultaneously [17].
Despite the large body of research related to multidirectional shaking for liquefaction and site response analyses, there is relatively little research related to the effect of multidirectional shaking on slope stability.Rutherford [19] developed a multi-directional direct simple shear device and found that shearing in a figure eight pattern developed more permanent deformations in fewer cycles than unidirectional cyclic shearing, and that the shear strains tended to accumulate in the same direction as the initial horizontal shear stresses (i.e.downslope), independent of the orientation of the cyclic loading.Anantanavanich et al. [20] used the advanced constitutive model MSimpleDSS (Anantanavanich et al. [21]) to perform nonlinear dynamic slope stability analyses for two generic offshore soft clay sites with depths of 20 m and 100 m and four ground motion pairs.They found that multidirectional shaking predicted about a 20%-40% increase in permanent displacements at the soil surface compared to unidirectional shaking.Carlton and Kaynia [7] conducted 27 three-dimensional seismic slope stability analyses in PLAXIS 3D, applying one, two, and three ground motion components.They found that accounting for two horizontal components rather than just one increases the predicted total displacements on the slope by 25%-50%, but that adding the vertical component has a negligible effect.Kayen [8] performed analyses for gently sloping marine slopes using a modified compliant multidirectional Newmark sliding block model.He found that adding the second component across the slope always gave a larger and more realistic shear stress vector and increased downslope displacements.
One argument for ignoring multidirectional shaking in 2D slope stability analyses is that 2D analyses are inherently more conservative than 3D analyses due to the constrained geometry, and therefore the increase in displacements due to multidirectional shaking will be cancelled out by the decrease in displacements due to 3D geometry effects.This assumption is based on the observation that the static factors of safety computed from 3D analyses are almost always higher than from 2D analyses [22].However, for dynamic analyses, this is not always the case.Ferrari [23] performed 3D seismic slope stability analyses for a slope dipping in two directions and compared the displacements to analyses of 2D cross sections of the slope.He found that the 3D analyses predicted larger permanent displacements than the 2D analyses due to the additional static shear stresses acting in the out of plane direction.Even for slopes dipping in only one direction 2D analyses cannot always be assumed to be more conservative than 3D analyses.For example, Azizian and Popescu [24] performed 3D seismic slope stability analyses by extruding 2D models in the out-of-plane direction to defined widths and applying fixed boundary conditions to the sides of the slope.They found that the difference in predicted permanent displacements at the slope crest depends on the width-to-height ratio (W/H) of the slope, with a reduction of less than 15% for 3 ≤ W/H ≤ 5, and negligible differences for W/H ≥ 5. Therefore, the increase in displacements due to multidirectional shaking cannot be assumed to be cancelled out by a reduction in displacements due to 3D geometry effects, and both should be considered in dynamic slope stability analyses.
The objective of this research is to understand better the effect of horizontal multidirectional shaking on slope stability and how it can be incorporated in a simple yet effective manner to modify the results from traditional unidirectional shaking slope stability analyses.To address this need, we performed 28,100 three-dimensional (3D) slope stability analyses using 230 ground motion record pairs, 48 soil profiles, and four different combinations of ground motion components.We then performed regression analyses and used a suite of machine learning techniques to estimate the permanent displacements due to multidirectional shaking applied in the downslope and across slope directions from traditional analysis results where only one ground motion component is applied downslope.The resulting models can then be used to modify the permanent displacements computed by traditional unidirectional slope stability analyses to account for multidirectional shaking.

Analysis model
We performed three-dimensional (3D) slope stability analyses in the time domain using the finite element program OpenSees [25].We used a modified version of a model originally developed by McGann and Arduino [26] to simulate 3D site response analysis of a system representing an infinitely wide uniform slope.The model consists of a column of 3D single integration point cubic brick elements with 1 m × 1 m x 1 m sides supported vertically at the base (Fig. 1).Tied boundary conditions [27] are applied in the two horizontal directions (x and y) at each elevation to constrain the column to deform as a shear beam.In the x and y directions at the base, we applied earthquake loading using the method of Joyner and Chen [28] and used viscous dampers (dashpots) to simulate the radiation of energy into the underlying infinite elastic half-space [29].In this way, the model is capable of exactly reproducing vertically propagating shear waves and pressure waves without any unwanted reflections from the boundaries.
To model sloping ground conditions in the column, we applied inclined gravity loads with components in the vertical (z) and downslope (x) directions according to the desired slope angle.This generates static shear stresses in the column in the downslope (positive x) direction that simulate an infinitely wide slope.We modelled slopes of finite height by adjusting the depth of the inclined gravity loads and using purely vertical gravity loads below this depth.The advantage of modelling the slope as a single column is that it requires very few elements compared to a "full" 3D model.As a result, the model is computationally extremely fast, with typical run times for a single analysis on an ordinary laptop computer between 10 and 100 s depending on the length of the applied ground motion.

Soil properties and material model
We modelled soil profiles 100 m deep with slope angles of 5, 10 and 15 degrees, slope heights (depths with static shear stresses) of 30 m and 100 m, dynamic shear strength ratios of SR = τ/σ ′ v = 0.2, 0.3 and 0.6, where τ is the dynamic shear strength and σ ′ v is the vertical effective stress, and shear wave velocity profiles for NEHRP (National Earthquake Hazards Reduction Program) site classes C, D and E according to the correlations of Carlton and Tokimatsu [30].Fig. 2 shows the variations in dynamic shear strength and shear wave velocity with depth.The total unit weight of the soil is constant at 19 kN/m 3 , and the water table is at the soil surface for all sites.The 100 m deep soil profile is situated on top of an infinitely deep elastic half-space (i.e. a compliant base).For sites with NEHRP D and E shear wave velocity profiles, the underlying elastic half-space had a shear wave velocity of 760 m/s, and for the NERHP C sites it had a shear wave velocity of 1211 m/s, the latter corresponding to the shear wave velocity at the bottom of the soil profile at 100 m depth.These geometry and soil property combinations result in 54 unique sites.However, six combinations that included the lowest strength profile (τ/σ ′ v = 0.2) and the largest slope angle (15 degrees) were not statically stable and were therefore excluded from the analyses, resulting in 48 unique sites.
For all sites and soil layers we used the hyperbolic elasto-plastic material model PressureIndependMultiYield (PIMY) in OpenSees [31,32].PIMY models plasticity using the multi-surface concept with an associative flow rule and Von Mises type yield surfaces.Plasticity is only modelled due to deviatoric stresses and strains, and volumetric response is linear elastic and independent of the deviatoric response.We assumed a failure strain of 1% and simulated small strain damping using conventional Rayleigh damping with a target damping ratio of 3% specified at frequencies of 1 Hz and 5 Hz.

Acceleration time series
We selected 230 ground motions from the PEER NGA West 2 online database [33].Each ground motion consists of two orthogonal horizontal acceleration time series (component A and B) and a vertical component (not used).The ground motions have time averaged shear wave velocities in the top 30 m (Vs 30 ) between 500 m/s and 1500 m/s to be consistent with the shear wave velocity of the underlying half-spaces used in the finite element model, peak ground accelerations (PGA) greater than 0.05 g for at least one component, and no pulse-like characteristics.Fig. 3 shows the moment magnitude (M w ), Joyner-Boore distance (R jb ), and PGA of both components of the ground motions.
One element of uncertainty in dynamic slope stability analyses is the orientation of the ground motion.When performing dynamic slope stability analyses, engineers do not know a priori which component or orientation of the design ground motion to apply downslope.Carlton et al. [34] found that for one ground motion, the estimated seismic slope displacements varied by 20-100% depending on the orientation of the ground motion and slope angle.While these findings were based on a small sample size, they illustrate the potential uncertainty due to ground motion orientation on predicted seismic slope displacements.Bray et al. [3] and Bray and Macedo [4] performed unidirectional slope stability analyses and accounted for ground motion orientation uncertainty by applying each component individually in its given orientation and in its opposite polarity.They then used the maximum of the two averages of the components in the different orientations in their regression analyses.In an attempt to capture some of the uncertainty related to ground motion orientation, we performed four analyses for each of the 230 ground motions: (1) component A applied in the x-direction (downslope); (2) component B applied in the x-direction; (3) component A in the x-direction (downslope) and B in the y-direction (across slope); and (4) component B in the x-direction and A in the y-direction.Based on the results of Carlton and Kaynia [7], who found that adding the vertical component had a negligible effect on the permanent displacements, we did not include the vertical component in any of the analyses.

Verification of analysis model
Carlton et al. [35] verified the OpenSees model described above against dynamic time domain simulations in PLAXIS 3D and FLAC 3D using actual slope geometries.PLAXIS 3D and FLAC 3D are commercial finite element and finite difference programs, respectively, that are widely used in geotechnical engineering practice.The slope models in PLAXIS 3D and FLAC 3D were 250 m wide, 300 m long, 68.4 m tall at the top of the slope, and 55 m tall at the foot of the slope.The slope itself was 13.4 m tall and 52 m long, with a slope angle of 15 degrees.Carlton et al. [35] found that all three programs (OpenSees, PLAXIS 3D, and FLAC 3D) gave similar predictions of permanent displacements due to earthquake shaking, and that none of the programs predicted consistently higher or lower permanent displacements.
We further verified the OpenSees model by comparing the results of unidirectional shaking with results from PLAXIS 2D for several analysis cases.We computed the response for a 60 m high slope with slope angles of 5, 10 and 15 degrees, and for a 100 m high slope with a slope angle of 10 degrees.All analyses used the same ground motion recorded at the Wonderland station during the 1994 Northridge earthquake (M w = 6.7,R jb = 15 km, PEER NGA West 2 record sequence number = 1011).For computational reasons only the first 10 s of the response was computed; this contains more than 95% of the Arias Intensity of the record.In PLAXIS 2D, we used the NGIADP material model [36], which has a hyperbolic stress-strain shape similar to the PIMY material model.The same soil and material input parameters were used for both programs, and isotropic behaviour was assumed in the NGIADP model.
The results presented in Fig. 4 show a reasonable agreement between the two computational models.There are some discrepancies, however, these are primarily due to inherent differences between the calculation models: (i) the NGIADP material model in PLAXIS 2D does not give identical stress-strain response as the PIMY material model used in OpenSees, and (ii) the way soil layering must be implemented in PLAXIS 2D leads to somewhat different vertical discretization of the soil profiles.Nevertheless, these results together with the verification presented in Carlton et al. [35] demonstrate the appropriateness of the model for seismic slope deformation analyses.As a result, we chose to use the OpenSees model in this work, due to the faster computation time, which allowed us to investigate a larger range of soil properties and ground motions.

Displacement database characteristics
We performed 44,160 3D slope stability analyses in OpenSees using the 48 sites, 230 ground motions, and four analysis cases described above.Table 1 lists each of the parameters varied in the analyses and their values.For approximately 1250 of the dynamic analyses the implicit numerical time-integration algorithm failed to converge and these analyses were discarded.To prevent analyses with negligible displacements from dominating the regression analysis, we also removed all analyses with maximum permanent displacements less than 0.5 cm and the corresponding other three analysis cases.For example, if case 1 (only  component A applied downslope) had a maximum permanent displacement less than 0.5 cm, then we excluded it as well as the corresponding cases 2-4.This resulted in a final database of 28,100 analyses, with 7025 analyses for each of the four cases.While the threshold permanent displacement that can be considered negligible would depend on the specifics of a given project, we chose a threshold permanent displacement value of 0.5 cm similar to Bray et al. [3] and Bray and Macedo [4], which is small enough for most practical applications.
We define the maximum permanent displacements due to seismic shaking as u 1 when resulting from analyses with only a single component (i.e.case 1 and case 2) and as u 2 when resulting from multidirectional shaking (i.e.case 3 and case 4).We define the ratio between the maximum permanent displacements as the multidirectional factor (MDF), where MDF = u 2 /u 1 .Fig. 5 shows histograms and lists the mean, median and standard deviations of u 1 , u 2 , and MDF for the final database (14,050 values of u 1 , u 2 , and MDF each).Fig. 5 shows that u 1 , u 2 , and MDF are roughly log-normally distributed, with mean values of 3.4 cm, 6.0 cm, and 1.76, median values of 2.74 cm, 5.04 cm, and 1.64, and maximum values of 281 cm, 393 cm, and 10.98, respectively.
The value of MDF is always greater than 1.0, illustrating that adding a second ground motion component results in larger expected permanent displacements than when only one component is used.This is most likely because by adding the earthquake component in the cross-slope direction the total energy applied to the slope increases and the loading pattern is more damaging to the soil.These results agree with the multi-directional direct simple shear tests of Rutherford [19], who found that shearing in a figure eight pattern increased permanent displacements compared to unidirectional cyclic shearing due to the combination of horizontal shear stress rotation as well as complete shear stress removal and reversal during loading.

Development of regression models
We developed regression models to estimate u 2 from u 1 and simplified ground motion and site parameters using two different approaches.Section 3.1 describes the first approach using least squares regression.Section 3.2 outlines the second approach that utilises machine learning techniques.Both approaches are applied to the same dataset described in Section 2.4.

Least squares regression
We investigated the influence of 13 different predictor parameters to predict u 2 .Table 2 lists the parameters and gives a brief description of each.We investigated four site parameters, eight ground motion parameters, and u 1 .Each of the site and ground motion parameters are readily available when conducting a seismic site response analysis, whereas u 1 can be calculated using one of the methods for predicting unidirectional seismic slope deformations presented in Section 1.The parameter with the greatest predicting power was u 1 , which is not surprising because u 1 is itself a result of an analysis that takes into account all of the other ground motion and site parameters considered.
Fig. 6 shows u 1 versus u 2 and u 1 versus MDF.As u 1 increases, u 2 increases and MDF decreases.In other words, for larger displacements, the ratio between predicted displacements from unidirectional and multidirectional shaking decreases.The effect of multidirectional shaking may be larger for smaller displacements because when u 1 is large the soil is already yielding and adding a second component in a perpendicular direction only increases the displacements by a minimal amount.In contrast, when u 1 is small, the applied cyclic shear stress in the downslope direction may not be large enough to overcome the shear strength and cause yielding, but when adding stresses also in the perpendicular direction the combined stress may exceed the yield limit.
We investigated a variety of functional forms including linear, log, natural log, exponential, power, and second and third order polynomials.We tried all the functional forms for each of the parameters individually against u 2 and the natural log of u 2 .Equation (1) presents Fig. 7 shows the residuals from equation ( 1) versus the 12 remaining predictor parameters, where we define the residuals as ln(u 2 ) from the numerical simulations minus ln(u 2 ) predicted from equation ( 1).The orange lines and error bars represent the mean and standard deviation for ten equally spaced bins for the ground motion parameters, and one bin for each unique value for the site parameters.Fig. 7 shows that there is no clear trend with any of the predictor parameters except possibly D 5- 95 .To determine if the increase in goodness of fit outweighed the increased complexity of adding another predictor parameter, we used Akaike information criteria (AIC) and Bayesian information criteria (BIC).Both measures apply scores to models based on their goodness of fit and penalties for a larger number of predictor parameters to prevent overfitting of models.AIC is based on information theory whereas BIC is based on Bayesian theory.When trying each of the functional forms described earlier with each of the 12 remaining predictor parameters against the residuals of equation ( 1), neither the AIC nor BIC metric supported adding a second predictor variable.This was further confirmed by the fact that adding a second or even third predictor variable did not increase R 2 by more than 0.01 or reduce the standard deviation by more than 0.02 natural log units.
Fig. 8 shows the predicted u 2 using equation ( 1) versus the value of u 2 from the numerical simulations, and u 1 versus the residuals of equation (1).The R 2 of equation ( 1) is 0.91, indicating a good fit to the data.The mean and standard deviations of the residuals for ten equally spaced bins, shown as the orange lines in Fig. 8b, show that the mean values are centred on zero, indicating no model bias with u 1 .However, the standard deviation tends to decrease with increasing u 1 .As a result, we fit equation ( 2) to the standard deviations of the ten bins: Equation ( 2) has an R 2 = 0.98 and u 1 is in cm.The standard deviation decreases as u 1 increases most likely for the same reason MDF decreases as u 1 increases.For non-pulse like motions, Bray and Macedo [4] found σ ln = 0.72 for slope displacements due to earthquakes from active crustal regions, and Bray et al. [3] found σ ln = 0.73 for slope displacements due to earthquakes from subduction zones.These values are larger than the maximum standard deviation from equation ( 2) that ranges from 0.45 for u 1 = 0.5 cm to 0.125 for u 1 = 300 cm.The reason for the smaller standard deviations found in this paper is most likely because equation (1) includes u 1 , which is a very strong predictor of u 2 and provides a better fit than when only using simplified ground motion and site parameters directly.

Machine learning approach
For the machine learning approach, we used the same database as described in Section 2.4 and investigated the same predictor parameters as listed in Table 2. Machine learning (ML) is an umbrella term that covers a range of empirical approaches for both regression and classification (supervised or unsupervised) analysis of nonlinear systems.Such systems can be massively multivariate involving anything from a few predictor parameters to more than one thousand.There are a wide range of regression algorithms available for ML.Common for the different ML approaches is their ability to use computational methods to automatically "learn" and improve directly from the data without being explicitly programmed.The algorithms adaptively improve their performance as the number of data points available for learning increases.
We first split the data into overall training and testing datasets.The training dataset contains the data the ML algorithm will use to tune hyperparameters and evaluate the model, while the testing dataset is used for final validation of the tuned models using previously unseen data.This helps to prevent overfitting of the model.We used a randomized 80/20 training/testing split to ensure balance between sufficient training data for optimizing the ML models, and sufficient data for validation of the final models.
Next, we evaluated the performance of each ML model using repeated k-fold cross validation with five splits and three repeats.This method repeats the cross-validation procedure multiple times and reports the scores from all runs as output.For each model, the mean and standard deviation of the R 2 are used as the performance metric to evaluate the goodness of fit of the models.We then performed model hyperparameter tuning using a randomized search cross-validation algorithm that randomly samples points in a bounded domain of input hyperparameters and uses cross validation to evaluate the model performance for each sampled combination of hyperparameters.Finally, we evaluated each model with the hold-out testing dataset, using R 2 to measure goodness of fit.
The ML regression analyses are implemented in the programming language Python using primarily the scikit-learn (also known as sklearn) machine learning library [37].Table 3 lists the six different ML algorithms that we tested.Fig. 9 shows the results from cross-validation model screening, where the mean and standard deviation of R 2 are plotted for each model using the full set of input parameters in Table 2 to train the models.All models perform well, with the ensemble and boostrap models (Bagging, Fig. 6.Predicted maximum permanent displacements due to all analyses with only a single component (u 1 ) versus a) maximum permanent displacements for all analyses with two components (u 2 ) and b) the multidirectional factor (MDF = u 2 /u 1 ).The orange line is the best fit line using equation (1).
Random Forest and Extreme Gradient Boosting) tending to perform the best, with R 2 values of 0.96 or higher, and XG Boost with the highest score of all with R 2 = 0.985.However, the multivariate linear least squares regression model also has a high R 2 value of 0.925.Part of the explanation for this can be found in Fig. 10, which shows the permutation feature importance of the input parameters (Table 2) used in the Extreme Gradient Boosting method.The permutation feature importance is defined to be the decrease in the model score when a single feature variable is randomly shuffled [38].This procedure breaks the relationship between feature and the target so a drop in the model score is indicative of how much the model depends on the given feature.Fig. 10 shows that the input variable u 1 is the most important input feature by a significant margin, indicating that a model that is primarily a function of u 1 may be adequate to accurately predict u 2 .This is further demonstrated in Fig. 11 where the mean and standard deviation of R 2 from cross validation are plotted for each model but using only u 1 as the input parameter to train the models.Fig. 11 shows that the performance scores are only slightly lower than when using the whole set of input parameters in the model training (e.g. for linear least squares regression the R 2 is now 0.91 vs. 0.925).
Finally, Fig. 12 shows the u 2 -predictions by the Extreme Gradient Boosting model using the full set of input parameters in Table 2, applied to the testing dataset, i.e., on the hold-out data not "seen" by the model in the training phase.The performance scores of the models on the Fig. 7. Residuals (in natural log units) versus parameters investigated (see Table 2 for a description of the parameters).The orange lines and error bars represent the mean and standard deviation for ten equally spaced bins for the ground motion parameters, and one bin for each unique value for the site parameters.testing dataset are in-line with the training scores (here approx.0.99), indicating that the models are not overfitted to the training data.

Limitations and applications
The intended use of the models developed in Section 3 are to give a quick, first order approximation of the increase in seismically induced permanent displacements due to multidirectional shaking based on the results of a traditional unidirectional seismic slope stability analysis.They are not intended to replace more detailed analyses for critical structures.We attempted to use a large enough parameter space that the multidirectional models should be applicable for most practical engineering design purposes.However, because the multidirectional models are empirical, they should not be extrapolated outside the range of parameters used to develop them.The valid parameter range of the multidirectional models is 0.5 cm ≤ u 1 ≤ 300 cm, 3.9 ≤ M w ≤ 7.6, 10 ≤ R jb ≤ 155 km, 0.2 ≤ SR ≤ 0.6, 5 • ≤ slope angle ≤15 • , and NEHRP site class C, D, and E. We also only used ground motions from active crustal regions, half-space properties representative of active crustal regions, and no pulse like motions.Therefore, the models should only be used with extreme caution in stable continental and subduction zone regions, and should not be used at all for ground motions with pulse like

Table 3
Machine Learning models investigated.

Linear regression
Multivariable linear regression that attempts to model the relationship between input and output variables by fitting a linear equation to the observed data.

Ridge Regression
An extension of linear regression that adds penalties to the loss function during training that encourages simpler models.

Bootstrap Aggregation Regression (Bagging)
An ensemble ML algorithm that combines many decision tree models to generate a higher performing model.

Random Forest Regression
An extension of bootstrap aggregation of decision trees.One of the most popular and widely used ML algorithms given its excellent performance across a wide range of modelling problems.

Multi-Layer Perceptron (MLP)
A "feed-forward" artificial neural network (ANN) that consists of multiple layers where each layer is fully connected to the following one through nonlinear activation functions.

Extreme Gradient Boosting
XGBoost is a gradient boosting algorithm that uses gradient descent optimization where trees are added one at a time to the ensemble to correct the prediction errors made by prior models.
Fig. 9. Performance metrics for different ML regression models for predicting the displacement u 2 when using the full set of input parameters listed in Table 2.The blue dots are individual results from cross-validation, and the orange lines and error bars represent the mean and standard deviation for each of the model results.To develop the database of seismically induced permanent displacements we used the hyperbolic elasto-plastic PIMY model in the finite element program OpenSees.The PIMY model does not consider excess pore pressure generation or strain softening.Therefore, the multidirectional models should not be used for liquefiable soils or soils with strain softening behaviour that could lose significant strength during or after earthquake shaking.More advanced constitutive models for estimating displacements due to multidirectional earthquake shaking exist (e.g.Anantanavanich et al. [21]; Yang et al. [39]).However, these constitutive models are not conducive to parametric studies, and the objective of this study was to develop a model that is applicable for a wide enough parameter space to be useful in practice.The method we used is similar to the state of the practice and is therefore consistent with how the multidirectional models are intended to be used.In addition, when comparing to similar types of parametric studies, the numerical and constitutive model used in this study are more realistic.For example, Newmark sliding block analyses or displacement models based on them (e.g.Jibson [40], Saygili and Rathje [41]) treat the sliding mass as a rigid block and neglect the dynamic response of the system.Therefore, these models are only valid for shallow stiff soil with a short fundamental period.The models of Bray et al. [3], Bray and Macedo [4] use a fully coupled, nonlinear, deformable stick-slip sliding model that takes into account the dynamic response of the system.However, this model uses an equivalent linear approximation that does not capture the change in dynamic soil properties with time, and it models permanent shear strains due to sloping ground through a coefficient of friction.
The main application of the multidirectional models is in conjunction with a 2D finite element or finite difference numerical simulation where only one component of ground motion is applied in the downslope direction.However, the multidirectional models can also be combined with results from simplified displacement models (e.g.Jibson [40], Saygili and Rathje [41]; Bray et al. [3], Bray and Macedo [4]) to provide a quick initial estimate of seismically induced permanent displacements for a single slope, or for regional mapping of landslide hazard for a given earthquake scenario [42].The multidirectional models could also be combined with a displacement model to conduct performance based probabilistic slope stability analyses [43,44], or combined with USGS ShakeMaps to generate SlideMaps [45], which provide near real time displacements based on ground motion conditions during the time of an earthquake event.
For uniform slopes with fixed ends in the lateral direction, one could use the model of Azizian and Popescu [24] to estimate the decrease in displacements due to 3D geometry effects and equation ( 1) to estimate the increase in displacements due to multidirectional shaking effects.For these types of slopes, the two phenomena combined predict u 2 = (0.8-1.40) * u 1 for 1 ≤ W/H ≤ 2 and 0.5 cm ≤ u 1 ≤ 20 cm.For W/H ≥ 2.5 and u 1 ≤ 100 cm, the effect of multidirectional shaking outweighs 3D geometry effects and u 2 is always larger than u 1 .This shows the importance of multidirectional shaking and that it cannot always be assumed to be cancelled out by a reduction in displacements due to 3D geometry effects.
While the ML models generally show better performance than the conventional least squares regression model (as measured by the correlation coefficient R 2 ), they are more difficult to implement in practice.As a result of this limitation, combined with the fact that the performance is only moderately better than the linear regression model (R 2 = 0.98 vs R 2 = 0.91), we recommend using the least squares regression model displayed in equation (1) for practical engineering applications.

Conclusions
This paper developed a model to estimate seismically induced permanent displacements due to multidirectional shaking on a slope using the results of a traditional unidirectional seismic slope stability analysis as input.To develop the model, we first created a simulated database of  slope displacements due to both multidirectional and unidirectional shaking.The database consists of 44,160 three-dimensional finite element slope stability analyses using 230 ground motion record pairs, 48 soil profiles, and four different combinations of ground motion components (see Table 1).The four combinations of ground motion components are (1) component A applied in the x-direction (downslope); (2) component B applied in the x-direction; (3) component A in the x-direction (downslope) and B in the y-direction (across slope); and (4) component B in the x-direction and A in the y-direction.To prevent analyses with negligible displacements from dominating the regression analysis, we removed all analyses with maximum permanent displacements less than 0.5 cm and the corresponding other three analysis cases, which resulted in a final database of 28,100 analyses.
The ratio of maximum permanent displacements calculated from applying two components to only applying one component (MDF) is always greater than one, which means that multidirectional shaking produces larger permanent displacements than unidirectional shaking.The mean value of MDF is 1.76 and the maximum value is about 11 for the ground motions and sites investigated in this study.However, MDF decreases as the expected displacements from unidirectional shaking increase.
We performed least squares regression analyses as well as applied machine learning techniques to the simulated database of displacements.The model from least squares regression (equation ( 1)) is only dependent on the predicted displacement from unidirectional slope stability analysis (u 1 ).It provides an efficient and straight-forward way to estimate displacement due to multidirectional shaking from analyses using only one component of a ground motion.The standard deviation (equation ( 2)) is also dependent on u 1 , and decreases as u 1 increases.The standard deviations are similar or smaller than the standard deviations for other ground motion models.We investigated six different machine learning algorithms and found that the Multi-Layer Perceptron (MLP) model gave the best fit to the data when using only u 1 as input and the Extreme Gradient Boosting Algorithm (XGBoost) when using all of the input parameters in Table 2.We have included the tuned MLP and XGboost models as an electronic supplement, however, we recommend using equation (1) for practical engineering applications due to its simplicity and only slightly lower predictive accuracy.
The intended use of the models presented in this paper are to give a quick, first order approximation of the increase in seismically induced permanent displacements due to multidirectional shaking based on the results of a traditional unidirectional seismic slope stability analysis.They can also be used with simplified slope displacement models (e.g.Jibson [40], Saygili and Rathje [41]; Bray et al. [3], Bray and Macedo [4]) in regional landslide hazard assessments [42] or performance based probabilistic slope stability analyses [43,44].However, the models are not intended to replace more detailed analyses for critical structures and should not be extrapolated outside the range of parameters used to develop them or for liquefiable soils or soils with strain softening behaviour that could lose significant strength during or after earthquake shaking.
Finally, the results of this study showed that the increase in displacements due to multidirectional shaking cannot be assumed to be cancelled out by a reduction in displacements due to 3D geometry effects, and both should be considered in dynamic slope stability analyses for critical slopes.

Fig. 2 .
Fig. 2. a) Dynamic shear strength and b) shear wave velocity with depth for soil profiles used in the analyses.

Fig. 3 .
Fig. 3. (a) moment magnitude (M w ) vs. the Joyner-Boore source to site distance (R jb ); and (b) PGA of component A (PGA A ) vs. PGA of component B (PGA B ) for the selected 230 ground motions.

Fig. 4 .
Fig. 4. Comparison of slope displacements with time estimated from PLAXIS 2D and OpenSees for slope heights of 60 m and slope angles of a) 5 degrees, b) 10 degrees and c) 15 degrees, and d) slope height of 100 m and slope angle of 10 degrees.

Fig. 5 .
Fig. 5. Histograms and mean (μ), median (x), and standard deviations in natural log units (σ ln ) of a) maximum permanent displacements due to all analyses with only a single component (u 1 ); b) maximum permanent displacements for all analyses with two components (u 2 ); and c) multidirectional factor (MDF = u 2 /u 1 ).
the model with the highest coefficient of determination (R 2 = 0.91) and lowest standard deviation in natural log units (σ ln ), where u 2 and u 1 are in cm.The orange line in Fig.6shows equation(1).ln(u 2 ) = 0.923 * ln(u 1 ) + 0.663

Fig. 8 .
Fig. 8. a) Predicted values of u 2 from equation (1) versus the values calculated from the numerical simulations and residuals (in natural log units) versus u 1 .The grey line is the 1:1 line representing equality.

Fig. 10 .
Fig. 10.Premutation feature importance for input variables used in ML regression, determined from the XG Boost model.

Fig. 11 .
Fig. 11.Performance metrics for different ML regression models for predicting the displacement u 2 when using u 1 as the only input parameter.The blue dots are individual results from cross-validation, and the orange lines and error bars represent the mean and standard deviation for each of the model results.

Fig. 12 .
Fig. 12. a) Predicted values of u 2 from the XG Boosting regression model using the full set of input parameters in Table 2 vs. the u 2 values calculated from the numerical simulations, and b) residuals (in natural log units) versus u 1 .The grey line is the 1:1 line representing equality.

Table 1
Overview of slope stability analyses.