ReliabilityBased Analysis of a Caisson Breakwater with the Application of Bayesian Inference
https://doi.org/10.1007/s11804021002378

Abstract
Caisson breakwaters are mainly constructed in deep waters to protect an area against waves. These breakwaters are conventionally designed based on the concept of the safety factor. However, the wave loads and resistance of structures have epistemic or aleatory uncertainties. Furthermore, sliding failure is one of the most important failure modes of caisson breakwaters. In most previous studies, for assessment purposes, uncertainties, such as wave and wave period variation, were ignored. Therefore, in this study, Bayesian reliability analysis is implemented to assess the failure probability of the sliding of Tombak port breakwater in the Persian Gulf. The mean and standard deviations were taken as random variables to consider dismissed uncertainties. For this purpose, the firstorder reliability method (FORM) and the first principal curvature correction in FORM are used to calculate the reliability index. The performances of these methods are verified by importance sampling through Monte Carlo simulation (MCS). In addition, the reliability index sensitivities of each random variable are calculated to evaluate the importance of different random variables while calculating the caisson sliding. The results show that the reliability index is most sensitive to the coefficients of friction, wave height, and caisson weight (or concrete density). The sensitivity of the failure probability of each of the random variables and their uncertainties are calculated by the derivative method. Finally, the Bayesian regression is implemented to predict the statistical properties of breakwater sliding with noninformative priors, which are compared to Goda's formulation, used in breakwater design standards. The analysis shows that the model posterior for the sliding of a caisson breakwater has a mean and standard deviation of 0.039 and 0.022, respectively. A normal quantile analysis and residual analysis are also performed to evaluate the correctness of the model responses.Article Highlights• Uncertainty of the number of waves and its influence on the maximum design wave height is modeled by random variables• A probabilistic procedure for sliding occurrence of caisson breakwaters is proposed for given geometry parameters.• The effect of different parameters on the sliding is approximated by the Bayesian linear regression.• Sensitivity analyses are carried out to assess the importance of each random variable• Reliability analysis is done to assess the sliding performance of a breakwater including the FORM and the Monte Carlo methods. 
1 Introduction
Breakwaters should resist random loads, such as waves and tsunamis. Therefore, their performance during such events is very crucial. Conventional design methods, which are also known as deterministic methods, cannot consider the real performance of such structures because of the unpredictable nature of random loads. In conventional design methods, structures are designed against specific return values of loads with safety factors. Here, safety factors are used to cover the innate unpredictable uncertainties of the loads and resistances of structures. Uncertainty can be epistemic and aleatory; the former is due to the lack of knowledge and can therefore be reduced by conducting more investigations or improving the models, whereas the latter cannot be reduced because of the randomness nature of this type of uncertainty. For example, the time of occurrence and maximum wave height over the next year or lifetime of the structure is unpredictable. Human and measurement errors and statistical uncertainties are epistemic and hence can be reduced (Faber, 2007). Uncertainty can be expressed through the concept of probability based on (1) external observations and repeats and (2) degree of belief in occurrences of an event. The first is the objective or classic approach, and the second is the subjective or Bayesian approach (Faber and Sørensen, 2002).
Caisson breakwater design is performed based on different criteria and environmental loads. Wave height and period are among the most important load factors that a breakwater encounters and must resist. Wind transfers its energy to a body of water, creating waves that move toward the coast. Wave height is affected by many phenomena before reaching and transferring its energy to breakwaters. The shoaling, breaking, diffraction, and refraction of waves due to water elevation and coastline morphology, tidal current, seabed bathymetry, and other unpredictable events can affect wave height. Taking these kinds of uncertainties into account is a major challenge in the reliability design (Goda and Takagi, 2000). The basis of the reliability design method is the randomness nature of loads and resistance of a structure or its components, including the epistemic and aleatory uncertainties and the decision on how to consider these randomness and uncertainties in the calculations. In most previous studies, the randomness nature of the mean and standard deviation of loads and resistances were ignored. In this study, these values are considered random variables to address some of the dismissed unpredictable nature of the mentioned values.
The reliability design procedure was developed in 1970 based on probability and statistics by Cornell (CORNELL CA, 1969). This method has been used since 1980 for breakwater design (Goda and Takagi, 2000). Toyama Toyama (1985), and Suzuki studied the consequences of the safety of breakwater sliding. Van der Meer applied another probabilistic procedure for the reliability analysis of rubblemound breakwaters (Meer (1988)). Wave height and return period were assumed to be known for the sake of assessments. This approach is one of the most significant shortcomings faced in previous studies. Goda et al. (2000) eliminated the safety factor before applying the reliability design method for breakwaters. They studied different wave heights with different return periods. Because sliding is the most probable type of failure in caisson and vertical breakwaters (Oumeraci, 1994; Oumeraci et al., 2001), in this study, the sliding failure mechanism was merely considered, and the randomness nature of wave heights is still not considered. Chaudhary et al. investigated the stability of a breakwater under the integrated effect of an earthquake and tsunami, and they inferred that sliding is the dominant failure mode (Chaudhary et al., 2017). Li et al. (2020) studied the Taichung harbor breakwater numerically and also reported sliding as the dominant failure mode. Loza and DeLeon used reliability analysis to assess the economic impact of port activities using the firstorder secondmoment (FOSM) approximation and Monte Carlo simulation (MCS). In their methodology, different failure probabilities can be obtained by the FOSM method because of the invariance nature of this method (DeLeón and Loza, 2019). The corrected firstorder reliability method (FORM) and importance sampling were used in this study to overcome invariance nature of FOSM method.
The effect of climate change due to the wavebreaking limit was negligible in the cases studied by Goda et al. (2000). They reported that sliding failure is likely to happen when the water elevation is more than 16 m in their studied breakwater sections. Water elevation was assumed as a random variable in this study to consider the variation. Shimosako and Takahashi (1999) proposed a procedure to assess the risk of sliding failure due to wave height with returned periods equal to the lifetime of the breakwater. To consider statistical inaccuracies and measurement errors, Goda et al. (2000) proposed 10% as the coefficient of variation (CoV) of wave height, and they ignored shoaling, breaking, diffraction, and refraction in their study for simplicity sake. Takayama and Ikeda (1993) studied the coefficient of friction between a concrete caisson and different bed types and reported a mean bias of 6% for the coefficient of friction based on a prototype analysis. Goda (2000) proposed a CoV of 10% for the friction coefficient (Goda and Takagi, 2000).
Different procedures are used to assess caisson sliding, such as statistical, experimental, and numerical methods. Ota et al. (2014) used a neural network to predict the breakwater performance and cumulative sliding failures of rubblemound breakwaters. Lee et al. (2012) investigated the sliding reliability of caisson breakwaters in the waters of Japan and Korea. They proposed a lognormal distribution for safety factors and obtained the reliability index using MCS. They also applied Chebyshev's inequality to determine the upper limit of the probability of failure by considering the mean value and CoV of safety factors. Kim and Suh (2012) studied the safety factors and reported the target reliability for caisson breakwaters in South Korea. They showed that conventional design methods have acceptable accuracy for designing caisson breakwaters. To obtain an acceptable performance of vertical breakwaters against sliding, Kim and Suh (2014) studied and evaluated different types of breakwater sections using MCS and a reliability index β of 2.33(Kim and Suh (2013), Kim and Suh (2011)). They also reported the effect of different wave heights and water depths that should be considered in the design of breakwaters (Suh et al., 2013). Based on their study on the 30year water level forecast in the port of Hitachinaka in Japan (2012), their reliability assessment indicated that the increases in water elevation and wave height due to climate change are acceptable in vertical breakwater design (Suh et al., 2012). Studies performed on breakwaters in different water elevations showed that the increasing wave height, unlike increasing water elevation, has a greater effect on the performance of vertical breakwaters against sliding failure (Lee et al., 2012). Mase et al. (2013) assessed the expected sliding of breakwaters by applying the level III reliability method. Their results showed that a 10% to 60% increase in sliding might occur after considering the effects of climate change. In a numerical study based on the smoothedparticle hydrodynamic (SPH) method, Akbari and Taherkhani assessed wave interactions with a composite breakwater located on a permeable bed. They reported that neglecting bed permeability gives rise to overestimated caisson sliding (Akbari and Taherkhani, 2019). Ghadirian and Bredmose also used the FORM to investigate the effect of bed slope on extreme waves' force (Ghadirian and Bredmose, 2019). M.koç et al. proposed a new reliabilitybased approach and its application in breakwaters (Koç and Imren Koç, 2020) (Koç, 2009). These methods were verified by the MCS and FORM. Recently, Radfar et al. (2021) introduced a probabilistic approach for the optimum design of rubblemound breakwaters under the joint probability density function (PDF) of wave heights, water levels, and different climate conditions.
In this study, different parameters involved in the sliding of breakwaters, such as wave period and height, coefficient of friction, water elevation, and caisson dimension, were considered random variables with random means and standard deviations. Some uncertainties were dismissed in the previous study due to the assumption of a constant mean and standard deviation. In this study, all means and standard deviations are random. Using this approach, some of the dismissed uncertainties are taken into account.
As mentioned in this chapter, due to the randomness in loads and resistance parameters of caisson breakwaters, a probabilistic method is a good approach for predicting their sliding response to incident waves. Accordingly, the procedure proposed in this paper evaluates the caisson sliding in a real case for Tombak port located on the southern coasts of Iran, one of the most important maritime transit areas in the world. In the following chapters, first, reliability methods, such as the FORM and MCS, are briefly explained, followed by sensitivity analysis and the Bayesian linear regression (BLR). These methods are used to assess caisson sliding due to the wave force. Then, based on the sensitive parameters determined from the sensitivity analysis, a direct solution is proposed for predicting caisson sliding using the BLR technique. The Bayesian inference results are then compared with conventional formulas through reliability analyses.
2 Reliability methods used in this study
Reliabilitybased design with the Bayesian approach is used in this research to consider the randomness and uncertainties of the parameters that participate in the limit state of sliding failure. Here, this method is explained first, and then it is employed in a real case study to assess its failure probability. The outcomes of this method and its performance are compared with conventional methods. Figure 1 shows the proposed methodology of the reliability assessment of caisson breakwaters. Each step in the methodology is discussed in the following sections. The reliability index obtained from each method are compared with one another to further investigate the method. The FORM design point is also used for importance sampling to reduce the number of samples. With this approach, the total number of sampling will be reduced dramatically:
• The FORM is implemented to evaluate the reliability index by mapping all random variables to standard normal space and to calculate the distance of the design point from the origin.
• The corrected FORM is used to consider the possible nonlinearity of the limit state function.
• Monte Carlo importance sampling is used to check the failure probability and accuracy of the FORM. When the limit states have a very steep curvature or more than one extremum, sampling methods are used to assess the answer.
The reliability indexes obtained from these methods are compared to determine the failure probability of Tombak caisson breakwaters against sliding failure:
• The sensitivity of the reliability index to each random variable and the contribution of their uncertainties in the probability of failure are examined using different sensitivity parameters.
• BLR is implemented to extract the responding model for the sliding of Tombak port's caisson breakwaters.
• Finally, the posterior prediction for noninformative priors can be used for different marine issues.
2.1 FirstOrder Reliability Method
The probability of the reliable performance for the expected lifetime of a component or structure (system) with the contribution of different types of random variables and uncertainties is known as reliability.
Loads, geometry, and other structural characteristics are considered random variables. Limit state functions g(x) separate the safety and failure region to calculate the failure probability. The general form of the reliability method can be defined by determining failure boundaries (Figure 2) between two sets of load and resistance. R presents the resistance of the structure, and S stands for the loads applied on the structure in Figure 2. fRQ presents the joint PDF of resistance and loads. Different probability distributions for loads and resistance of a structured form the joint probability distribution models (fRQ) for the two sets of random variables.
The probability for which the combination of random variables vector (X) lies in the failure region is defined as
$$ {P}_{f}=P({\boldsymbol{X}}\in F) $$ (1) $$ {P}_{f}={\int }_{F}{f}_{x}\left({\boldsymbol{x}}\right){\mathrm{d}}x={\int }_{g(x)\le 0}{f}_{x}\left({\boldsymbol{x}}\right)\mathrm{d}x $$ (2) F stands for the failure region, and X stands for the vector of random variables. Because this integral does not have an analytical answer (Oumeraci et al., 2001), other alternative methods are used to solve and estimate the answer.
One of the most common methods for this purpose is the FORM, which finds a design point and linearizes the limit state curve at this point to calculate the distance from the origin (see Figure 2). The probability of failure shown in Eq. (2) (Der Kiureghian, 2005) is predicted by transforming all random variables $ {{\boldsymbol{X}}}_{i} $ to standard normal space $ {{\boldsymbol{U}}}_{i} $ and then linearizing the limit state curve with Taylor approximation, where the distance of the limit state from the origin is presented as the reliability index (see Figure 2b). If all random variables have a normal distribution and the limit state function is linear, then the exact answer will be obtained (Jiao and Moan, 1990).
Standard normal space is symmetric, and the probability density function exponentially decreases from the origin of this space. The transformation to the standard normal space for statistically independent random variables is performed for each random variable separately (Figure 3). Here, the Nataf transformation is used to convert the independent variables to correlated random variables (A Der Kiureghian, 2005). The correlation between the two standard normal variables $ {{\boldsymbol{\rho}}}_{z.ij} $ is derived by solving Eq. (6), where $ {{\boldsymbol{\rho}}}_{ij} $ is the correlation between two nonnormal random variables i and j (Liu and Der Kiureghian, 1986). Thus, when two random variables with a nonnormal distribution transfer to a standard normal space, the correlation between them changes, and the Nataf transformation is implemented to obtain the correlation in the new space (Eqs. (6) and (7)).
$$ {F}_{i}({x}_{i})={\Phi }_{i}({u}_{i}) $$ (3) $$ {u}_{i}{=\Phi }^{1}\left({F}_{i}\left({x}_{i}\right)\right)\text{ and}\,\,\,{x}_{i}={F}^{1}\left(\Phi \left({u}_{i}\right)\right) $$ (4) $ \Phi ({u}_{i}) $ indicates the standard normal cumulative distribution function (CDF).
$$ {{z}_{i}=\Phi }^{1}\left({F}_{i}\left({x}_{i}\right)\right) $$ (5) $$ {{\boldsymbol{\rho}}}_{ij}=\underset{\infty }{\overset{+\infty }{\iint }}\left(\frac{{x}_{i}{\mu }_{i}}{{\sigma }_{i}}\right)\left(\frac{{x}_{j}{\mu }_{j}}{{\sigma }_{j}}\right){\varphi }_{2}\left({z}_{i}. {z}_{j}.{{\boldsymbol{\rho}}}_{z.ij}\right){\mathrm{d}z}_{i}{\mathrm{d}z}_{j} $$ (6) where $ {x}_{i} $=$ {F}^{1}\left(\Phi \left({z}_{i}\right)\right);\,{x}_{j} $=$ {F}^{1}\left(\Phi \left({z}_{j}\right)\right) $; $ \mu\,\text{and}\,\sigma\,\text{are the mean and standard deviation} $, respectively; and $ {\varphi }_{2}\left({z}_{i}. {z}_{j}.{\rho }_{z.ij}\right) $, as shown in Eq. (7).
$$ {\varphi }_{2}\left({z}_{i}. {z}_{j}.{\rho }_{z.ij}\right)= \frac{1}{2\pi \sqrt{1{\rho }_{z.ij}^{2}}}\exp (\frac{{z}_{i}^{2}+{z}_{j}^{2}2{\rho }_{z.ij}{z}_{i} {z}_{j}}{2(1{\rho }_{z.ij}^{2})}) $$ (7) Density calculation out of the area split by the hyperplane has a relation with the distance of the plane from the origin. The design point u* (Eq. (8)) is located on the limit state function G(u) = 0 with a minimum distance from the origin (Der Kiureghian, 2005). The firstorder Taylor approximation is used to linearize the limit state function G(u) at the design point (Eq. (9)). The convergence criterion to check the adequacy of iterations is defined as the ratio between the limit state value at the design point and the limit state value at means. The threshold of the convergence ratio is usually selected as approximately 0.001 (Haldar and Mahadevan, 2000) (Haukaas and Der Kiureghian, 2003). Finding the design point is an iterative procedure, and different methods can be used to do so. The improved Hasofer–Lind–Rockwitz–Fiessler (iHLRF) method (1990) is the latest method used to determine the design point on the limit state function (Der Kiureghian, 2005). The convergence problem of the iHLRF method is the unit step size used for the iterative procedure. The unit step size may be unsuitable for different iterations. The Armijo rule is adopted in the iHLRF method to solve the step size problem. This rule cuts the step size in half until the Merit function (Eq. (13)) accepts the step size where the penalty parameter C should be selected at each step to satisfy the condition of Eq. (14) (A Der Kiureghian, 2005).
$$ {{\boldsymbol{u}}}^{*}=\text{argmin}\left\{\Vert {\boldsymbol{u}}\Vert G\left({\boldsymbol{u}}\right)=0\right\} $$ (8) $$ G(\mathit{\boldsymbol{u}}) \cong G({\mathit{\boldsymbol{u}}^ * }) + \nabla {G^{\rm{T}}}({\mathit{\boldsymbol{u}}^ * })(\mathit{\boldsymbol{u}}  {\mathit{\boldsymbol{u}}^ * })\mathop \to \limits^{G({\mathit{\boldsymbol{u}}^ * }) = 0} G(\mathit{\boldsymbol{u}}) \cong \nabla {G^{\rm{T}}}({\mathit{\boldsymbol{u}}^ * })(\mathit{\boldsymbol{u}}  {\mathit{\boldsymbol{u}}^ * })$$ (9) $$ \alpha =\frac{\nabla G\left({\boldsymbol{u}}\right)}{\Vert \nabla G\left({\boldsymbol{u}}\right)\Vert } $$ (10) $$ G\left(u\right)\cong \Vert \nabla G\left({{\boldsymbol{u}}}^{*}\right)\Vert {\boldsymbol{\alpha }}^{\mathrm{T}}\left({\boldsymbol{u}}{{\boldsymbol{u}}}^{*}\right)=\Vert \nabla G\left({{\boldsymbol{u}}}^{*}\right)\Vert \left({{{\boldsymbol{\alpha }}^{\mathrm{T}}{\boldsymbol{u}}}^{*}\boldsymbol{\alpha }}^{\mathrm{T}}{\boldsymbol{u}}\right) $$ (11) $$ G\left({\boldsymbol{u}}\right)\cong \Vert \nabla G\left({{\boldsymbol{u}}}^{*}\right)\Vert \left({\beta \boldsymbol{\alpha }}^{\mathrm{T}}{\boldsymbol{u}}\right) $$ (12) $$ m\left({\boldsymbol{u}}\right)=\frac{1}{2}{\Vert {\boldsymbol{u}}\Vert }^{2}+c\left(G\left({\boldsymbol{u}}\right)\right) $$ (13) $$ C\ge \frac{\Vert\boldsymbol{u}\Vert }{\Vert \nabla G\left({\boldsymbol{u}}\right)\Vert } $$ (14) $$ {P}_{f}=\Phi (\beta ) $$ (15) The FORM calculation algorithm is shown in. Figure 4
2.2 Reliability Sensitivity Analysis
Some metrics are developed to determine the sensitivity of each random variable or other constants or variables, such as the mean or standard deviation of one variable. These importance factors provide information to define the random variable with more influence on failure probability or reliability index. When the problem and limit state function have many variables that contribute to the reliability index value β, it is hard to determine how many variables act as loads or resistance. It is also complicated to know which variable uncertainties should be reduced to produce a higher reliability index, i.e., a lower probability of failure. Sensitivity is derived through the differentiation of the failure probability or reliability index with respect to the parameters (shown as θ in Eq. (16)) that we want to determine its effect on the reliability index or probability of failure.
$$ \frac{\partial \beta }{\partial \theta }=\frac{\partial \beta }{\partial {u}^{*}}\frac{\partial {u}^{*}}{\partial \theta }=\frac{\nabla {G}^{Italic{T}}}{\Vert \nabla G\Vert }\frac{\partial {u}^{*}}{\partial \theta }={\boldsymbol{\alpha }}^{Italic{T}}\frac{\partial {u}^{*}}{\partial \theta } $$ (16) The $ \boldsymbol{\alpha } $ vector extracted from the tangent line equation on the limit state function is one of the important byproducts of the FORM analysis. Alpha vector elements represent the importance of each standard normal variable ui that is transformed from the original space to calculate the probability of failure. The absolute value of each element of the $ \boldsymbol{\alpha } $ vector shows the importance of that random variable. The positive value represents the load, whereas the negative value represents the capacity action of that parameter in the limit state function. As mentioned earlier, $ \boldsymbol{\alpha } $ shows the contribution of the mapped standard normal variables to the total variance of the limit state function. If the basic random variables are correlated, then this importance vector cannot demonstrate the effect of their correlations on the reliability index or other assessment criteria. The correlation problem in the alpha vector is solved by the Gamma vector. The $ {\boldsymbol{\gamma}} $ vector (Eq. (17)) terms arise from the contribution of variances to measure the importance of basic random variables.
$$ \gamma =\frac{\boldsymbol{\alpha }{{\boldsymbol{J}}}_{u.x}\widehat{{\boldsymbol{D}}}}{\Vert \boldsymbol{\alpha }{{\boldsymbol{J}}}_{u.x}\widehat{{\boldsymbol{D}}}\Vert } $$ (17) $$ \text{where}\,J\,\text{is the Jacobian of the transformation}={{\boldsymbol{J}}}_{u.x}={{L}}_{0}^{1}\text{diag}[{{\boldsymbol{J}}}_{ii}] $$ (18) $$ {{\boldsymbol{J}}}_{ii}=\frac{{f}_{i}({{\boldsymbol{x}}}_{i})}{{\varphi }({{\boldsymbol{u}}}_{i})} i=\text{random variables} $$ (19) $$ {{L}_{0}\text{ is the Cholesky decomposition of the correlation matrix }{R_0},{ i}.{e}., R}_{0}={L}_{0}{L}_{0}^{T} $$ (20) $ \widehat{\boldsymbol D}=\mathrm{diag}\left[\widehat{{\boldsymbol\sigma}_{\boldsymbol i}}\right] $ is the diagonal matrix of the standard deviation of the equivalent normal of random variables. $ {\boldsymbol{\delta}} $ and $ {\boldsymbol{\eta}} $ vectors indicate the importance and the contribution of the mean and standard deviation of variables or uncertainties of random variables in the reliability index.
2.3 Monte Carlo Importance Sampling
In this method, the failure probability integral (Eq. (2)) is solved through the sampling of random variables that contribute to the limit state function. For this purpose, a sufficient number of realizations are simulated according to their probability density function. These random samples are used to calculate the limit state function. The ratio between the number of failed samples (N_{f}) and the total number of realizations (N) indicates the probability of failure.
$$ {N_f} = \sum\limits_{j = 1}^N {1\left( {g\left( {{x_j}} \right)} \right)} $$ (21) $$ {P}_{f\text{MC}}\approx \frac{{N}_{f}}{N} $$ (22) The accuracy of the answer increases with the number of simulations. The total number of simulations is a great challenge in the MCS. Some criteria are used to determine the total number of realizations required for a good approximation. This parameter is obtained by minimizing the CoV (Eq. (23)) with respect to the acceptable error (Eq. (24)). This minimization yields the total number of sampling (Eq. (25)). Clearly, a small CoV leads to a higher number of samples.
$$ V_{P_{f\mathrm{MC}}}\approx\frac1{\sqrt{P_{f\mathrm{MC}}N}} $$ (23) $$ \varepsilon \approx \frac{\frac{{N}_{f}}{N}{P}_{f\text{MC}}}{{P}_{f\text{MC}}} $$ (24) $$ N\approx \frac{1}{{{V}_{{P}_{f}}}^{2}}\frac{1{P}_{f\text{MC}}}{{P}_{f\text{MC}}} $$ (25) where $ V_{P_{f\mathrm{MC}}} $ is the coefficient of failure, N is the total number of samples, $ {N}_{f} $ is the total number of failed samples, $ P_{f\mathrm{MC}} $ is the probability of failure estimated by the MCS, and $ {\boldsymbol{\varepsilon}} $ is the error. The MCS algorithm is shown in Figure 5.
The probability of failure calculated by the MCS has a randomness nature, and the CoV of this random variable shows the reliability of this value. If the total number of simulations is infinite, then the CoV will tend to be zero. A CoV value of less than 2% is recommended in the literature (Haldar and Mahadevan, 2000). The failure probability is often small in the design of marine structures and breakwaters. As illustrated in Eq. (25), a large number of samples are required for a small probability of failure with a reasonable CoV. The simulation of a large number of random variables requires a high computational cost and is also a timeconsuming procedure. As explained above, the high number of simulations is the weakness of the direct MCS. The variance reduction method can be used to reduce the number of required simulations. In this method, sampling points are concentrated in the most probable region instead of spreading random variables among the whole domain. In this study, sample reduction is achieved by concentrating the simulation near the design point of the FORM instead of a wide possible range of random variables.
3 BLR Procedure for Random Variable Estimation
BLR models are implemented to predict the posterior of random variables. This method is also known as noninformative priors for multiple parameters (Box and Tiao, 1992). Many parameters are involved in the posterior of the model that is to be extracted from the random variables. If parameters X are the n × k matrix of known constants, then n equations can be written for the regressands, or the model responds y. The model parameters or regression coefficient θ will be obtained by minimizing the square of residuals (minimizing the least square in Eq. (28)) with a variance of $ {{\boldsymbol{\sigma}}}^{2} $. Epistemic and aleatory uncertainties and unknown parameters that are not taken into consideration are involved in the residuals. The mathematical notation of the BLR is explained as follows (Haldar and Mahadevan, 2000; Box and Tiao, 1992):
$$ E\left({\boldsymbol{y}}\right)={\boldsymbol{X}}{\boldsymbol{\theta}} {\Rightarrow }{\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{\theta}}+{\boldsymbol{\varepsilon}} $$ (26) where $ {\boldsymbol{\varepsilon}} $ is the residual vector
$$ {\Vert {\boldsymbol \epsilon} \Vert }^{2}\approx {\varepsilon }_{1}^{2}+{\varepsilon }_{2}^{2}+\cdots +{\varepsilon }_{n}^{2} $$ (27) $$ \widehat{{\boldsymbol{\theta}}}=\text{argmin}\left\{{\Vert {\boldsymbol{\varepsilon}}\Vert }^{2}\right\} $$ (28) $$ {\boldsymbol{l}}({\boldsymbol{\theta}}{\boldsymbol{\sigma}}.{\boldsymbol{y}})\propto \exp \left[\frac{1}{2{{\boldsymbol{\sigma}}}^{2}}{\left({\boldsymbol{y}}{\boldsymbol{X}}{\boldsymbol{\theta}}\right)}^{Italic{T}}({\boldsymbol{y}}{\boldsymbol{X}}{\boldsymbol{\theta}})\right] $$ (29) $$ {\left({\boldsymbol{y}}{\boldsymbol{X}}{\boldsymbol{\theta}}\right)}^{Italic{T}}\left({\boldsymbol{y}}{\boldsymbol{X}}{\boldsymbol{\theta}}\right)={\left({\boldsymbol{y}}\widehat{{\boldsymbol{y}}}\right)}^{Italic{T}}\left({\boldsymbol{y}}\widehat{{\boldsymbol{y}}}\right)+{\left({\boldsymbol{\theta}}\widehat{{\boldsymbol{\theta}}}\right)}^{Italic{T}}{{\boldsymbol{X}}}^{Italic{T}}{\boldsymbol{X}}\left({\boldsymbol{\theta}}\widehat{{\boldsymbol{\theta}}}\right) $$ (30) where $ \widehat{{\boldsymbol{\theta}}}={\left({{\boldsymbol{X}}}^{Italic{T}}{\boldsymbol{X}}\right)}^{1}{{\boldsymbol{X}}}^{Italic{T}}{\boldsymbol{y}} $ is the vector of the least square estimate of θ and $ \widehat{{\boldsymbol{y}}} $ is the vector of the fitted data. The posterior distribution is determined via Eq. (31):
$$ {\boldsymbol{p}}\left({\boldsymbol{\theta}}{\boldsymbol{\sigma}}.{\boldsymbol{y}}\right)=\frac{\sqrt{{{\boldsymbol{X}}}^{Italic{T}}{\boldsymbol{X}}}}{{\left(2{\boldsymbol{\pi}}{{\boldsymbol{\sigma}}}^{2}\right)}^{{~}^{{\boldsymbol{k}}}\!/\!\!{~}_{2}}}Italic{exp}\left[{\frac{1}{2{{\boldsymbol{\sigma}}}^{2}}\left({\boldsymbol{\theta}}\widehat{{\boldsymbol{\theta}}}\right)}^{Italic{T}}{{\boldsymbol{X}}}^{Italic{T}}{\boldsymbol{X}}\left({\boldsymbol{\theta}}\widehat{{\boldsymbol{\theta}}}\right)\right],i=1, \cdots , k $$ (31) 4 Study Area for Breakwater Assessment
In this study, the reliability analysis is performed for the caisson breakwater of Tombak port, located on the southern coasts of Iran. Tombak port is located at 27.702°N, 52.203°E (Figure 6), and is of great importance for service and export purposes in the northwest part of the Persian Gulf shoreline (Limits of oceans and seas (1953)).
Figure 7 represents the section of Tombak port's caisson breakwater. As reported in previous studies, a 16m water depth is critical for caisson sliding (Goda and Takagi, 2000). The Tombak port's caisson breakwater is located at a water depth of 35 m. Because of the high water depth and uncertainty of marine structures due to the presence of environmental loads, reliability studies have been implemented for the sliding failure mode to assess the performance of this breakwater. According to previous studies and failures that happened in the caisson breakwaters, sliding is reported as the most probable failure mode (Oumeraci et al., 2001).
Figure 7 shows the section of breakwater assessed in this study. The schematic view of the dynamic pressure on the caisson is shown in Figure 8. Goda's formulation is the most used formula for calculating wave force on the vertical wall and the amount of sliding (Goda, 2010; OCDI, 2002). Therefore, Goda's formulation is implemented in this study to extract the total induced forced on the caisson as in the following equations:
$$ {P}_{1}=\frac{1}{2}\left(1+\cos \beta \right)\times $$ (32) $$ (\left(0.6+{\left[\frac{4\pi h/L}{\sinh \left(4\pi h/L\right)}\right]}^{2}\right)+\left(\min \left\{\frac{{h}_{b}d}{3{h}_{b}}{\left(\frac{{H}_{\max }}{d}\right)}^{2}\right\}.\frac{2d}{{H}_{\max }}\right){\cos }^{2}\beta )\rho g{H}_{\max } $$ $$ {P}_{2}=\frac{{P}_{1}}{\cosh \left(2\pi h/L\right)} $$ (33) $$ {P}_{3}=\left(1\frac{{h}^{^{\prime}}}{h}[1\frac{1}{\cosh \left(2\pi h/L\right)}]\right){P}_{1} $$ (34) $$ {P}_{u}=\frac{1}{2}(1+\cos \beta )\times \left(0.6+{\left[\frac{4\pi h/L}{\sinh \left(4\pi h/L\right)}\right]}^{2}\right)\left(1\frac{{h}^{^{\prime}}}{h}[1\frac{1}{\cosh \left(2\pi h/L\right)}]\right)\rho g{H}_{\max } $$ (35) $$ {\eta }^{*}=0.75(1+\cos \beta ){H}_{\max } $$ (36) $$ L=\frac{g{T}^{2}}{2\pi }\tanh \frac{2\pi d}{L} $$ (37) T is the wave period, $ {H}_{\max } $ is the maximum wave height, β is the wave obliquities degree, and ρ is the density of the seawater (Goda, 2010; Takahashi, 2010). The horizontal wave load is calculated by summing the wave pressure acting on the walls. The pressures can be calculated as presented in Eqs. (32), (33), and (34). The uplift load is calculated according to Eq. (35) (see Figure 8). The concept of the limit state function for sliding is as follows: the horizontal load should be smaller than the buoyant weight of the caisson, from which the uplift pressure is deducted by considering the friction coefficient.
(Buoyant caisson weight − Uplift pressure) × Friction coefficient ≥ Wave horizontal force.
In the following section, all parameters of the limit state function are explained and considered random variables with random means and standard deviations.
5 Input parameters
The value of caisson sliding depends on several geometric, material, and environmental parameters. Among these parameters, the friction coefficient between the caisson and bed layer, water and concrete density, caisson width and its wet zone, total water depth, wave height, and wave period is taken into account in this study, as listed in Table 1. The most challenging parameter is the wave height because the mean value of other parameters can be easily defined based on the caisson dimension and material.
Random var Distribution type Mean Standard deviation Friction coefficient Lognormal Normal (0.5, 0.1) Normal (0.1, 0.1) Water density Lognormal Normal (1030, 10.3) Normal (10.3, 1) Concrete density Lognormal Normal (2000, 100) Normal (200, 20) Caisson width (B) Lognormal Normal (14, 0.5) Normal (0.28, 0.1) Caisson wet zone (h_{W}) Lognormal Normal (17.7, 1.75) Normal (0.708, 0.1) Total water depth (h) Lognormal Normal (35, 1.75) Normal (1.89, 0.1) Wave height (H_{S}) Lognormal Normal (4.64, 0.5) Normal (0.464, 0.1) Total number of waves (N) Lognormal Normal (2335, 233) Normal (233, 10) Wave period (T) Lognormal Normal (9.25, 1) Normal (0.4625, 0.1) To obtain the wave characteristics, first, longterm analysis is performed using 27year recorded wave data at a depth of 70 m in front of the studied site. To determine the probability of nonexceedance of any significant wave height, the lognormal and Weibull distributions are used (Holthuijsen, 2007). Then, by modeling the wave propagation, the significant wave height is evaluated as 4.64 m at the breakwater location. The maximum wave height is required to obtain the sliding occurrence that defines the sliding or nonsliding condition of the caisson. According to the Rayleigh distribution, the maximum wave height can be calculated based on the significant wave height and number of waves during a storm via Eq. (38), where N is the total number of the wave.
$$ {mod}({H}_{\mathrm {max}})\approx {H}_{S}\sqrt{{ln}({{N}})/2} $$ (38) Because the mean wave period is 9.25 s and the estimated storm duration is 6 h in the area studied, the mean value for the total number of the wave should be 2335, so the maximum wave height will be nearly twice the significant wave height (Holthuijsen, 2007). Meanwhile, to include the uncertainty of this value, a lognormal distribution is implemented with the normal mean and standard deviation.
Finally, the random parameters considered in the sliding assessments for the Tombak port are shown in Table 1. The negative value of the physical parameters, such as width and density, is irrational, so the lognormal distribution is assumed for these parameters. The reliability calculation is performed by considering this assumption. All statistical characteristics of the random variables have a constant value in classic inference, but Bayesian inference is used in this study. Therefore, random variables are considered for the statistical properties of each parameter, where all means and standard deviations are taken as random variables with a normal probability distribution.
6 Results and discussion
The FORM beta index is calculated according to the procedure explained in the previous sections. The nonlinearity of the limit state function affects the calculated failure probability by the FORM. This nonlinearity can change the real probability of failure because of the type of hyperplane fitted at the design point in this method instead of the real hyperbolic limit state function. This problem can be solved by implementing the secondorder method. However, when the limit state has many random variables, the calculation of the secondorder reliability index will be complicated. The correction of the FORM with the first principal curvature is a good solution with enough accuracy, and in this study, the reliability index obtained through this method is approved by the Monte Carlo importance sampling method.
The Monte Carlo importance sampling method is also used for the verification of the FORM. Sampling is performed until the CoV is small enough (approximately less than 2%). Importance sampling is conducted by shifting the origin to the design point of the FORM in the Monte Carlo method. This method significantly reduces the total number of samples. The total number of sampling is more than 600000 if the sampling starts from the mean of each random variable. However, by starting the sampling around the design point of the FORM, the total number is reduced to less than 8000. All reliability analyses were performed using R_{t} (Mahsuli and Haukaas, 2013), a computer program for probabilistic analysis. The reliability index and probability of failure are illustrated in Table 2.
Methods β P_{f} No FORM 2.611 4.52×10^{3} Corrected FORM 2.624 4.34×10^{3} Monte Carlo 2.605 4.59×10^{3} 7443 Table 2 reports the failure probability obtained from the reliability methods, where the performance of the Tombak breakwater against sliding is in the acceptable range. The reliability indexes of the FORM and Monte Carlo method validate each other's values. The corrected FORM reliability index shows that the nonlinearity of the limit state function is slight, so the FORM results did not significantly change.
The sensitivity of each random variable is obtained by the derivation of the reliability index with respect to each random variable. Table 3 shows the sensitivities of each random variable. The results show that the reliability index is sensitive to the coefficient of friction, wave height, and caisson weight (or concrete density). The coefficient of friction is one of the most dominant variables in the sliding of caisson breakwaters, according to the results in Table 3. The probability of failure also shows the sensitivity to this random variable and its uncertainties (standard deviation).
Random variable $ \alpha$ $\gamma $ $ \delta$ $\eta $ $\beta $ sensitivity to the mean $\beta $ sensitivity to the standard deviation Caisson width (B) −0.078 −0.078 0.078 −0.017 0.280 −0.062 H_{s} 0.175 0.175 −0.172 −0.071 −0.740 −0.307 N 0.022 0.022 −0.023 0.001 0.000 0.000 Wave period (T) 0.166 0.166 −0.163 −0.064 −0.353 −0.138 h −0.110 −0.110 0.112 −0.038 0.059 −0.020 h_{W} 0.080 0.080 −0.080 −0.014 −0.113 −0.019 Coefficient of friction −0.774 −0.774 1.118 −1.682 11.176 −16.822 Concrete density −0.560 −0.560 0.648 −0.870 0.003 −0.004 Water density 0.056 0.056 −0.056 −0.008 −0.005 −0.001 y sampling and using each random variable as inputs in the limit state function, the sliding failure CDF and PDF are obtained. The CDF and PDF diagrams of the sliding occurrence are presented in Figure 9a and b. The CDFs of the Gumbel, gamma, and normal distributions are plotted beside the CDF of the sliding occurrence in Figure 9c. This function tends to be a Gumbel distribution with a mean of 1.03 and CoV of 0.57. The probability density function can be used instead of the limit state function in the probabilistic study of caisson sliding to simplify and speed up the calculation and can also be utilized as an initial estimation.
BLR is implemented (Gardoni et al., 2002; Box and Tiao, 1992) to analyze the data accumulated through the sliding limit state sampling. These data are used to extract the reduced model of the sliding occurrence model. The analysis shows that the model posterior for the sliding of caisson breakwaters has a mean and standard deviation of 0.039 and 0.022, respectively. By using the proposed regression model, the safety factor against sliding is easily obtained via Eq. (39) by substituting the coefficients of Table 4. Based on the input parameters, when the equation result is more than zero, it presents a nonsliding condition, and values less than zero present a sliding occurrence. $ {{\boldsymbol{X}}}_{{\boldsymbol{i}}} $ is the random variable value involved in the limit state function. Their coefficients and correlations are shown in Table 4.
$$ \text{sliding occurrence}={\boldsymbol{Y}}={{\boldsymbol{\theta}}}_{i}\times {{\boldsymbol{X}}}_{i}+{\boldsymbol{\varepsilon}}\ {\boldsymbol{\theta}}=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\theta 1\\ \theta 2\end{array}\\ \theta 3\end{array}\\ \theta 4\end{array}\\ \theta 5\end{array}\\ \theta 6\\ \theta 7\end{array}\\ \theta 8\end{array}\right]\ {\boldsymbol{X}}=\left(\begin{array}{ccc}{x}_{11}& \cdots & {x}_{18}\\ \vdots & \ddots & \vdots \\ {x}_{n1}& \cdots & {x}_{n8}\end{array}\right)\ {\boldsymbol{\varepsilon}}=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{\varepsilon }_{1}\\ {\varepsilon }_{2}\end{array}\\ {\varepsilon }_{3}\end{array}\\ {\varepsilon }_{4}\end{array}\\ {\varepsilon }_{5}\end{array}\\ {\varepsilon }_{6}\\ {\varepsilon }_{7}\end{array}\\ {\varepsilon }_{8}\end{array}\right] $$ (39) The adequacy of the proposed regression model, which represents the relationship between the regressors and response variable, is evaluated using the coefficient of determination or R2 factor. The R factor has a value between 0 and 1, and the closer it gets to 1 implies that the regression model provides good predictions. The information of this coefficient should be used with caution when the model has many regressors. In this case, it is preferable to perform a residual analysis. Plotting residuals are very informative. If the residuals have a horizontal band on both sides of zero, then their mean values are approximately zero, and the variance is constant around the regression line (Haldar and Mahadevan, 2000).
The R factor for the newly proposed formulas is close to 1. Although this is one of the proper conditions for an accurate model, it is not sufficient, and further controls are required to approve the model adequacy. In Figure 10, the Bayesian regression sliding occurrence model is plotted versus the sliding occurrence based on Goda's formulation prediction. The good compatibility between the results shows the competent prediction of the proposed model. Of note, the new model is simpler to be utilized in comparison to the conventional model. The Q–Q plot in Figure 11 shows the quantiles of the new formula and conventional method and approves the adequacy of the new model. The normal quantile plot or Q–Q plot is a graphical method used to compare two probability distributions by plotting their quantiles against each other. Residual analysis is another method used in this research to assess the appropriateness of BLR data (Figure 12). The plotted data in this figure represent a uniform distribution of residuals on both sides of zero, which means that all random variables' residuals are homoscedastic.
7 Conclusions
Breakwater performance is crucial for marine transportation, and there are many uncertainties in the marine environment. Reliability analysis is an approach used to assess each criteria's performance. The method proposed here can be used to design and assess existing breakwaters and detect failure probability. The FORM and MCS were used in this study to examine the sliding occurrence of the Tombak caisson breakwater in the Persian Gulf. The probability of failure was obtained due to the longterm maximum wave height. By means of the sensitivity analysis, the parameters used in the calculation were categorized to the load and resistance, and the sensitivity of the reliability index to each parameter was procured. Considering the beta index calculated with the Monte Carlo method, the Bayesian reliability analysis showed a small nonlinearity in the sliding limit state. Some of these nonlinear properties were taken into account in the FORM by considering the first principal curvature in the limit state function. The Monte Carlo importance sampling method for the sliding of caissons was also used to assess the reliability index, and its results confirmed that the beta index was secured by the corrected FORM. The beta index also demonstrated an acceptable probability of failure for the sliding of caissons in Tombak port (OCDI, 2002). The reliability sensitivity analysis illustrated the importance of the coefficient of friction and the effect of its uncertainties (aleatory and epistemic uncertainties) on the reliability index. The reliability sensitivity analysis of the standard deviation showed that the failure probability was highly sensitive to uncertainties, such as wave height and friction coefficient. The alpha sensitivity vector revealed the importance of the wave heights and weight of concrete blocks after the coefficient of friction. Based on the sensitivity analysis performed in this study, bed friction is one of the major design parameters of caisson breakwaters. Caisson sliding is very sensitive to the mean value and standard deviation of the coefficient of friction, so an accurate evaluation of this parameter in the design process is recommended to reduce uncertainties.
Bayesian reliability assessment and Bayesian regression, which are noninformative priors for multipleparameter procedures, are also explained in this paper. The regression coefficient of determination and residuals analysis demonstrated good accuracy in the model response. Residuals were normally distributed around zero with no heteroskedastic observation. This method was implemented to predict the posterior of the sliding occurrence. The regression model used in this study was unbiased, and it explicitly calculated the sliding occurrence of the caisson breakwater to facilitate its application. With this method, the design procedure can be significantly shortened. In this study, Goda's formula, which has a longer computational process than the proposed regression formula, was summarized to obtain sliding probability. For various engineering topics, this procedure can be used in conjunction with laboratory studies to reduce the computation time.
Acknowledgements: Tombak breakwater data was provided by Pars Geometry Consultants, and their support is fully acknowledged. 
Figure 3 Mapping independent random variables to standard normal space (Der Kiureghian, 2005)
Figure 8 Parameter used in the calculation of the wave's force acting on the vertical wall (CEM, 2002)
Table 1 Random variables used for the probabilistic study
Random var Distribution type Mean Standard deviation Friction coefficient Lognormal Normal (0.5, 0.1) Normal (0.1, 0.1) Water density Lognormal Normal (1030, 10.3) Normal (10.3, 1) Concrete density Lognormal Normal (2000, 100) Normal (200, 20) Caisson width (B) Lognormal Normal (14, 0.5) Normal (0.28, 0.1) Caisson wet zone (h_{W}) Lognormal Normal (17.7, 1.75) Normal (0.708, 0.1) Total water depth (h) Lognormal Normal (35, 1.75) Normal (1.89, 0.1) Wave height (H_{S}) Lognormal Normal (4.64, 0.5) Normal (0.464, 0.1) Total number of waves (N) Lognormal Normal (2335, 233) Normal (233, 10) Wave period (T) Lognormal Normal (9.25, 1) Normal (0.4625, 0.1) Table 2 Results for the failure probability
Methods β P_{f} No FORM 2.611 4.52×10^{3} Corrected FORM 2.624 4.34×10^{3} Monte Carlo 2.605 4.59×10^{3} 7443 Table 3 Sensitivities of the random variables in the FORM
Random variable $ \alpha$ $\gamma $ $ \delta$ $\eta $ $\beta $ sensitivity to the mean $\beta $ sensitivity to the standard deviation Caisson width (B) −0.078 −0.078 0.078 −0.017 0.280 −0.062 H_{s} 0.175 0.175 −0.172 −0.071 −0.740 −0.307 N 0.022 0.022 −0.023 0.001 0.000 0.000 Wave period (T) 0.166 0.166 −0.163 −0.064 −0.353 −0.138 h −0.110 −0.110 0.112 −0.038 0.059 −0.020 h_{W} 0.080 0.080 −0.080 −0.014 −0.113 −0.019 Coefficient of friction −0.774 −0.774 1.118 −1.682 11.176 −16.822 Concrete density −0.560 −0.560 0.648 −0.870 0.003 −0.004 Water density 0.056 0.056 −0.056 −0.008 −0.005 −0.001 
Akbari H, Taherkhani A (2019) Numerical study of wave interaction with a composite breakwater located on permeable bed. Coast Eng 146: 1–13. https://doi.org/10.1016/j.coastaleng.2018.12.006 Box GEP, Tiao GC (1992). Bayesian inference in statistical analysis. John Wily and Sons, Inc., 25–60 https://doi.org/10.1002/9781118033197 CEM (2002). Coastal engineering manual. Part VI  Chapter 5  Fundamentals of Design. U. S Army Corps Eng. Chaudhary B, Hazarika H, Ishibashi I, Abdullah A (2017) Sliding and overturning stability of breakwater under combined effect of earthquake and tsunami. Ocean Eng 136: 106–116. https://doi.org/10.1016/j.oceaneng.2017.03.021 Cornell CA (1969). Probabilitybased structural code. Am Concr. InstJ, 66, 974–985. https://doi.org/10.14359/7446 DeLeón D, Loza L (2019) Reliabilitybased analysis of Lázaro Cárdenas breakwater including the economical impact of the port activity. Int J Disaster Risk Reduct 40: 101276. https://doi.org/10.1016/j.ijdrr.2019.101276 Der Kiureghian A (2005) Firstand secondorder reliability methods. CRC Press, Boca Raton, USA, Engineering design reliability. https://doi.org/10.1201/9780203483930 Faber MH (2007). Risk and safety in civil engineering. Lecture notes. Inst. Struct. Eng., 335. https://doi.org/10.3929/ethza004230964 Faber MH, Sørensen JD (2002). Reliability based code calibration joint committee on structural safety. Proceedings of the 9th International Conference on Applications of Statistics and Probability, San Francisco, 1–17 Gardoni P, Der Kiureghian A, Mosalam KM (2002) Probabilistic capacity models and fragility estimates for reinforced concrete columns based on experimental observations. J Eng Mech 128: 1024–1038. https://doi.org/10.1061/(ASCE)07339399(2002)128:10(1024) Ghadirian A, Bredmose H (2019) Investigation of the effect of the bed slope on extreme waves using first order reliability method. Mar Struct 67: 102627. https://doi.org/10.1016/j.marstruc.2019.05.005 Goda Y (2010). Random seas and design of maritime structures. In: Advanced Series on Ocean Engineering, Advanced Series on Ocean Engineering. World Scientific. https://doi.org/10.1142/3587 Goda Y, Takagi H (2000) A reliability design method of caisson breakwaters with optimal wave heights. Coast Eng J 42: 357–387. https://doi.org/10.1142/s0578563400000183 Haldar A, Mahadevan S (2000) Probability, reliability, and statistical methods in engineering design. John Wiley Sons, New York, pp 138–272 Haukaas T, Der Kiureghian A (2003). Finite element reliability and sensitivity methods for performancebased earthquake engineering. Pacific Earthquake Engineering Research Center, Technical report No. 2003/14. Holthuijsen LH (2007). Waves in oceanic and coastal waters. Cambridge University Press, New York, 73104. https://doi.org/10.1017/CBO9780511618536 International Hydrographic Organization (1953). Limits of oceans and seas. International Hydrographic Organization (IHO), Bremerhaven, Pangaea. 2021. Jiao G, Moan T (1990) Methods of reliability model updating through additional events. Struct Saf 9: 139–153. https://doi.org/10.1016/01674730(90)90005A Kim SW, Suh KD (2011) Evaluation of target reliability indices and partial safety factors for sliding of caisson breakwaters. J Coastal Res 64: 622–626 Kim SW, Suh KD (2013) Performance analysis of vertical breakwaters designed by partial safety factor method. J Coast Res 65: 296–301. https://doi.org/10.2112/si65051.1 Koç ML (2009) Risk assessment of a vertical breakwater using possibility and evidence theories. Ocean Eng 36: 1060–1066. https://doi.org/10.1016/j.oceaneng.2009.07.002 Koç ML, Imren Koç D (2020) A cloud theory based reliability analysis method and its application to reliability problems of breakwaters. Ocean Eng 209: 107534. https://doi.org/10.1016/j.oceaneng.2020.107534 Lee CE, Kim SW, Park DH, Suh KD (2012) Target reliability of caisson sliding of vertical breakwater based on safety factors. Coast Eng 60: 167–173. https://doi.org/10.1016/j.coastaleng.2011.09.005 Li MS, Hsu CJ, Hsu HC, Tsai LH (2020) Numerical analysis of vertical breakwater stability under extreme waves. J Mar Sci Eng 8: 1–15. https://doi.org/10.3390/jmse8120986 Liu PL, Der Kiureghian A (1986) Multivariate distribution models with prescribed marginals and covariances. Probabilistic Eng Mech 1: 105–112. https://doi.org/10.1016/02668920(86)900330 Mahsuli M, Haukaas T (2013) Computer program for multimodel reliability and optimization analysis. J Comput Civ Eng 27: 87–98. https://doi.org/10.1061/(ASCE)CP.19435487.0000204 Mase H, Tsujio D, Yasuda T, Mori N (2013) Stability analysis of composite breakwater with wavedissipating blocks considering increase in sea levels, surges and waves due to climate change. Ocean Eng 71: 58–65. https://doi.org/10.1016/j.oceaneng.2012.12.037 Nowak A, Collins K (2012). Reliability of structures. McGrawHill High. Educ., 50–60 OCDI (2002). Technical standards and commentaries for port and harbor facilities in Japan. The Overseas Coastal Area Development Institute of Japan. Ota T, Matsumi Y, Kawamura H, Ohno K (2014) Prediction model for change in performance of rubble mound revetment under damage progression. Coastal Engineering Proceedings 1(34): 63. https://doi.org/10.9753/icce.v34.structures.63 Oumeraci H (1994) Review and analysis of vertical breakwater failures  lessons learned. Coast Eng 22: 3–29. https://doi.org/10.1016/03783839(94)900469 Oumeraci H, Kortenhaus A, Allsop W, de Groot M, Crouch R, Vrijling H, Voortman H (2001). Probabilistic design tools for vertical breakwaters, CRC Press, Vol2d 1–17. Radfar S, Shafieefar M, Akbari H, Galiatsatou PA, Mazyak AR (2021) Design of a rubble mound breakwater under the combined effect of wave heights and water levels, under present and future climate conditions. Appl Ocean Res 112: 102711. https://doi.org/10.1016/j.apor.2021.102711 [34] Takahashi S (2010). Design of vertical breakwaters. Port and Airport Research Institute (PARI). Japan report No. 34. https://doi.org/10.1142/9789814282413_0004 Toyama S (1985) Application of the reliability design method to the safety of caisson breakwaters against sliding. Technical Note of Port and Harbour Research Institute 540: 49 Shimosako K, Takahashi S (1999). Application of deformationbased reliability design for coastal structures  espected sliding distance method of composite breakwaters. Proceedings of the International Conference Coastal Structures, 363–371. Suh KD, Kim SW, Kim S, Cheon S (2013) Effects of climate change on stability of caisson breakwaters in different water depths. Ocean Eng 71: 103–112. https://doi.org/10.1016/j.oceaneng.2013.02.017 Suh KD, Kim SW, Mori N, Mase H (2012) Effect of climate change on performancebased design of caisson breakwaters. J. Waterw. Port. Coast Ocean Eng 138: 215–225. https://doi.org/10.1061/(ASCE)WW.19435460.0000126 Takayama T, Ikeda N (1993) Estimation of sliding failure probability of present breakwaters for probabilistic design. Report of the Port and Harbour Research Institute 31(5): 3–32 van der Meer JW (1988) Deterministic and probabilistic design of breakwater armor layers. J. Waterw. Port. Coast Ocean Eng 114:66–80. https://doi.org/10.1061/(ASCE)0733950X(1988)114:1(66)