Calibration of Load and Resistance Factors for Breakwater Foundations under the Earthquake Loading

: This study investigates the system stability of breakwater foundations subjected to earthquakes from a probabilistic point of view. A fully probabilistic approach, i.e., a combination of the Monte Carlo simulation and Bishop’s simpliﬁed method, has been developed to evaluate the system failure probability of foundation damage, one of the prevailing failures encountered during earthquakes. Twelve sections of perforated caisson breakwaters located around Korea were chosen as case studies. First, the reliability analysis was performed for all the breakwaters at existing conditions; then, the calibration process involving the estimation of load and resistance factors was conducted for 12 breakwaters at three levels of the target reliability index. As the performance function, used in the stability analysis of breakwater foundations, is deﬁned based on an implicit shape with a high-dimensional space of variables, the calibration process of load and resistance factors becomes cumbersome and complicated. Therefore, this study has proposed a sensitivity analysis to be implemented prior to the calibration process to elicit the effects of variables on the stability of each breakwater, which, thereafter, effectively directs the calibration process. The results of this study indicate that the failures in the foundation of breakwaters frequently occur in different modes. Therefore, the failure probability should be estimated considering all possible failure modes of the foundation. The sensitivity results elucidate that the soil strength parameters are the dominant variables, contributing to the stability of foundations, whereas the seismic coefﬁcient presents the negative effect, causing the insecurity of breakwaters. In particular, the deadweights, though directly contributing to the seismic forces, show a small effect on the stability of foundations. The calibration shows that the load factors slightly vary with an increase in the target reliability index and set 1.10 for three safety levels. In contrast, the resistance factor exhibits an inverse relationship with the speciﬁed reliability index. Especially when the load factor equals 1.10, the resistance factors are 0.90, 0.85, and 0.80, corresponding to the reliability index of 2.0, 2.5, and 3.0, respectively. Eventually, it is proved that the sensitivity analysis prior to the calibration process makes the procedure more Accordingly, the iteration of simulation execution is diminished, and the convergence is quickly accomplished.


Introduction
Reliability-based design concepts and their applications to load and resistance factor design (LRFD) have proven to be efficient for cost-saving and material use [1,2]. In general, there are three levels of reliability-based design. The level III method, the fully probabilistic approach e.g., Monte Carlo simulation (MCS), derives the failure probability based on the achieved. In other words, the RI of the designed structures will be larger than the RI that forth stipulated in the calibration process of LRFs [11].
Currently, factors in Equation (1) are calibrated using a general closed-form solution (CFS) of an approximated linear equation comprising load and resistance terms. The load terms are commonly assumed as the normal distribution, and the resistance component is assumed as the lognormal type [12]. Statistical insights of load and resistance are represented throughout their bias factors, which account for the difference between the actually measured and the predicted (nominal) values. The mentioned LRFs were calibrated using the level II probabilistic approach, e.g., MVFOSM or FORM [5,12,13]. The calibrated LRFs using CFS are dependent on the approximation process adopted for estimating statistical properties of the load and resistance bias factors [12]. This approach is applicable for the cases where the measured loads and resistances are available as structural components [5], spread footings [14], soil nail walls [11,12,15], or driving piles [1,16]. However, in the case of the bearing capacity of the BRW's foundations, the measured bearing capacity in terms of circular slip is rarely available and, sometimes, impossible to implement [17]. In this context, the MCS become the applicable choice to derive reliable LRFs.
With regard to port facilities, there are two types of environmental loading: wave loadings and earthquakes. The damages of BRWs under the environmental loading are classified as the overall and local failure modes [18]. The sliding and overturning damages are representative for the overall failures of the superstructure. Whereas, settlements and shear failure due to the insufficient bearing capacity are common overall failure shapes occurred in the foundation. The dynamic stability of BRWs during earthquakes and wave impacts was further investigated [19][20][21][22]. The sliding of caisson structures along the rubble mound surface caused by breaking wave loads was investigated using Newmark's method [20]. In their study, caissons were modeled as a mass-spring dash-pot linked to the ground. The effect of parameters, such as the mass, stiffness, and damping coefficient, on the dynamic response, as well as horizontal and rotational displacement of caissons were investigated and validated based on the large-scale model test. In addition, the dynamic response of the upright BRW subjected to impulsive breaking wave forces was studied and verified with the small-scale test of the pendulum [21]. Moreover, the dynamic response of caisson structures, located in the array, was studied [19]. Their study examined the response of caissons considering the effects of adjacent caisson blocks under the wave impact loading. Furthermore, the sliding and overturning stabilities of caisson blocks under the combined effect of earthquake and tsunami were investigated [22].
Furthermore, the historical failure observations revealed that BRWs were mainly damaged owing to their foundation failures when they were found to be vulnerable during earthquakes [23]. However, behaviors of a BRW under the wave effect were more attractive to study than those relating to the seismic excitation. In addition, previous studies mainly concentrated on the stability of caisson blocks, i.e., sliding or overturning failure modes within the CDAs. In contrast, the foundation stability of BRWs subjected to earthquake has been barely examined.
To evaluate the stability of the slopes, several methods can be applied, e.g., the strength reduction method, and limit equivalent method (LEM). Recently, the method of shear band propagation has proved to be the most appropriate approach to evaluate the safety factor of the cut slopes that consist of weak clay layers [24][25][26]. Oumeraci et al. [7], developed several equations corresponding to the predefined failure surfaces in the BRW foundation using the upper bound theory. These equations with simplifying hypotheses were expressed in the explicit definitions implying the convenient estimation of the derivatives of the performance function in the reliability analysis. However, the predefined failure surfaces might not be representative for all possible failure surfaces that occurred in the foundation.
Among the LEMs, Bishop's simplified method (BSM) is recommended as the suitable approach for verifying the stability of the BRW foundation [8,17].
Based on the aforementioned concepts, this study aims to investigate the foundation stability of BRWs subjected to seismic loading using the fully probabilistic approach. Safety in terms of RI or failure probability (PF) was investigated for 12 BRWs around Korea. Moreover, LRFs were calibrated for three levels of target RIs (2.0; 2.5, and 3.0) using the MCS. The fully probabilistic procedure proposed in our previous study [27], which comprises BSM and the MCS, evaluated the probability of failures in the reliability analysis. Based on the fully probabilistic approach, the calibrated LRFs are expected to provide an efficient solution in the design process of new BRWs. The rest of the paper is set as follows: First, the seismic analysis approaches are discussed to opt for the appropriate method. The stability of the BRW foundation should be considered as a system, including all possible failure modes as described in Section 3. In Section 4, the sensitivity analysis is carried out to discover the prominent variables showing dominant effects on the stability of each BRW. Section 5 illustrates the MCS for the PF evaluation of the BRW foundation-induced earthquakes. Section 6 surveys the target reliability index. Section 7 presents the calibration process of LRFs by incorporating a sensitivity analysis. Section 8 discusses the results of the study. Section 9 concludes the paper.

Seismic Analysis in Reliability Assessment
There are some cases in which the earthquake-induced stability verification is omitted during the performance assessment of BRWs. However, if BRWs are located in deep water and the design wave height is small, then the actions associated with seismic excitation can become the dominant factor [8]. For deep-water sites, perforated caisson BRWs, including wave absorbing chambers, can considerably mitigate wave impacts, becoming a reasonable design choice. However, owing to the presence of wave chambers, the structures become more eccentric with regard to its deadweight, which result in a more unstable state compared to the upright caisson sections [28]. Therefore, the investigation of the overall stability in terms of foundation stability supporting perforated caisson sections has been conducted in this study. Figure 1 illustrates the typical section of the perforated caisson. Their overall design conditions are listed in Table 1. To analyze the structure behavior-induced earthquakes, three methods of earthquake engineering, which are widely accepted and applied to the port facility structures, are the pseudo-static method, Newmark's method, and dynamic analysis. In the pseudo-static approach, or seismic coefficient method, the cyclic nature of the earthquake is ignored, and seismic forces are treated statically, i.e., the equivalent static components of the seis-mic force are applied to the centroid of the sliding mass that acts along the out-of-slope direction. The seismic coefficient is converted based on the acceleration time history of an earthquake, considering the seismic response of the ground. Owing to the transient nature of earthquakes, the effective seismic coefficient, used in the pseudo-static approach, should be a fraction of the ratio between the peak ground and gravity accelerations rather than the full ratio. This fraction depends on the peak ground acceleration, and the average ratio of 0.6 was proposed [29,30]. Based on Newmark's method, known as simplified dynamic analysis, permanent deformations of BRWs are evaluated for the cases where pseudo-static safety factors during an earthquake are less than unity. The acceleration at which the failure occurs is defined as the threshold seismic value. During an earthquake, the deformations of the sliding mass are accumulated, at which the acceleration is larger than the threshold value. In general, the dynamic analysis uses finite element or finite difference approaches in which the local site effects are often treated as part of the soil-structure interaction analysis. First, the soil foundation is considered as equivalent linear, non-linear total stress, or effective stress model. Then, the seismic performances are evaluated. Fairly comprehensive results can be obtained from dynamic analysis, such as the failure modes of soil-structure systems, extent of deformation, stress-strain states in soil, and structural components. However, seismic responses based on Newmark's method and dynamic analysis are sensitive to the input accelerations. Moreover, BRW behaviors depend on the characteristics of earthquakes such as intensity, magnitude, duration of excitation, and frequency, as well as the features of structures, i.e., the stiffness or self-weight of structures.
Historical observation reveals that the highest magnitude of the earthquake that occurred in Korea is about 7.5. In the Korean standard [31], there are two basic specified seismic zone factors, i.e., 0.11 and 0.07 with 500 years of return period to account for the different local sites. Based on the specified seismic zone factors, the seismic design spectrum will be determined with respect to the ground soil type, and the importance of the designed structures. When earthquake is considered, seismic forces can cause dynamic effects in three directions: two horizontal and one vertical. The BRW structures are known to be comparatively rigid structures, leading to small vibrations during earthquakes compared with seismic excitations. In addition, their fundamental period is short, and the damping factor is large. Therefore, except for BRWs located near the seismic center, the influence of the vertical seismic component is not as high as that of the horizontal component [30]. Hence, in this study, the seismic responses owing to horizontal excitation are further investigated using the seismic coefficient method within the time-consuming approach of the MCS.

System Reliability Analysis
As described in the previous studies, the location of critical failure slip surfaces in slopes depends on the soil stratigraphy and soil strength properties [32,33]. Cho [33] analyzed the stability of two soil slopes consisting of two clay layers using multipoint firstorder approximations. His study demonstrated that although there is one deterministic failure surface, the results of the probabilistic analysis present two different failure surfaces that can be encountered: one crossing path at the upper layer and the other at the upper and lower layers. Moreover, several soil slopes were examined by Zhang et al. [32]. In their study, the most critical failure slip surfaces were located as the representative failure modes based on the reduced shear strength parameters. Thereafter, the safety factor (FS) of each failure surface was depicted by a second-order polynomials using the response surface approach combined with Hassan and Wolf algorithm [34]. Then, the system PF was evaluated using the MCS based on the fitted second-order polynomials.
In this study, the system reliability analysis is investigated using the MCS rather than the approximated multi-point FORM or response surface approach where the achieved PFs are mainly dependent on the identified major failure modes. The results of the MCS show that there are some cases where the failure surface is represented by one slip surface. Alternatively, some problems show several representative slip surfaces. The system PF is, hence, calculated as the ratio of the total number of failures of all representative slip surfaces to the total size of the simulation. Figure 2 presents the results of the MCS for Port Nine. The foundation of Port Nine consists of four soil layers where layer one is 7 m of the artificial rubble mound arranged above three sediment sand layers. The thicknesses of the three soil layers are 2 m, 8 m and 11.6 m, and the friction angles of the layers are 32 • , 33 • , 34 • , respectively. The shear strength properties of the rubble mound were comprehensively investigated in the study of Kobayashi et al. [17]. By following that, OCDI 2009 [8] specified the shear strength parameters as 20 kN/m 2 for the apparent cohesion and 35 • for the internal friction angle as the standard values of rubble mounds. The failure modes are affected by the soil properties of the soil foundation. Generally, the rupture surfaces intend to occur in the soil having relatively low shear strength [24,33]. In the foundation of Port Nine, the shear strengths of sand layers are quite similar, hence, the failure modes would potentially cross path any layers. The dashed red line in Figure 2a illustrates the critical failure surface using CDA, whereas the blue and black lines depict the possible failure modes resulted through the MCS. There are two potential failure modes in the MCS. The first failure surface, the black line, crosses four upper layers, and the second failure mode, the blue line, occurs in two upper layers. The two bundles of driving loads and resistances (shown as scattered in Figure 2b) are associated with two failure modes. It should be noted that if one considers only the most representative failure mode, i.e., mode one (gray circles in Figure 2b), then the PF would be merely determined by the number of failure points that belong to the mode one realization. It implies that the failure cases belonging to mode two were neglected. Therefore, the result of PF, based on the representative failure mode, would be smaller than the actual system PF. Although the PF of a representative failure can be a good approximation for a homogenous slope, the system PF can be significantly larger than that of the most critical failure mode in the foundation including layered soils.

Sensitivity Analysis
To make the calibration process effective, the effect of random variables on the stability of the structure should be surveyed. The sensitivity of the variables exhibits the extent of their influence on the safety level. Several methods can be applied in the structure sensitivity analysis [35][36][37]. Especially, the perturbation technique was applied to derive the boundary value problems and the system of equations based on which of the sensitivities, i.e., the partial derivatives of the objective function with respect to all the parameters, can be analytically obtained [35].
Alternatively, for a problem including a large number of variables, the tornado diagram method is commonly used. The tornado diagram depicts the relative sensitivities of input variables on performance functions as "swings," where the larger the swings are, the effect of the variables will be more significant on the performance function. The swing of each variable is constructed based on its upper bound and lower bounds. For slope problems, the safety factor, defined as the ratio of the shear strength to the stress along the slip surface, reflects the stability of the slope under driving forces, such as surcharge load, self-weight, and seismic load. Few well-known approaches proposed by Bishop, Spencer, and Morgenstern-Price based on the LEM have been applied for slope stability analysis [38]. However, those methods define the safety factor in an implicit manner. Therefore, it is difficult to directly calculate the analytical derivatives of the objective function with respect to variables.
In 1999, Hassan and Wolf proposed the algorithm (HWA) to estimate the RI of slopes based on the MVFOSM [34]. In HWA, the partial derivative of the safety factor with respect to each variable is numerically approximated by determination of FS using parameter values that are greater or less than its mean by several times (m times) of its standard deviation. The sensitivity of variable X i is determined using Equation (2). Derivatives estimated using Equation (2) reflect the slope at the mean point of the performance function with respect to each variable in the hyperspace of input variables.

∂FS ∂X
where FS i + and FS i − correspond to safety factors estimated considering the upper and lower values of the variable X i , while the other variables are maintained at their means, µ. The increment times, m, larger than or equals to unity are concluded to have a minor influence on relative sensitivities; hence, it is considered as 1.0 in this study. The safety factors in Equation (2) are determined using BSM as shown in Equation (3) in Section 5.
The sensitivity analysis results are plotted in Figure 3 for Port One. The foundation of Port One includes the 11 m thickness of rubble mound sitting on the sand layers (with 6.1 m of thickness) above another sandy gravel layer. The effective friction angles of sand and sandy gravel layers are 30 • and 40 • , respectively. In Figure 3a, the dashed black line, or the baseline, shows the mean values of safety factors that are calculated using the mean of all the input variables. The input variables are presented later in Section 5. The blue and red lines denote the safety factors estimated at the upper and lower bounds of each variable, respectively. It is perceived that when the lower bound of the soil strength parameters or the upper bound of the seismic forces is applied, the safety factor will be reduced and vice versa. The variation of safety factors in the sensitivity analysis depicts the extent of the sensitivity of the variable as demonstrated in Figure 3b. The red square markers in Figure 3b illustrate the variables having adverse effects, whereas the blue diamond markers demonstrate those having positive effects on the stability. The larger the sensitivity of the variable, the more would be the influence of the variable on the foundation stability. Furthermore, the weight components exhibit small adverse sensitivities, whereas the seismic coefficient presents a considerable influence on the insecurity of the BRW foundation, although they both contribute to the seismic force.

MCS Approach
The MCS is a robust and consistent method to simulate the problems by considering uncertainty variables to evaluate the PF. The MCS is the best solution for problems in which the objective function is modeled based on a large dimension, exhibits nonlinear behavior or is defined based on the implicit aspect. Uncertainties are introduced into the MCS to reflect the inherent randomness of a natural process. In the case of soils, even in the nominally homogeneous soil layers, soil properties may exhibit considerable variation from one point to another. This phenomenon is inherent and known as the spatial variability of soil. In reality, stratigraphy and soil parameters are mainly dependent on the overall process comprising four stages: weathering, erosion, transport, and deposition. Therefore, uncertainties reflect the lack of insights owing to the lack of necessary information. In this study, the geometry of the sections of BRWs is assumed to be similar to the designed manner. The concerned uncertainties are related to soil properties; tide variation; seismic excitation; and construction work, i.e., the unit weights of materials and armor layers. The statistical parameters of the input variables are listed in Table 2. The shear strength parameters of soils are more variable compared to the weights. The coefficient of variation (COV) of 6% is specified for the shear strength properties in OCDI 2009 [8] to develop the PSFs. In principle, the COV of variables reflects the degree of the site understanding, and the low variability of the shear strength of soils is reported in the range of 10% to 30% [39]. In this study, the COV of 10% is proposed in the reliability analysis. The filling material are mostly dependent on the construction work at the field, hence, its weight presents the most variation among the considered participation weights. BSM is a widely accepted method based on LEM, in which the slip surface is assumed to be circular, and the safety factor is determined with the equilibrium of moments around the slip center [8,17]. The procedure to combine the MCS and BSM to estimate the safety factor, as well as RIs, was successfully investigated in the previous study [27]. In this study, the same procedure is applied considering the horizontal seismic force in the denominator of Equation (3) to estimate the safety factor, FS. The horizontal seismic force is the product of the seismic coefficient and net deadweight without deducting buoyancy. Then, the performance function, g is defined in Equation (4).
where w' and w are the effective weight and total weight of the comprising slice segment, respectively. q Ex is the surcharge load distributing on the slice along the vertical direction. F Eq and a correspond to the horizontal seismic force and its lever arm about the failure center. R is the radius of the considered failure surface. c is the apparent cohesion of the soil (or undrained shear strength in the case of cohesive soil), and ϕ is the friction angle of the soil along the base surface of each slice. s is the width of each slice. Finally, θ is the inclination angle, horizontal to the base surface. The horizontal seismic force increases the stress along the failure slip; hence, the safety factor is decreased. In Equation (4), vector X includes the considered constant and uncertainty variables. For each sampling, the performance function g is determined, and the number of failure cases is accumulated (for the case that g is less than zero) as N fail . Then, the probability of failure, P f , is determined as the ratio of total failure cases to the size of simulation, N (Equation (5)). The associated RI and probability of failure are expressed in Equation (6).
The result of P f in the MCS is known to be sensitive to the simulation size; therefore, the number of realizations should be surveyed before conducting a simulation. In general, the size of simulation should be greater than ten times the reciprocal of P f (i.e., N > 10/P f ), so that the coefficient of the variation of P f (COV Pf ), estimated based on the different MCS runs, is less than 0.3 [27]. The relative error, ε, with a 95% confidence interval of the binomial distribution of each simulation cycle is evaluated using Equation (8) [40].
For example, if the maximum target RI is set as 3.0 (or equivalent P f of 0.0013), then the minimum size of simulation, N, estimated by Equation (7), should be 2500 to reach a COV lower than 0.3. In this study, PFs are evaluated using the simulation of 200,000 samplings. Considering the simulation size of 200,000, the maximum error with 95% confidence interval is expected to not larger than 0.12%. In this study, this size in the MCS has been used to evaluate the RI of the existing structures and to calibrate LRFs with respect to the specified target RIs. Within the MCS of the existing structures, if there is no failure event encountered, the PF is estimated based on the approximated distribution function of FS that is fitted based on its MCS outcome.

Target Reliability Index
It is given that LRFD, if applied, would provide a uniform and consistent safety level. As a result, LRFD is assured to generate a safe and economical class of structures. However, prior to the calibration of LRFs, the target RI, β T , should be appropriately specified. In principle, the specified safety levels may directly affect the design process, i.e., the type of structure, applied material, and cost of investment. Therefore, the specified safety level should be linked to the potential consequences of failure. In general, for more important structure, the higher safety level should be specified. The safety level used in LRFD codes should be at least specified to yield those applied in earlier codes. Pertaining to AASHTO LRFD Design Specifications, the LRFs for strength I limit state have been calibrated for a target RI of 3.5 during a 75-year design life of a bridge [10,15]. However, the stability design of bridge substructures (i.e., slopes and foundations) is categorized into the service I limit state, leading to a load factor of 1.00. The associated resistance factors are calculated by fitting into FS in CDA codes.
The range of the target PF for slopes was set between 0.001 and 0.01 (related to RIs of 2.3-3.1). A value of 0.01 is for low potential risk of designed slopes, whereas 0.001 is for high potential risk [41]. The sliding failure mode of caisson blocks along the rubble mound surface was investigated [42]. The target RI was obtained using the statistical information of the safety factors based on the deterministic approach. The target RIs, close to the specified values of 2.367 and 2.433 in OCDI 2009 [8] corresponding to the consideration of the normal or lognormal distribution of the target RI, were derived. Besides, several internal failure modes of BRWs where the deep mixing method was applied to improve foundations were investigated [43]. The value of 3.1 was specified as the system target RI with respect to internal failure modes. Moreover, sliding and overturning failures of gravity quay walls were examined for 30 different design conditions [9]. The target RI was set as the minimum RI among 30 considered quay walls. To investigate the trend, although the minimum RI is 2.32 for Port Five in this study, three levels of target RIs are chosen as 2.0, 2.5, and 3.0 to calibrate LRFs.

MCS for the Calibration Process
Previously, LRFs were mainly calibrated using FORM based on a general closed-form solution [1,11,16]. However, FORM cannot be directly implemented in this study because highly non-linear implicit function evaluations are needed in the process [9]. Moreover, there is no analytical function that can exactly reflect all the failure modes of a system to adopt FORM. Alternatively, MCSs were successfully applied in the calibration process of load and resistance factors in the slope and retaining wall problems [44,45]. In the calibration process of LRFs using the MCS, system characteristics are adjusted so that the safety level of the system meets the given target RI. In other words, the nominal values of the variables of the system are iteratively changed in their spaces, and then, the PF is evaluated and compared with the target PF to verify the convergence [9,44]. In principle, several combinations of input variables produce the target RI, especially for the case of a high-dimensional and non-linear problem, such as the BRW foundation stability that includes uncertainties in the multilayered subsoil and superstructure components.
To effectively perform the MCS for the calibration of LRFs, the sensitivity analysis has been proposed to analyze the insight of variables on the performance function, as described in Section 4. Based on the sensitivity analysis, the input variables that have minor effects on the safety factor are retained for the initial condition. Meanwhile, variations have been iteratively applied to variables that have a dominant effect on the objective function. Then, the PF is evaluated for each searching iteration. The searching converges until the difference of the PF of the adjusted system and the target PF falls within 5% [44]. This procedure neglects variations applied to the variables that have a minor effect, and results in a more efficient calibration process thanks to the insight about variables affecting the performance function.
The calibration process for Port 12 at a target reliability index of 2.0 is depicted in Figure 4. The load and resistance components are plotted as gray circles for the initial condition. It can be seen that within the 200,000 sampling of the MCS no failure event was found in the figure. It means the resistance values are always greater than those of load in the MCS at the initial condition. In this situation, the PF is approximated using the analytical probability density function of the safety factor, which is fitted to the results of the MCS. During the calibration process, the nominal values of the dominant input variable (the horizontal seismic coefficient in this case) has been increased, and the others are retained. Load and resistance values resulting from the variation of dominant variables for the target RI of 2.0 are scattered, as shown by blue dots in Figure 4. The realizations in the calibration process are mainly transmitted toward the vertical direction, which implies that the loads are considerably enlarged, whereas resistances are slightly increased. Both loads and resistances correlated to the weight of the sliding mass have been increased, implying that the critical slip surface is extended in the soil foundation when the seismic coefficient is increased. Furthermore, the statistical characteristics, i.e., the mean and variance of the load, as well as the resistance, are varied during the calibration process. This phenomenon reflects the non-linear behavior of the performance functions represented by the implicit function [27]. The number of failure cases, i.e., the number of points in the failure zone, is 4496 considering the variations; the system PF become 0.0225 (i.e., 4496/200,000). Thus, PF falls within a 5% difference compared with the target PF of 0.02275 (corresponding to the RI of 2.0), implying that the searching process is converged. The limit state points (points approximately lying on the limit margin (R = Q), presented as brown points in Figure 4) are singled out to evaluate LRFs using Equation (9) [9,44].
where Q and R correspond to loads and resistances. LF and RF are the calibrated load and resistance factors, respectively. The subscript n is abbreviated for nominal components, whereas LS stands for the limit state points. The statistical features of limit state points (similar to the most probable point or the design point in FORM) are discussed in Section 8. Then, the calibration process is implemented for all ports at three levels of target RIs, and the calibrated LRFs are presented in the next section.

Results and Discussion
In Figure 5, the initial safety factors based on CDA and corresponding RIs resulted from the MCS are presented. The smallest CDA safety factor is 1.274 associated with the lowest RI of 2.32 in the case of Port Five. The existing structures were designed to meet the requirement of the minimum stipulated safety factor of 1.0, in the current standard [31], for earthquake verification. Hence, the minimum safety factor of all the existing BRWs is not lower than unity. The linear regression line is approximated at a high coefficient of determination (R 2 = 0.863). Based on the regression line, the associated RI corresponding to the CDA safety factor of the existing BRW can be approximated. For example, if the safety factor based on CDA is 1.0, the associated RI can be estimated as 1.05. The COV at the initial condition of safety factors and RIs are 0.17 and 0.36, respectively. This implies that the safety level of the existing structures presents a large variation in the probabilistic point of view.  Figure 6 shows the distributions of safety factors for Ports Five and Six. Although the CDA safety factors do not differ much (~1.15 times), the PF (PF when safety factor is lower than unity, the area of the histogram on the left side of FS = 1 margin) for Port Five are 340 times larger than that for Port Six. This implies that pertaining to the probabilistic point of view, the results of the safety factors using CDA do not reflect a uniform safety level, which should be an accomplished objective to reach an economical and consistent design of structures. Hence, the traditional CDA should be progressively converted to probabilistic approaches, e.g., LRFD, to achieve a consistent and uniform safety level. Figure 7a depicts the distribution of load and resistance at limit states (LSs) where the safety factors closely equal unity (brown points in Figure 4). The bottom part of Figure 7a illustrates that the safety factor at the limit state exhibits a uniform distribution. Moreover, the major overlap at the top part of Figure 7a exhibits the most similar distributions of load and resistance at the limit state. This proves that the chosen range of safety factor as the limit state (FS within 0.99 and 1.01) is sufficient enough to approximate limit points. In addition, the histograms of loads and resistances lying on the limit state margin are identical to the normal distribution; hence, the estimated LRFs using the means of load and resistance denote the most probable points. Figure 7b depicts the histograms of loads and resistances pertaining to the MCS. The limit state points of load and resistance represented by their mean values are identically placed at the same value denoted by the solid black line. The limit state points are located between the nominal load and resistance. Therefore, the load factor will be greater than unity, and the resistance factors will be less than unity (Equation (9)). The calibrated LRFs for 12 ports at three safety levels are listed in Table 3. Based on the small variation of the achieved RIs (~0.2% for three levels) during the calibration process, the variation of calibrated LRFs is less than 7%, implying a good convergence during the calibration process.  In Figure 8, the resistance factors, normalized with associated load factors, are plotted with nominal safety factors for all ports at all three target safety levels. The nominal safety factors, shown in Figure 8, are safety factors evaluated with respect to obtained nominal values that meet the requirement of the specified target RIs. The normalized resistance factors are linearly inversely related to the nominal safety factor. The high regression coefficient of 0.982 implies a strong correlation between the nominal safety factors and normalized resistance factors after the calibration process. Based on the regression line, the equivalent normalized resistance factors associated with the specified safety factors can be derived. For example, if the existing safety factor is specified as 1.0, then the associated normalized resistance factor should be 0.96. During the calibration process, the nominal FS varied as the input variables were adjusted to meet the specified target RIs. For instance, the RI at the initial condition of Port Five is 2.32, which is smaller than the target RI of 3.0. Therefore, the nominal safety factor is increased from 1.274 at the initial condition to 1.373 to reach the target RI of 3.0. Alternatively, the initial RI estimated as 4.05 for Port Six is larger than the target RI of 3.0; hence, the safety factor should be reduced from 1.459 to 1.343 to reach the target RI of 3.0.  Safety factors based on CDA at the initial condition are sorted and plotted as a dashed black line in Figure 9. In addition, the resistance factors, normalized with the associated load factors, are illustrated for three safety levels. If a higher safety level is specified, a lower resistance factor should be used. The safety factors show a larger variation (COV of 0.17) compared with the results of the resistance factors (highest COV of 0.04). Moreover, although safety factors, which are sorted, present the increasing trend, the corresponding normalized resistance factors do not exhibit a clear trend. It implies that the resistance factors slightly rely on the characteristics of the chosen initial structures. In other words, the chosen structures comprising different nominal parameters do not influence the calibrated resistance factors. The reason for this can be demonstrated using Equation (9), where the load and resistance forces at the limit state are normalized with their nominal values, making the calibrated LRFs slightly dependent on the initial nominal safety factors. This phenomenon illustrates the same behavior with LRFs for bridge design, in which the resistance factors are not dependent on the magnitude of live load and dead load components, but the ratio of the components [5,13]. To propose effective factors for the design of new structures, LRFs should be validated. First, LRFs were rounded to the nearest value of 0.05 [46]. Then, the calculated resistance factor was verified for the initial structures, and the safety condition was again evaluated. The proposed LRFs considers the average of 12 ports and are listed in Table 4. Figure 10a,b illustrate the verification results for Port Five using proposed LRFs at target reliability values of 2.0 and 2.5, respectively. The red color depicts the load effects; the blue color represents the resistance components. Nominal results are derived from the stability analysis with regard to the nominal (design) input values. Factored components are products of terms with their factors, and the limit states were verified using Equation (1). As shown in Figure 10a, the limit state equation is satisfied for the target RI of 2.0. However, when the proposed LRFs at the safety level of 2.5 are applied, the limit state equation is not satisfied as demonstrated in Figure 10b. This implies that the failure probability of foundations may be greater than 0.0062 and smaller than 0.023. This situation is figured because the initial RI of Port Five lies between 2.0 and 2.5.

Conclusions
This study considered the earthquake-induced stability of BRW foundations using a fully probabilistic approach, i.e., a comprehensive approach combining BSM and the MCS. The stability of BRW foundations with regard to the PF or reliability index was examined as a system safety assessment. The sensitivity of variables was analyzed using HWA, where the sensitivities were numerically evaluated based on the relative variations of safety factors corresponding to the increment or reduction of one standard deviation from the mean of each variable. The sensitivity analysis illustrates that the shear strength properties of the soil foundation are the most important factors keeping the structure stable, whereas the seismic coefficient is the most adverse component causing the failure of foundation.
The minimum system reliability index of 2.32 was evaluated for Port Five. Based on that, three reliability indexes of 2.0, 2.5, and 3.0 were recommended, and the calibration process was performed to obtain the associated load and resistance factors using the MCS. The calibrated LRFs, validated with the CDA, were proposed to design new BRWs using the LRFD approach. The results of this study implied that LRFD provides a consistent and uniform safety level, demonstrated by small COVs of RFs compared with those belonging to the safety factors of CDA. Therefore, the LRFD should be progressively applied rather than CDA while designing new BRWs. It is expected that the proposed LRFs would be applied in the design process of new BRWs.
Furthermore, the failure modes in BRW foundations can be altered with the variation in soil properties or external excitations. Therefore, the system PF should be investigated using the MCS to properly examine all the possible failure modes instead of FORM or general closed-form solutions. Finally, with the implementation of sensitivity analysis prior to the calibration process to elicit the dominant input variables, the calibration process has become directional and efficient owing to the reduction of time for the analysis of the MCS.