Level III Reliability Design of an Armor Block of Rubble Mound Breakwater Using Probabilistic Model of Wave Height Optimized for the Korean Sea Wave Conditions and Non-Gaussian Wave Slope Distribution

: In this study, a Level III reliability design of an armor block of rubble mound breakwater was developed using the optimized probabilistic wave height model for the Korean marine environment and Van der Meer equation. To demonstrate what distinguishes this study from the others, numerical simulation was ﬁrst carried out, assuming that wave slope follows Gaussian distribution recommended by PIANC. Numerical results showed that Gaussian wave slope distribution overpredicted the failure probability of armor block, longer and shorter waves, and on the contrary, underpredicted waves of the medium period. After noting the limitations of Gaussian distribution, some efforts were made to develop an alternative for Gaussian distribution. As a result, non-Gaussian wave slope distribution was analytically derived from the joint distribution of wave amplitude and period by Longuet–Higgins using the random variables transformation technique. Numerical results showed that non-Gaussian distribution could effectively address the limitations of Gaussian distribution due to its capability to account for the nonlinear resonant wave–wave interaction and its effects on the wave slope distribution that signiﬁcantly inﬂuences the armor block’s stability. Therefore, the non-Gaussian wave slope distribution presented in this study could play an indispensable role in addressing controversial issues such as whether or not enormous armor blocks like a Tetrapod of 100 t frequently mentioned in developing countermeasures against rough seas due to climate change is too conservatively designed.


Introduction
Water waves generated offshore by various mechanisms propagate toward nearby beaches, and during this process, waves become irregular due to the variability inherent in the marine environment. When such irregular waves enter the shallow water, waves undergo drastic changes due to shoaling, refraction, and diffraction, and as a result, wave height grows, leading to the breaking down of waves at the final stage of shoaling process [1]. Port's outer facilities are deployed in such an irregular and harsh marine environment and, therefore, should be designed to secure sufficient durability even in a harsh marine environment. Nonetheless, armor blocks of the port's outer facilities have been designed using a deterministic design method based on the Hudson equation [2] in Korea. Considering its simple structural form, which makes its practical applications more amenable, the Hudson equation's popularity is understandable. However, it has flaws such that the structural form of Hudson's equation and its parameters cannot reflect sea wave conditions accurately enough for the optimal design of armor block. It is worthy of note that the wave period of incoming waves and storm waves' duration period were overlooked in Hudson's equation as well. The limitation of the deterministic design method is 2 of 24 on full display on the enormous armor blocks such as a TTP of 100 t frequently mentioned in developing countermeasures against rough seas due to climate change, leading to a critique raised in the coastal engineering community that the port's outer facilities in Korea are too conservatively designed.
An alternative that can address the limitations of the deterministic design method seems to be reliability design. In the reliability design, the degree of safety can systematically be evaluated by statistically describing the uncertainty inherent in the armor block's hydraulic stability strength and wave force acting on the armor block. Moreover, by allowing a reasonable failure probability of a armor block, design factors that can be conservatively evaluated in a deterministic design can be easily specified, leading to significant cost savings. Considering the irregularities inherent in sea wave conditions, reliability design should be carried out based on long-term situ wave data. However, in Korea, situ wave data continuously measured over a long period are not available in most cases. As a result, reliability design was performed based on the probability model preferred in the United States [3][4][5]. However, according to Cho and Kim [6] and Cho et al. [7], the Korean marine environment has a significant structural difference from Europe or the United States. This difference is on full display in the East Sea, where enough fetch length required for the full development of wind waves is secured only in a limited direction due to its counterclockwise twisted shape. The rationale for this reasoning can be found in the fact that, unlike the United States, where swells of a long period are prevailing since short waves are left behind the wave packet due to its slow celerity over the course of propagation, swells of a wide variety of periods are occurring on the east coast of Korea [6][7][8].
In light of the facts mentioned above, the reliability design of a port's outer facilities should be implemented based on the optimized probability model for the Korean marine environment, but these efforts are very rare. Based on the Van der Meer equation [9][10][11], wave slope and wave height dictate the stability of the armor block of rubble mound breakwater. For short-term wave height, the Rayleigh distribution for the offshore, the twoparameter Weibull distribution [12] and modified Glukhovskiy distribution [13] for shallow waters are most frequently mentioned in the current literature [1,7,14]. On the other hand, for the long-term wave height, the tri-variate Weibull distribution is preferred [1,7,14]. Lately, some efforts to develop a probabilistic model optimized for the Korean marine environment was made by Choi and Cho [1]. In this study, Choi and Cho [1] first hindcasted the significant wave heights and peak periods off the Ulsan harbor every hour from 1 January 2003 to 31 December 2017 based on the meteorological data by Japan Meteorological Agency (JMA) and National Oceanic and Atmospheric Administration (NOAA), and SWAN. Then, Choi and Cho [1] proceeded to derive a long-term probabilistic model from the numerically simulated wave data using maximum likelihood method (MLM) and showed that the agreements were more remarkable in the modified Glukhovskiy distribution than in the three variates Weibull distribution, which has been preferred in the literature.
Unlike wave height distribution, for which considerable research achievement has been made over the last decades, wave slope distribution has drawn less attention despite its great engineering value. Due to the lack of a robust probabilistic model, in Korea, wave slope has been empirically estimated in reliability design practice such as the following: Wave height is first simulated using Monte Carlo simulation based on the probabilistic model mentioned above, and then wave period is evaluated from simulated wave height using a simple empirical relationship. Finally, wave slopes are estimated from these simulated wave heights and empirically evaluated periods. In doing so, wave period has been assumed to be linearly proportional to wave height [3][4][5]. However, according to Cavanie et al. [15], Longuet-Higgins [16], Park and Cho [17], as sea wave conditions become harsh, the correlation between wave height and its associated wave period is progressively weakened. As a result, wave height and its associated wave period behave as a mutually independent random process. These behavioral characteristics imply that reliability design based on empirical relationship telling that wave period is linearly proportional to wave height can cause significant errors. On the contrary, wave slope is assumed to follow Gaussian distribution [18] in Europe or the United States. However, Gaussian distribution presents limitations, such as imposing a no-negligible probability of occurrence even for negative wave slope. Therefore, it is subject to extensive test in the near future whether Gaussian distribution preferred in Europe or the United States can be applied to the Korean marine environment and harsh sea wave conditions due to climate change.
In this rationale, this study aims to develop a Level III reliability design method for the armor block of rubble mound breakwaters based on the Van der Meer equation and a probabilistic model optimized for the Korean marine environment. In doing so, the armor block's failure probability is analytically derived from the Van der Meer equation using the transformation technique of random variables [19]. Following Choi and Cho [1], the modified Glukhovskiy distribution was used as a probabilistic model for wave height. In the case of wave slope, in an effort to develop an alternative to Gaussian distribution, non-Gaussian wave slope distribution is analytically derived using the random variable transformation technique from the joint distribution of wave amplitude and period by Longuet-Higgins [16], which is the most frequently mentioned in the current literature. The non-Gaussian wave slope distribution presented in this study was verified using the situ wave data collected between 26 April 2017-20 April 2018, using an ultrasonic wave height meter (Nortek) near the Mang-Bang beach (129 • 13 34.56 E, 37 • 24 11.22 N) (see Figure 1) on the east coast of Korea where the water depth was 26.5 m [6,7]. J. Mar. Sci. Eng. 2021, 9,223 3 of 28 independent random process. These behavioral characteristics imply that reliability design based on empirical relationship telling that wave period is linearly proportional to wave height can cause significant errors. On the contrary, wave slope is assumed to follow Gaussian distribution [18] in Europe or the United States. However, Gaussian distribution presents limitations, such as imposing a no-negligible probability of occurrence even for negative wave slope. Therefore, it is subject to extensive test in the near future whether Gaussian distribution preferred in Europe or the United States can be applied to the Korean marine environment and harsh sea wave conditions due to climate change.
In this rationale, this study aims to develop a Level III reliability design method for the armor block of rubble mound breakwaters based on the Van der Meer equation and a probabilistic model optimized for the Korean marine environment. In doing so, the armor block's failure probability is analytically derived from the Van der Meer equation using the transformation technique of random variables [19]. Following Choi and Cho [1], the modified Glukhovskiy distribution was used as a probabilistic model for wave height. In the case of wave slope, in an effort to develop an alternative to Gaussian distribution, non-Gaussian wave slope distribution is analytically derived using the random variable transformation technique from the joint distribution of wave amplitude and period by Longuet-Higgins [16], which is the most frequently mentioned in the current literature. The non-Gaussian wave slope distribution presented in this study was verified using the situ wave data collected between 26 April 2017-20 April 2018, using an ultrasonic wave height meter (Nortek) near the Mang-Bang beach (129 13'34.  This article is structured as follows: for the sake of self-containment, Section 2 provides details about the deterministic design of armor block of rubble mound breakwaters based on the Van der Meer equation and Level I, II reliability design method. In Section 3, a probabilistic wave height model optimized for the Korean marine environment by Choi and Cho [1] is presented. The derivation procedure of non-Gaussian wave slope distribution using the random variable transformation technique from the joint distribution of wave amplitude and period by Longuet-Higgins [16] is also presented in Section 3. In Section 4, Level III reliability-based design of armor block of rubble mound breakwater is presented, and in Section 5, numerical results are provided. This article is structured as follows: for the sake of self-containment, Section 2 provides details about the deterministic design of armor block of rubble mound breakwaters based on the Van der Meer equation and Level I, II reliability design method. In Section 3, a probabilistic wave height model optimized for the Korean marine environment by Choi and Cho [1] is presented. The derivation procedure of non-Gaussian wave slope distribution using the random variable transformation technique from the joint distribution of wave amplitude and period by Longuet-Higgins [16] is also presented in Section 3. In Section 4, Level III reliability-based design of armor block of rubble mound breakwater is presented, and in Section 5, numerical results are provided.

Deterministic Design Based on the Van der Meer Equation
Van der Meer [18] carried out a series of hydraulic model tests at Delft Hydraulics, retrieving the pioneering study of Thompson and Shuttler [20]. These hydraulic model tests were the most comprehensive ones, including various wave conditions and permeability characteristics of the armor and core block. Based on the data accumulated in the process, Van der Meer [9][10][11]18] proposed an empirical equation for plunging and surging breaker.
In the case of rubble mound breakwater of 1:1.5 slope, the size of tetrapod (TTP) is given by (see Figure 2):  Thompson and Shuttler [20]. These hydraulic model tests were the most comprehensive ones, including various wave conditions and permeability characteristics of the armor and core block. Based on the data accumulated in the process, Van der Meer [9][10][11]18] proposed an empirical equation for plunging and surging breaker.
In the case of rubble mound breakwater of 1:1.5 slope, the size of tetrapod (TTP) is given by (see Figure 2): In Equation (1) (2) where random variables are denoted by bold letters.
As shown in Equation (2), when the armor block's hydraulic stability strength R is less than the wave force S , the armor block of rubble mound breakwater is unstable, and otherwise, it is stable.

Reliability Design
In the hierarchy of reliability design, the deterministic design discussed in Section 2.1 can be classified as a Level I reliability design. According to Van der Meer [9,10], it is the ratio of wave height to the nominal diameter of armor block that determines the stability of armor block, and this point of view leads to the Level I partial safety factor method described below.

Partial Safety Factor Method (Level I)
For the sake of convenience in the following analysis, Van der Meer [18] first reclassified six mutually independent random processes into the armor block's hydraulic stability strength and wave force. Then, Van der Meer [18] proceeded to introduce the hydraulic stability strength reduction factor and load factor that account for the irregularities In Equation (1), ∆ = ρ rock /ρ water − 1 is the relative density, D n = (W TTP /ρ rock ) 1/3 is the nominal diameter of armor block, W TTP is the weight of armor block, N z is the relative damage level, S om is the wave slope, N od is the number of waves during a storm, ρ rock and ρ water are the density of armor block and water, respectively.
After rearranging terms in Equation (1) to clarify the physical mechanisms underlying the Van der Meer equation, it can be rewritten as follows: where random variables are denoted by bold letters. As shown in Equation (2), when the armor block's hydraulic stability strength R is less than the wave force S, the armor block of rubble mound breakwater is unstable, and otherwise, it is stable.

Reliability Design
In the hierarchy of reliability design, the deterministic design discussed in Section 2.1 can be classified as a Level I reliability design. According to Van der Meer [9,10], it is the ratio of wave height to the nominal diameter of armor block that determines the stability of armor block, and this point of view leads to the Level I partial safety factor method described below.

Partial Safety Factor Method (Level I)
For the sake of convenience in the following analysis, Van der Meer [18] first reclassified six mutually independent random processes into the armor block's hydraulic stability strength and wave force. Then, Van der Meer [18] proceeded to introduce the hydraulic stability strength reduction factor and load factor that account for the irregularities inherent in the armor block's hydraulic stability strength and wave force rather than describing the irregularities inherent in random variables such as N od , N z , ∆, D n , H S and S om that dictate the armor block's stability by introducing a probability model for each random variable. In this case, the reliability function Z defined as the difference between the armor block's hydraulic stability strength and wave force acting on the structure can be written as: With the given values of the reliability function Z, the state of armor block of rubble mound breakwater can be classified such as the following: In Equation (3), the hydraulic stability strength reduction factor γ Z and load increase factor γ H S can be written as [21]: where T L is a design life cycle of the structure, H S T L is the central estimate of H S with a return period equal to the design life cycle of the structure, H S 3T L is the central estimate of H S with a return period equal to three times the design life cycle T L of the structure, H S T P f is the central estimate of H S with a return period corresponding to the permissible failure probability of the structure, P f is the permissible failure probability, F H s is a random coefficient introduced to account for the short-term variability of H S , having a mean value of 1.0 and the variational coefficient σ F H S [18], T P f is the a return period corresponding to the permissible failure probability of the structure, N is the number of wave data used in the evaluation of probability coefficients of wave height probability distribution, k α , k β , and k s are the coefficients obtainable by carrying out optimization for each failure mode [18].

Level II Reliability Design
While the reliability function nonlinearly depends on the armor block's hydraulic stability strength and wave force, in the Level II reliability design method, the reliability index β around the most probable failure point (MPFP) or design point is first calculated, from which the probability of failure can be obtained. In doing so, the reliability function Z is assumed to follow the Gaussian distribution and linearized around the MPFP such that only the first-order term is retained after expanding the reliability function around the MPFP using a Taylor Series. Here, β denotes the distance from the mean of the reliability function to the limit state and can be written as (see Figure 3): J. Mar. Sci. Eng. 2021, 9, 223 6 of 28 MPFP such that only the first-order term is retained after expanding the reliability function around the MPFP using a Taylor Series. Here, β denotes the distance from the mean of the reliability function to the limit state and can be written as (see Figure 3):

Wave Height Distribution
Rayleigh distribution is the most frequently referred probabilistic model of wave height in the current literature, and despite being a linear model, it shows good agreement with situ data in the case of deep-water random waves [22]. As irregular waves enter the shallow water, wave height distribution undergoes drastic changes due to shoaling, re- MPFP denotes a point where the reliability function's probability density function reaches its peak along the line Z = 0 (limiting state).

Wave Height Distribution
Rayleigh distribution is the most frequently referred probabilistic model of wave height in the current literature, and despite being a linear model, it shows good agreement with situ data in the case of deep-water random waves [22]. As irregular waves enter the shallow water, wave height distribution undergoes drastic changes due to shoaling, refraction, diffraction, and wave breaking. It was Glukhovskiy [23] who first developed wave height distribution at a finite depth. Later, the uncertainty concerned with the evaluation of probability parameters of the early model [24] was addressed by Klopman and Stive [13] and Klopman [25]. As of now, modified Glukhovskiy distribution has been regarded as a representative shallow-water wave height distribution [1]. However, modified Glukhovskiy distribution is known to give conservative values over the surf zone, the last stage of the shoaling process [26]. Lately, Battjes and Groenendijk [12] proposed composite Weibull distribution by classifying the sample space of wave height into two groups to effectively explain wave height attenuation due to wave breaking to address the deficiency of modified Glukhovskiy distribution mentioned above. Later, composite Weibull distribution developed by Battjes and Groenendijk [12] becomes the preferred probabilistic model of wave heights over the surf zone. On the other hand, the three variates of the Weibull distribution were frequently mentioned in the current literature as the long-term wave height distribution.
For the sake of self-containment, Rayleigh distribution, Glukhovskiy distribution, and composite Weibull distribution were presented in Sections 3.1.1-3.1.3, respectively. In Sections 3.1.4 and 3.1.5, the details about the tri-variates Weibull distribution and optimized wave height distribution for the Korean marine environment were provided.

Rayleigh Distribution
Rayleigh distribution F H (h) can be written as [22]: where P[H < h] denotes a probability that wave height H is less than h.

Modified Glukhovskiy Distribution
Modified Glukhovskiy distribution F H (h) can be written as: where A and κ are probability coefficients depending on the relative depth d [= d/h] and can be written as: As shown in Equation (9), modified Glukhovskiy distribution is converging to Rayleigh distribution as d → ∞ .

Composite Weibull Distribution
Composite Weibull distribution F H (h) can be written as [12]: The probability density function f H (h) of composite Weibull distribution can be written as [27]: H RMS is the root mean squared wave height, and the transition wave height H tr = H tr /H RMS can be written as [12]: In Equation (14), d is the water depth, and tan α is the bottom slope. As shown in Equations (12) and (13), composite Weibull distribution is defined by the probability parameters such as κ 1 , κ 2 , H 1 , H 2 , and H tr , and the evaluation method of these parameters can be found in Battjes and Groenendijk [12] and Choi and Cho [1].

Tri-Variates Weibull Distribution
The tri-variates Weibull distribution is known to provide a satisfactory description of the statistical distribution of long-term wave height, and its distribution function F H (h) can be written as [14]: where µ L , H L , and κ L are the probability parameters that are to be determined from the situ wave data using maximum likelihood method (MLM).

Optimized Probabilistic Model of Wave Height for the Korean Sea Wave Conditions
Lately, Choi and Cho [1] developed a probabilistic model optimized for the marine environment in Korea. In this study, Choi and Cho [1] first hindcasted the significant wave heights and peak periods off the Ulsan harbor every hour from 1 January 2003 to 31 December 2017 based on the meteorological data by Japan Meteorological Agency (JMA) and National Oceanic and Atmospheric Administration (NOAA), and SWAN. Then, Choi and Cho [1] derived a long-term probabilistic model of wave height using MLM from the numerically simulated wave data. It was shown that the agreements were more remarkable in the probability distribution in line with modified Glukhovskiy distribution than in the three parameters Weibull distribution, which has been preferred in the literature (see       The wave height probability density function f H (h) by Choi and Cho [1] can be written as: The corresponding wave height distribution function F H (h) can be written as: According to Choi and Cho [1], the parameters of modified Glukhovskiy distribution given in Equation (17) have a value of A p = 15.92, H p = 4.374 m, κ p = 1.824 with a confidence level of 95%.

Gaussian Distribution
For the case of wave slope that significantly influences the armor block's hydraulic stability strength, the Gaussian distribution was preferred in Europe or the United States [18] and can be written as [18,28]: where µ S is the mean wave slope, and σ S is the standard deviation of wave slope. Even though the Gaussian distribution in Equation (18) has been preferred due to its simple structural form, the Gaussian wave slope distribution is known to have flaws such as the non-negligible occurrence probability of negative wave slope, overpredicted longer and shorter waves, and underprediction of waves of the medium period.

Non-Gaussian Wave Slope Distribution
Upon noting the limitations that Gaussian distribution presents, some efforts were made in this study to develop an alternative for Gaussian distribution preferred in Europe or the United States. As a result, the non-Gaussian distribution of wave slope is analytically derived from the joint distribution of wave amplitude and period by Longuet-Higgins [16] using the standard transformation technique of random variables [19]. In doing so, what improvements non-Gaussian wave slope distribution could make over the Gaussian distribution will also be examined.
The derivation of non-Gaussian wave slope distribution begins with the joint distribution of amplitude A and its associated period T by Longuet-Higgins [16], which can be written as: where ν is the bandwidth parameter and is defined as follows [16]: In Equation (20), m i denotes the i th moment of the wave spectrum S ςς (ω) and can be written as: (21) where ω is the angular frequency. After carrying out the deep-water approximation, wave slope can be written as: where L o is the deep water wave length. Non-Gaussian wave slope distribution can be derived from Equation (22) using the standard random variables transformation technique [19]. However, as of now, only one relationship is available, and as a result, an auxiliary random variable should be introduced to solve these ill-posed problems, which can be written as: The derivation proceeds by first obtaining the joint distribution of S om and S D , f S om S D (S om , S D ), from Equations (22) and (23), and the joint distribution of amplitude and its associated period f A T (A, T) defined in Equation (19) using the standard technique of transformation of random variables [19].
The transformation mentioned above can be written as: In Equation (24), Jacobian J|·| denotes the ratio of the area involved in mapping Lastly, the wave slope distribution can be obtained by carrying out the integration of f S om S D (S om , S D ) over the entire sample space of auxiliary random variable introduced to solve the ill-posed problem (marginal integration) and f S om (S om ), obtained in this way, can be written as: Figure 6 shows f S om (S om ) for ν = 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8. For the sake of comparison, Gaussian distribution with the same mean and standard deviation as non-Gaussian wave slope distribution presented in this study is included in Figure 7. In Figure 6, the enhanced bandwidth of the wave spectrum implies that sea wave conditions become harsh, and as a result, the number of component waves appearing in the random wave field increases. As resonant wave-wave interaction is initiated from a resonant triad formed among the component waves, the wave energy shift to longer and shorter waves, and these wave energy redistributions are intensified as the sea conditions become harsher [ ν → 1.0 ], which is well known in the coastal engineering community after the studies of Hasselmann [29] and Phillips [30] on the development of wind waves. Interestingly, as the sea conditions become rough [ ν → 1.0 ], the probability mass moves toward relatively steep and mild slope waves in Figure 6, which complies with the wave energy redistributions under resonant wave-wave interaction mentioned above. On the contrary, in Gaussian distribution, some anomalies such as the non-negligible occurrence probability of negative wave slope and overshooting near the mean wave slope can be found. In non-Gaussian wave slope distribution, when wave slope is relatively small, the probability increases rapidly from zero. After the rapid growth zone, it reaches the plateau, where the variation of wave slope distribution becomes mild. This trend becomes more pronounced with increasing ν (see Figures 6 and 7).
wave slope distribution presented in this study is included in Figure 7. In Figure 6, the enhanced bandwidth of the wave spectrum implies that sea wave conditions become harsh, and as a result, the number of component waves appearing in the random wave field increases. As resonant wave-wave interaction is initiated from a resonant triad formed among the component waves, the wave energy shift to longer and shorter waves, and these wave energy redistributions are intensified as the sea conditions become harsher [ 1.0 ν → ], which is well known in the coastal engineering community after the studies of Hasselmann [29] and Phillips [30] on the development of wind waves. Interestingly, as the sea conditions become rough [ 1.0 ν → ], the probability mass moves toward relatively steep and mild slope waves in Figure 6, which complies with the wave energy redistributions under resonant wave-wave interaction mentioned above. On the contrary, in Gaussian distribution, some anomalies such as the non-negligible occurrence probability of negative wave slope and overshooting near the mean wave slope can be found. In non-Gaussian wave slope distribution, when wave slope is relatively small, the probability increases rapidly from zero. After the rapid growth zone, it reaches the plateau, where the variation of wave slope distribution becomes mild. This trend becomes more pronounced with increasing ν (see Figures 6 and 7).   These unique features of non-Gaussian wave slope distribution are concerned with long waves appearing in the wave field due to the nonlinear sub-harmonic resonant interaction. The presence of probability plateau can be foreseeable from the joint distribution of amplitude and period by Longuet-Higgins [16] shown in Figure 8. The increase in ν implies that nonlinearity is more pronounced, and as a result, long waves appear more frequently in the wave field due to the resonant wave-wave interaction (see Figure 9), and the probability of a relatively small wave slope increases, leading to the formation of probability plateau in the wave slope distribution mentioned above. In light of the Van der Meer equation, it can be easily perceived that the deformation of wave slope distribution enforced by the resonance interaction significantly affects the armor block's hydraulic stability strength and its failure probability. Therefore, the non-Gaussian wave slope distribution presented in this study, which can account for the nonlinear sub-harmonic reso- These unique features of non-Gaussian wave slope distribution are concerned with long waves appearing in the wave field due to the nonlinear sub-harmonic resonant interaction. The presence of probability plateau can be foreseeable from the joint distribution of amplitude and period by Longuet-Higgins [16] shown in Figure 8. The increase in ν implies that nonlinearity is more pronounced, and as a result, long waves appear more frequently in the wave field due to the resonant wave-wave interaction (see Figure 9), and the probability of a relatively small wave slope increases, leading to the formation of probability plateau in the wave slope distribution mentioned above. In light of the Van der Meer equation, it can be easily perceived that the deformation of wave slope distribution enforced by the resonance interaction significantly affects the armor block's hydraulic stability strength and its failure probability. Therefore, the non-Gaussian wave slope distribution presented in this study, which can account for the nonlinear sub-harmonic resonant interaction and its effects on the wave slope distribution, is expected to play a significant role in the reliability design practice of armor block. long waves appearing in the wave field due to the nonlinear sub-harmonic resonant interaction. The presence of probability plateau can be foreseeable from the joint distribution of amplitude and period by Longuet-Higgins [16] shown in Figure 8. The increase in ν implies that nonlinearity is more pronounced, and as a result, long waves appear more frequently in the wave field due to the resonant wave-wave interaction (see Figure 9), and the probability of a relatively small wave slope increases, leading to the formation of probability plateau in the wave slope distribution mentioned above. In light of the Van der Meer equation, it can be easily perceived that the deformation of wave slope distribution enforced by the resonance interaction significantly affects the armor block's hydraulic stability strength and its failure probability. Therefore, the non-Gaussian wave slope distribution presented in this study, which can account for the nonlinear sub-harmonic resonant interaction and its effects on the wave slope distribution, is expected to play a significant role in the reliability design practice of armor block. Putting together the facts mentioned above, PIANC's [18] recommendation that wave slope follows Gaussian distribution should be revised. Moreover, armor block reliability design based on the non-Gaussian wave slope distribution presented in this study  Putting together the facts mentioned above, PIANC's [18] recommendation that wave slope follows Gaussian distribution should be revised. Moreover, armor block reliability design based on the non-Gaussian wave slope distribution presented in this study could effectively address too conservatively designed armor block problems frequently occurred in developing the countermeasures against rough seas due to climate change.

Verification of Non-Gaussian Wave Slope Distribution
In Figure 10, the times series data of measured peak wave period T P , significant wave height H S , and the wave slope S om between 26 April 2017-20 April 2018 using an ultrasonic wave height meter (Nortek) near the Mang-Bang beach [129 • 13 34.56 E, 37 • 24 11.22 N] (see Figure 1) on the east coast of Korea were depicted. Figure 11 showed the frequency analysis results of the situ wave slope data, and for the sake of comparison, Gaussian distribution of the same mean and standard deviation was also included. As discussed in Section 3.2.2, the Gaussian distribution limitation, such as providing a nonnegligible occurrence probability even on the negative wave slope, can be re-affirmed. The Gaussian distribution flaws such as under-prediction of longer and shorter waves and over-prediction of the medium-period waves can also be found.

Level III Reliability-Based Design of an Armor Block of Rubble Mound Breakwater
Based on the Van der Meer equation, the reliability function Z can be written as: With the given statistical properties of the reliability function Z, the reliability of the structure defined as the probability of Z > 0 can be written as: Even though the reliability function Z depends on six mutually independent random processes such as the N od , N z , ∆, D n , H S and S om ; hereafter, the reliability function Z is assumed to depend on the armor block's hydraulic stability strength R and wave force S acting on the structure for the sake of convenience in the following analysis of armor block's failure probability.
In this case, the typical behavior characteristics of the reliability function Z in the [R, S] plane can be depicted such as in Figure 12. In Figure 12, the joint distribution of R and S is also included for comparison. In Figure 12, f R (r) and f S (s) denote the probability density function of the armor block's hydraulic stability strength and wave force, respectively. Red dotted line denotes the limiting state, and black solid line denotes the linear approximation of the limiting state Z = 0.
When Level III reliability-based design is applied, the failure probability is analytically derived by carrying out the integration of the probability density function of Z over the sample space where Z < 0, which can be written as: In this case, the typical behavior characteristics of the reliability function Z in the [ , ] R S plane can be depicted such as in Figure 12. In Figure 12, the joint distribution of R and S is also included for comparison. In Figure 12, When Level III reliability-based design is applied, the failure probability is analytically derived by carrying out the integration of the probability density function of Z over the sample space where 0 < Z , which can be written as: Assuming that the armor block's hydraulic stability strength and wave force acting on the structure are mutually independent, the joint distribution of R and S, f RS (R, S), can be written as: Further assuming armor block's hydraulic stability strength and wave force to depend on the wave slope S om and H S , R can be written as: where A is given by: As the probabilistic model of wave force S = H T S , Weibull distribution is preferred in Europe and the United States, which can be written as [27]: where H W is the scale factor, and κ W is the shape factor. Assuming that wave force follows the Weibull distribution, the mean E H T S and variance E H T S 2 of wave force can be written as: In Equations (34) and (35), E[·] denotes the expected value of the random variable in the bracket, and Γ[·] is the Gamma function [31].
Applying the standard transformation technique of random variables [19] to Equation (31), the probability density function of armor block's hydraulic stability strength R can be written as: where Jacobian J|·| denotes the ratio of the area involved in mapping [S om ] → [R] . After completing the transformation involved in Equation (36); finally, the probability density function of R can be written as: From Equations (29), (30), (33), and (37), the failure probability of armor block P[R < S] can be written as: Here, upon assuming that the armor block's hydraulic stability strength and wave force are mutually independent, Equation (38) can be rewritten as: After substituting the Weibull distribution into Equation (39) and carrying out the integration, the following relationship can be obtained: After inserting Equation (40) into Equation (39) and carrying out the integration (see Figure 13); lastly, the failure probability of armor block P[R < S] can be written as: where f (ς) is given by: Even though the armor block's failure probability P[R < S] in Equation (41) cannot be given in closed form, using the stationary phase method [32] based on the fact that the integrand in Equation (41) is concentrated near the peak ς = ς o , P[R < S] in Equation (41) can be written as: where f (ς o ) is given by: Stationary point ς o in Equations (43) and (44) can be evaluated from the following relationship:

Level III Reliability-Based Design of an Armor Block Using the Probabilistic Model of Wave Height Optimized for the Korean Sea Wave Conditions and Gaussian Wave Slope Distribution
The modified Glukhovskiy distribution by Choi and Cho [1] can be rewritten as follows:

Level III Reliability-Based Design of an Armor Block Using the Probabilistic Model of Wave Height Optimized for the Korean Sea Wave Conditions and Gaussian Wave Slope Distribution
The modified Glukhovskiy distribution by Choi and Cho [1] can be rewritten as follows: Therefore, it can be easily perceived that probability parameter in Weibull distribution can be written as To demonstrate the robustness of the probabilistic model presented in this study, the probability density functions of the armor block's hydraulic stability strength for varying D n [2.0, 2.05, · · · , 2.5 m] were depicted in Figure 14. It was shown that with the increase in D n , a significant probability mass shifts to the higher hydraulic stability strength of armor block as expected, and in doing so, the standard deviation also increases. As D n is increasing by 0.05 m, the armor block's hydraulic stability strength increases by as much as 2.5%. In Figure 15, the joint distribution of R and S was depicted for D n = 1.9 m and 2.3 m, respectively. It was shown that as D n is decreasing, a significant probability mass shifts to the weaker hydraulic stability strength of armor block. n D , a significant probability mass shifts to the higher hydraulic stability strength of armor block as expected, and in doing so, the standard deviation also increases. As n D is increasing by 0.05 m, the armor block's hydraulic stability strength increases by as much as 2.5%. In Figure 15, the joint distribution of R and S was depicted for 1 .9 n D m = and 2.3 m , respectively. It was shown that as n D is decreasing, a significant probability mass shifts to the weaker hydraulic stability strength of armor block.   Figure 16 shows armor block failure probability as n D varies, and the failure probability evaluated based on the Weibull distribution is also included for comparison. In the Weibull distribution, the armor block's failure probability is overestimated compared to the probabilistic wave height model optimized for the Korean marine environment, which means that the armor block has been designed too conservatively. In Tables 1 and 2 Figure 16 shows armor block failure probability as D n varies, and the failure probability evaluated based on the Weibull distribution is also included for comparison. In the Weibull distribution, the armor block's failure probability is overestimated compared to the probabilistic wave height model optimized for the Korean marine environment, which means that the armor block has been designed too conservatively. In Tables 1 and 2 Table 2).   Figure 16. Failure probability of TTP as D n varies for different wave height distribution (see Table 2).  In this section, the effect of the non-Gaussian wave slope distribution (see Section 3.2.2) on the joint distribution of armor block's hydraulic stability strength R and wave force S, the probability density function of armor block's hydraulic stability strength, and armor block's failure probability will be examined. Figure 17 shows the joint distribution of the armor block's hydraulic stability strength and wave force for D n = 2.3 m and ν = 0.4, and the joint distribution of the armor block's hydraulic stability strength and wave force based on the Gaussian wave slope distribution was also included for comparison. The most discernable difference is that in the case of non-Gaussian wave slope distribution, the center of probability mass shifts toward the weaker hydraulic stability strength of armor block [3.1 < R < 3.4], and in doing so, the sample space in the [R − S] plane where non-negligible probability mass resides shrinks as well. Figure 18 shows the variation of the probability density function of the armor block's hydraulic stability based on the non-Gaussian wave slope distribution as wave conditions become harsher. force S , the probability density function of armor block's hydraulic stability strength, and armor block's failure probability will be examined. Figure 17 shows the joint distribution of the armor block's hydraulic stability strength and wave force for , and the joint distribution of the armor block's hydraulic stability strength and wave force based on the Gaussian wave slope distribution was also included for comparison. The most discernable difference is that in the case of non-Gaussian wave slope distribution, the center of probability mass shifts toward the weaker hydraulic stability plane where non-negligible probability mass resides shrinks as well. Figure 18 shows the variation of the probability density function of the armor block's hydraulic stability based on the non-Gaussian wave slope distribution as wave conditions become harsher.
(b) Gaussian.   It is worthy of note that as sea wave conditions are becoming harsher, the armor block's weakened hydraulic stability strength is quantitatively verified. Moreover, as sea wave conditions are becoming harsher [ 1.0 ν → ], the probability density mass shifts toward stronger and weaker hydraulic stability strength of armor block, which also complies with the deformation characteristics of the non-Gaussian wave slope distribution enforced by the nonlinear resonant wave-wave interaction (see Section 3.2.2).
In Table 3, the most probable armor block strength for varying bandwidth of wave spectrum was summarized. In Figure 19, the probability density function of R based on the non-Gaussian wave slope distribution is compared with the one based on the Gaussian distribution for   It is worthy of note that as sea wave conditions are becoming harsher, the armor block's weakened hydraulic stability strength is quantitatively verified. Moreover, as sea wave conditions are becoming harsher [ ν → 1.0 ], the probability density mass shifts toward stronger and weaker hydraulic stability strength of armor block, which also complies with the deformation characteristics of the non-Gaussian wave slope distribution enforced by the nonlinear resonant wave-wave interaction (see Section 3.2.2).
In Table 3, the most probable armor block strength for varying bandwidth of wave spectrum was summarized. In Figure 19, the probability density function of R based on the non-Gaussian wave slope distribution is compared with the one based on the Gaussian distribution for D n = 2.3 m, and ν = 0.6. It was shown that while armor block's hydraulic stability strength R should have a larger value [3.8 < R < 5.0] in the mild seas where moderate swells are prevailing, in Gaussian wave slope distribution, mild seas are occurring less frequently, leading to the under-estimation of armor block's hydraulic stability strength.  In order to examine how the non-Gaussian wave slope distribution affects the failure probability of armor block P[R < S], in Figure 20, P[R < S] as D n varies was depicted, and P[R < S] based on the Gaussian wave slope distribution was also included for comparison. In the Gaussian wave slope distribution, the armor block's failure probability was overestimated, and these trends are enhanced when D n is small. On the other hand, the over-estimation problem is diminished as D n is becomes larger.

Conclusions
The armor block of rubble mound breakwater built in a rough sea should be designed to be robust enough to survive even in harsh wave conditions. However, most armor block of a rubble mound breakwater in Korea has been designed using a deterministic design method based on Hudson equation [2]. Lately, in developing countermeasures against rough seas due to climate change, the Hudson equation's stability coefficient was readjusted to have smaller values, and as a result, it was not difficult to find the enormous armor blocks such as a TTP of 100 t. These design practices have raised a critique in the coastal engineering community that the port's outer facilities in Korea are too conservatively designed. An alternative that can address the limitations that the deterministic design method presented seems to be reliability design. If the irregularities inherent in the marine environment are considered, reliability design should be carried out based on the long-term wave data. However, situ wave data continuously measured over a long period are very rare in Korea. Therefore, reliability design was performed based on the probability model preferred in the United States [3][4][5].
In this rationale, this study developed a reliability design method of armor block of rubble mound breakwaters based on the Van der Meer equation and a probabilistic wave height model optimized for the Korean marine environment. In this process, the armor block's failure probability is analytically derived from the Van der Meer equation using the random variable transformation technique [19]. Now that the armor block's failure probability is analytically derived, this study's design method could be classified as a Level III reliability design. The design method presented in this study does not depend on the specific probability distribution of the random variables involved in the van der Meer equation and can be extended to any kinds of probability distributions available in the literature [18].
To demonstrate what distinguishes this study from the others, numerical simulation was first carried out, assuming that the wave slope follows the Gaussian distribution recommended by PIANC [18]. In doing so, the modified Glukhovskiy distribution by Choi and Cho [1] was used as a probabilistic wave height model. It was shown that Gaussian wave slope distribution has flaws such as the over-estimated failure probability of armor block and the non-negligible occurrence probability of the negative wave slope. It was also shown that Gaussian wave slope distribution over-predicted longer and shorter waves and under-predicted waves of the medium period. Upon noting the limitation of Gaussian distribution, some efforts were also made to develop an alternative for Gaussian distribution. As a result, the non-Gaussian wave slope distribution was analytically derived from the joint distribution of wave amplitude and period by Longuet-Higgins [16] using the random variable transformation technique.
Numerical results showed that non-Gaussian distribution could effectively address the limitations of Gaussian distribution due to its capability to account for the nonlinear resonant wave-wave interaction and its effects on the wave slope distribution. In the case of Gaussian wave slope distribution, the armor block's failure probability was overestimated, and these trends were enhanced when the armor block is small. On the other hand, these over-estimation problems were diminished as the armor block became larger. These discrepancies occurred since the armor block should be more robust and of immense hydraulic stability strength in the mild seas where mild slope waves prevail. However, Gaussian wave slope distribution predicted that mild slope waves occurred less frequently in the mild seas due to its lack of apparatus accounting for the nonlinear sub-harmonic resonant wave-wave interaction.
In light of the Van der Meer equation, it can be easily perceived that the deformation of wave slope distribution enforced by the presence of long waves appeared in wave field due to the nonlinear sub-harmonic wave-wave interaction has a significant effect on the armor block's hydraulic stability strength and its failure probability. Therefore, the non-Gaussian wave slope distribution presented in this study could play an indispensable role in resolving the controversial issue such as whether or not enormous armor blocks like a Tetrapod of 100 t frequently mentioned in developing the countermeasures against the rough seas due to climate change is too conservatively designed.
However, it is worth mentioning that even though the Level III reliability design method of an armor block presented in this study is not limited to a specific site, to extend its application to other sites such as the west coast of Korea, it is a prerequisite that the probability parameters of modified Glukhovskiy wave height distribution should be re-evaluated using the hindcasted wave data for at least 15 years.