Damage in Rubble Mound Breakwaters. Part I: Historical Review of Damage Models

Abstract

The term “damage” in rubble mound breakwaters is usually related to the foremost failure mode of this kind of coastal structures: their hydraulic instability. The characterization of the breakwater response against wave action was and will be the goal of hundreds of studies. Because of the large amount of information, the present review on damage in rubble mound breakwaters is divided in two papers, which are closely linked but conceptually different; whereas Part II is focused on the various approaches for defining and measuring damage, Part I summarizes the diverse strategies for modelling damage development and progression. The present paper compiles 146 references on this topic, chronologically discussed over almost a century of history: from 1933 to 2020. It includes 23 formulations of hydraulic stability models and 11 formulations of damage progression models, together with main advances and shortcomings up to date. The future of rubble mound design is linked to risk-based tools and advanced management strategies, for which deeper comprehension about the spatial and temporal evolution of damage during the useful life of each particular structure is required. For this aim, damage progression probabilistic models, full-scale monitoring and standardization will presumably be some of the key challenges in the upcoming years.

1. Introduction

Damage in rubble mound breakwaters can be defined as the partial or total loss of its functionality and it is usually related to the hydraulic instability of the armor layer [1,2]. Rubble mound breakwaters behave as a static granular system until wave energy is higher than the one needed for the initiation of movement. At that moment, wave forces might provoke movements in the armor layer such as rocking, displacements, sliding of a whole blanket or settlement because of compaction. In this way, pieces close to SWL (still water level) tend to be moved toward the toe of the structure. This causes the geometrical initial shape to evolve with wave attack, up to the point of dealing with erosion and accretion volumes that can be especially high in the case of berm breakwaters. When underlayer units start to be removed, the progression of the geometrical rearrangement is much quicker, as the energy needed to move lighter units is obviously lower. Once the core material is exposed, the destruction of the structure is imminent.

Hydraulic instability of armor layers is a complex process because of the stochastic nature of both wave loading, initiation of movement, and damage progression. The highly non-linear flow over the slope, involving wave breaking, together with the variable shape of armor units and their random placement, make unfeasible to achieve analytical expressions for the calculation of actions and reactions. Thus, determining instantaneous armor unit stability also seems unachievable. This is the reason why stability formulae are historically empirical or semi-empirical, associated with experimental results and ideally improved with prototype observations. Indeed, as stated in Melby and Hughes [3], it is recommended to utilize physical models if at all possible to verify stability and it is a common requisite in design tenders. These physical models are aimed at relating the response of the armor layer to parameters of the incident wave train and breakwater characteristics.

Since the pioneering works of de Castro [4] and Iribarren [5], several stability models have been proposed and many hundred studies on breakwater armor layer stability have been published. This gives an idea of the relevance of this failure mode; however, the achievement of a unique formula seems to be rather difficult, partly because of the wide range of typologies and designs. Regarding the large amount of data on this topic, this paper is aimed to present a historical overview of the most relevant models and advances up to date, focused on statically stable rubble mound breakwater’s trunks. Breakwater’s heads and other rubble mound typologies such as berm breakwaters are not considered herein.

The review on damage in rubble mounds breakwaters is divided into two separated papers, which are closely linked but conceptually different. Part I is focused on how damage development and progression is modelled, i.e., which are the available hypothesis and strategies for characterizing the structural response. This characterization is especially relevant, not only from the point of view of the breakwater design, but also regarding conservation and maintenance plans during the useful life. Part II [6] is centered on how damage is measured and parametrized, i.e., which are the instrumental techniques available for monitoring the structural response and, based on these instrumental records, which are the strategies for defining damage descriptors. Note that damage progression models presented in Part I directly consider these descriptors within their formulations.

2. Historical Review of Damage Models for Rubble Mound Breakwaters

Most hydraulic stability formulae were derived for rock armoring rather than for artificial blocks. According to Allsop et al. [7], rock armored breakwaters dominate in many areas of the world, although concrete armored breakwaters have probably a more detailed database because of the records from armor unit licensees. The present manuscript is mainly focused on quarry stones and parallelepiped armor units. Special shaped concrete units are roughly outlined in Appendix A.

In Table 1 most hydraulic stability formulae are summarized, with some slight modifications in order to be consistent with the symbols used herein and to fit the following general structure (adapted from Hald [8]):

$$N_S=\frac{H}{\mathsf{\Delta }{D}_{n50}}<f(K,p_1,p_2,\dots ,p_n)$$

(1)

Thus, the stability of the armor layer was found to be reached when the stability number (N_S) is lower than a certain function (f). This function depends on the n parameters (p₁, p₂, …, p_n) influencing stability and an empirical coefficient (K) determined by the parameters not directly accounted for in the stability equation (see Appendix B for the rest of symbols).

While hydraulic stability formulae were originally aimed at characterizing the initiation of movement of the armor layer, damage progression models were later designed to predict the evolution of rubble mounds’ geometry by means of a quantitative damage descriptor. Consequently, damage progression models represent a much more useful design tool and are indispensable in any reliable maintenance or conservation program. However, because of the complexity of the problem, damage progression models were developed just in the past decades. More research is still needed, especially considering the adequacy of the different models to prototype measurements, which is not extensive nowadays. A review of damage progression models is presented in Table 2. Note that some models, mainly based on experimental results, define damage through the dimensionless erosion area (S) presented in Broderick and Ahrens [9], while others, mainly built under theoretical assumptions, define the damage descriptor in a generic way (D).

Before expounding on a historical synthesis on armor layer stability from the first hydraulic stability models to the latest contributions, the timeline of Table 3 is presented. It shows a selection of the most relevant hydraulic stability models and damage progression models from Table 1 and Table 2. This scheme permits, not only arranging chronologically the different proposals in the study of armor layer stability, but also having an idea of the parameters/properties accounted for and the innovations introduced by each of them.

2.1. First Approaches to Hydraulic Stability Characterization. Models from Iribarren and Hudson

The first published formula for the calculation of the design weight of armor units evidenced by the authors was proposed in 1933. In this year, de Castro [4] reported a brief paper about rubble mound breakwaters using a modest approach based on his experience and highlighting the prevailing importance of rundown forces on instability.

The recommendations of de Castro inspired Iribarren [5] to develop another model in 1938, considered by many authors as the pioneer work on armor stability. In 1950, after 12 years of analyzing the consistency of the previous formula, Iribarren and Nogales [10] ratified the empirical coefficients presented in 1938. They also generalized the formula by introducing some modifications in the wave height parameter in order to account for water depth and wave period. Their final model was published in 1965 [11], summarizing the previous work and studying the limitations on the application of the formula by analyzing different slopes and types of armor units and wave breaking conditions. Also, a simple stability curve was suggested, relating the percentage of damaged slope and a coefficient between wave height and the wave height causing total destruction of the structure.

The initial works of Iribarren were translated and published in English in 1949 [12] and, after that, several models were proposed all around the world. In France, Larras [13] developed in 1952 a model considering water depth and wave length and, three years later, Beaudevin [14] published a simpler one. In the United States, Mathews (unpublished report presented in Bruun [15]) and Epstein and Tyrrel [16] built the first North-American hydraulic stability formulations in 1948 and 1949 respectively. Hickson and Rodolf [17] proposed a model for jetties based on the mining experience in 1950 and Hudson and Jackson [18] presented a formula in 1953 consisting of some modifications on Iribarren’s. In Sweden, Hedar [19] ratified the Spaniard’s results in 1953 reaching to the same stability criterion during downrush despite considering rotation, and not sliding, to be the dominant mode of motion. In Norway, Svee [20] ended up in 1962 with a model considering just a lifting force normal to the slope. In the Soviet Union, there were many publications on this topic between 1959 and 1967 such as Goldschtein and Kononenko [21], the design code SN-92-60 [22], Rybtchevsky [23], or Metelicyna [24].

However, the most widely known empirical stability model corresponds to Hudson [25,26]. After trying to characterize the friction coefficient for different armor units, Hudson noticed that the experimental coefficients of Iribarren’s model could not be determined accurately from small scale breakwater stability experiments because they suffer variations from test to test for the same experimental conditions. Thus, Hudson proposed in 1958 a simpler formula after carrying out new stability tests with normal non-overtopped regular waves and using quarry stones and tetrapods (it was extended a few years later to tribars, tetrahedron and other special-shape armor units):

$$W=\frac{{\gamma }_s\cdot H^3}{K_D{\left(S_r-1\right)}^3\cot \left(\alpha \right)}$$

(2)

where K_D is tabulated for each armor unit type as a function of the H/H_D=0 ratio. H_D=0 is the limit wave height that produces no damage, considered as less than 1% volume of units eroded relative to the total volume of stones in the active armor layer (see Appendix B for the rest of symbols). At the same time, the H/H_D=0 ratio was related to a certain percentage of damage, which means that Hudson’s model was probably one of the firsts providing quantitative information about the level of damage. The latter was possible owing to a standardized method for damage profiling, developed by the U.S. Army Engineer Waterways Experiment Station (WES). Despite not considering water depth or wave period, Hudson recommended his model as an initial approximation of the major forces from both breaking and non-breaking waves. However, he pointed out that “there is some doubt as to which of the various wave heights in natural wave trains should be selected as the design wave”. In fact, the Shore Protection Manual (SPM) suggested first H_1/3 [27] and, afterwards, H_1/10 [28] as an alternative wave height parameter in Equation (2) regarding irregular waves, without offering a clear justification for this modification.

Probably one of the first laboratory tests on quarry stones stability using irregular waves corresponds to Carstens et al. [29]. The experiments were carried out in 1966 at the Technical University of Norway, using a magnetic tape for sending the electric signal to the piston paddle, and testing wide and narrow spectra.

The difficulty to select an appropriate design wave for the stability models was also pointed out by Font [30] in 1968, together with empirical evidence on the influence of storm duration on rubble mound breakwater stability.

In 1974, Battjes [31] introduced Iribarren’s number (ξ), also known as surf similarity parameter, in the study of smooth and impermeable slopes. This parameter was also found to be practical for characterizing the breaker type, run-up and run-down on both smooth and permeable slopes by Bruun and Günbak [32] in 1976.

In the 70s, there was an effort to compile the formulae available up to date and to provide some design recommendations. Bruun [33] summed up the most relevant stability formulae and presented some clues for rubble mounds design, such us increasing block size on breakwater’s heads. This was also stated by Iribarren in 1965 [12]. In addition, the first English speaking design manual, the already-referred SPM, was published in 1973, and re-edited in 1975, 1977, and 1984. Also, a complete list of the available hydraulic stability models (mainly developed using regular waves) was provided in the final report of PIANC [34] in 1976.

2.2. 70s to 80s: Intensive Research on the Stability of Rubble Mound Breakwaters. Van der Meer’s Formulae

At the end of the 70s and beginning of the 80s, catastrophic failures were experienced by a series of large rubble mound breakwaters, as it was reported by Oumeraci [35]. This shock to the profession had two effects. First, for the rehabilitation of some damaged rubble mounds the old berm breakwater concept was re-discovered. Second, there was an extensive increase on research activities toward improving the design and construction of rubble mound breakwaters. Because of the ever-growing dimensions of the structures and the need to move into more hostile environments, reliable design formulae were demanded by the coastal engineering community. Van Hijum [36] studied in 1976 the equilibrium profile of coarse material under wave attack, which was the basis for subsequent research on berm breakwaters. Bruun [37] analyzed in 1978 the common reasons for damage or breakdown of mound breakwaters, which were ratified from a practical point of view by the 25-years-experienced engineer Kjelstrup [38] in 1979. In addition, many experimental works on breakwater stability were carried out:

(1)

Ahrens [39] and Ahrens and McCartney [40] conducted in 1975 several large-scale tests using regular waves and stated, among others, the influence of wave period on riprap stability.

(2)

Thompson and Shuttler [41] presented a detailed study on riprap stability in 1975. They concluded that “the erosion damage caused by irregular waves on a riprap slope is itself a random variable” and that the method for positioning the stones highly affects damage evolution. Long-term and short-term experiments were tested with impermeable core, finding no relation with depth or wave period, but with mean number of zero crossing waves, among others. In addition, it was suggested that damage tends to reach equilibrium or, in other words, that damage curves are meant to be asymptotic.

(3)

Bruun and Günbak [32] assured in 1976 that, regarding breakwater stability, “the significance of wave period is clearly demonstrated.” They started to work on a risk-based approach in the design of sloping structures considering Iribarren’s number.

(4)

Whillock and Price [42] presumably coined the concept of “fragility” in 1976 when indicating that elements with high void ratio and designed to interlock, such as dolosse, were associated with a reduction in the safety margin as failure is approached. They also supported the idea of quarry stones being more stable to oblique wave attack than to normal attack. However, this assumption was denied for armor units that are susceptible to drag forces. Indeed, they demonstrated that the overall stability for a particular test with dolosse decreased up to an angle of incidence of 60°. In line with the concept of “fragility”, Magoon and Baird [43] underlined in 1977 the importance of rocking movements in the breakage of armor units designed to interlock.

(5)

Losada and Giménez-Curto [44] proposed in 1979 a stability model for the initiation of damage by means of design curves depending on armor unit type and Iribarren’s number. They used Iribarren’s data to fit their model, and compared it satisfactorily with the results from Hudson [26] and Ahrens and McCartney [40]. However, they stated the difficulty to establish a comparison between experimental results undertaken in different laboratories because of divergences in both experimental process and damage criteria. For this reason, it was remarked that “To obtain general criteria on the behaviour of rubble mound breakwaters under wave action, laboratories should establish uniform experimental procedures and criteria.”

(6)

Losada and Giménez-Curto [45] studied the influence of oblique incidence in 1982. They concluded that, for gravity armor units, the stability of steep slopes under oblique wave attack is not worse than for normal incidence. On the contrary, they found that for high interlocking armor units (such as dolosse or tetrapods) oblique wave attack is hazardously worse than normal incidence, in agreement with Whillock and Price [42].

(7)

Broderick and Ahrens [9] presented in 1982 a technical paper on scale effects using the large-scale tests from Ahrens [39]. These tests, with wave heights up to 1.83 m and periods up to 11.3 s were compared with a 1:10 Froude scaled model. They found a 20% reduction in the zero-damage stability numbers from the large-scale tests whereas run-up increased about 20%. This was identified to be due to the improper modeling of the flow regime within the filter layer for the small-scale experiments and to the lack of penetration of the wave run-up into the filter layer. Scale effects were less severe at high levels of damage and the shapes of the damaged profiles for tests with the same relative depth were very similar. Furthermore, at the zero-damage level, wave period had less influence on the small scale. In their study, the widely used dimensionless erosion area (S) was firstly proposed.

(8)

Jensen [46] published in 1984 a thorough monograph on rubble-mound breakwaters. He suggested the parameter H_13.6% instead of H_S as a descriptor for wave height in both deep waters (where wave height is usually characterized by a Rayleigh distribution) and shallow waters (where a certain percentage of waves break). Also, in order to reach a stable damage level, scaled storms representing at least 8 to 10 hours in prototype were recommended. Furthermore, possibly the first formula on rear slope stability was proposed.

(9)

Hedar [47] developed in 1986 a complete stability model taking into account the water depth, wave height at breaking, rubble mound slope, the internal friction angle, a permeability function, and two empirical coefficients.

This intensive research on breakwater stability evidenced many shortcomings of the simple and widely used Hudson’s formula. In this context, Van der Meer [48,49,50,51] referred as VdM in this manuscript, developed in 1985 a revolutionary stability model with a dissemination similar to Hudson’s. It considered the influence of wave height (H), wave period (by means of Iribarren’s number, ξ₀), number of waves (N_w), equivalent cube length (D_n50), relative excess specific weight (Δ), breakwater slope (cotα), and permeability of the core (P). Furthermore, the aforementioned dimensionless erosion area (S) was introduced in the formulae as follows:

$$\frac{H}{\mathsf{\Delta }{D}_{n50}}=6.2{\xi }_0^{-0.5}{P}^{0.18}{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}\hspace{1em}\to for plunging waves$$

(3)

$$\frac{H}{\mathsf{\Delta }{D}_{n50}}=1.0{\xi }_0^P\sqrt{\cot \alpha } P^{-0.13}{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2} \to for surging waves$$

(4)

The formulation distinguished between plunging waves and surging waves, and even included the possibility of being applied in a probabilistic design. For this purpose, the coefficient 6.2 in Equation (3) can be assumed to be normally distributed with a standard deviation of 0.4. Similarly, the coefficient 1.0 in Equation (4) can be assumed to be normally distributed with a standard deviation of 0.08. In the formulation, minimum stability is reached when ξ ≈ 3. This situation corresponds to collapsing waves: while plunging waves are meant to cause damage during run-up and surging waves during run-down, collapsing waves are supposed to erode armor layer both during run-up and run-down.

VdM performed 262 tests and also considered the 300 tests from Thompson and Shuttler [41] finding, contrary to them, a damage relation with wave period after reanalyzing their data. He tested impermeable, conventional, and homogeneous rubble mound models concluding that stability is straightly related to permeability. Despite using riprap and selected natural rocks, he found no influence of armor grading on stability. This is why just the parameter D_n50 was considered in the formulation. Moreover, after trying two different Pierson-Moskowitz spectra (a wideband and a narrowband), no relation was found between damage and either spectral shape or wave groupness. Strong differences between monochromatic and irregular tests were pointed out, together with a damage dependence on the number of waves. Indeed, as it was reported by different authors, these formulae are only valid for 1000 to 7000 waves and it tends to overestimate the damage for more than 8000 waves. In 1998, VdM [51] related the Hudson’s “no-damage” criteria and filter exposure (failure criteria) to different values of S depending on breakwater slope. Also, an adaptation of his single storm formulae to calculate the damage caused by more than one storm event was proposed.

VdM [52] also developed in 1988 a formulation for cubes, tetrapods and accropode, singling out the study for the most common slope for each armor type. A different damage descriptor was used in this case, firstly introduced by Hedar [53] in 1960: N_o, the number of units displaced out of the armor layer within a strip width of one equivalent cube length. This equivalent cube length was assumed to be equal to D_n50 for cubes, 0.65 h for tetrapods and 0.7 h for accropode.

$$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=\left(6.7\frac{N_o^{0.4}}{N_w^{0.3}}+1\right){\left(\frac{H_S}{L_0}\right)}^{-0.1}\hspace{1em}for cubes$$

(5)

$$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=\left(3.75\frac{N_o^{0.5}}{N_w^{0.25}}+0.85\right){\left(\frac{H_S}{L_0}\right)}^{-0.2}\hspace{1em}for tetrapods$$

(6)

$$Failure, N_0>0.5\to \frac{H_S}{\mathsf{\Delta }{D}_{n50}}=4.1 for accropode$$

(7)

2.3. Studies Based on VdM’s. Van Gent’s Formula

The results from VdM were analyzed by many authors and further conclusions and improvements on armor stability were settled down. Among the extensive research based on VdM approaches, the following studies are highlighted:

(1)

Using exactly the same rocks tested by VdM, Latham et al. [54] demonstrated in 1988 the dependency of damage regarding armor shape. An additional coefficient to consider the armor shape effect was suggested: as rock units become more rounded they also become more unstable, especially under surging waves.

(2)

Medina and McDougal [55] highlighted in 1990 some shortcomings of VdM formulae, such as the complexity of the equation or the overestimation of damage for more than 7000 waves. They were especially critical with the independence of VdM’s model on wave groupness. In fact, a small but systematic higher stability for random waves from narrowband spectra was detected in VdM tests. A simpler stability model was alternatively proposed.

(3)

Kaku et al. [56] found in 1991 that VdM’s damage levels were not accurate enough for forecasting models. Therefore, they proposed a new empirical model assuming that the damage level approaches to an asymptotic equilibrium for the same wave energy. They also included the initiation of armor movement based on the similarity between the stability number and the Shields parameter used for sediment transport. However, in 1992, Smith et al. [57] indicated the difficulties of both empirical formulae in accounting for the complex wave and structural interactions affecting breakwater reshaping, mainly because they were built up under static stability conditions.

(4)

Melby and Hughes [3] utilized in 2003 part of the small-scale laboratory data of VdM to fit a stability equation derived from basic principles for uplift, sliding, and rolling incipient motion. It was based on the assumption that the maximum wave momentum flux at the toe of the structure is proportional to the maximum wave forces on armor units.

(5)

Vidal et al. [58,59] proposed in 2004 a modification of VdM’s formulae after accomplishing a comparison between the results from Thompson and Shuttler [41], Losada and Giménez-Curto [44], and VdM [50]. This modification consists of using H₅₀ instead of H_S. H₅₀ is defined as the average wave height of the 50 highest waves reaching the rubble-mound breakwater. This new way of describing wave height parameter is further discussed in Section 2.4.

(6)

Mertens [60] attempted in 2007 to transform the datasets of Thompson and Shuttler [41] and VdM [50] into comparable information with the one generated by Van Gent et al. [61] in 2003. He also reported some deviations in VdM data because of the influence of stone roundness. Following this line, Verhagen and Mertens [62] proposed in 2009 a methodology for unifying the formulation for both deep and shallow waters. To accomplish this, a correction factor based on the Iribarren’s number was added in order to incorporate the effect of the foreshore. For the correct application of this method they claimed for an accurate calculation of the wave height and wave period at the toe of the breakwater, including a precise determination of H_2% and T_m−1,0.

Among these works, probably the most renowned corresponds to Van Gent et al. [61], who reanalyzed VdM’s results and added some new experimental data in shallow water conditions. Initially, the parameter H_2%/H_S was added. However, after noticing that the influence of this parameter was small, as so was the influence of wave period regarding the data scatter due to other reasons, a single and simpler formula was proposed. In the new formula the permeability was incorporated in a direct way by the nominal diameter of the core material:

$$\frac{H}{\mathsf{\Delta }{D}_{n50}}=1.75{\left(\frac{S}{\sqrt{N}}\right)}^{0.2}\sqrt{\cot \alpha }\left(1+\frac{D_{n50,core}}{D_{n50}}\right)$$

(8)

2.4. Wave Groupness and Wave Height Parameter Discussion. H₅₀ Parameter by Vidal et al.

The correlation between wave groups and armor layer damage was a controversial topic for the researching community especially during the 80s. Some authors suggested that wave groups with short wave runs were more damaging, others stated that long wave runs were more damaging and others, such as VdM, did not find any major damage dependence on wave grouping. Medina et al. [63,64] recommended in 1990 avoiding the use of traditional wave grouping parameters, such as spectral peakedness or mean run length, after not finding significant correlation between them and armor damage. Alternatively, using linear theory (and therefore excluding breaking design conditions), they proposed an envelope exceedance coefficient accounting for the variability of the energy flux, related to a newly defined groupness factor. Both coefficients were proved to be capable of explaining mean damage variability higher than 50% for a given design sea state and storm duration.

The new parameters for taking into account wave groupness in breakwater damage were found to be related to the maximum wave heights of the design sea state by Medina et al. [63,64]. This was the trigger for restarting a debate on finding an adequate wave height descriptor, i.e., a representative parameter of the spectrum/sea state regarding breakwater’s damage. Wave height is raised to the third power in most stability models and, thus, damage is extremely sensitive to the way it is defined. As it was mentioned before, Hudson [26] pointed out the uncertainty of selecting a “design wave” from a natural wave train. Also, the SPM did not offer a clear justification for replacing the design wave from Hudson, obtained using regular waves, by the equivalent H_1/3 (in the third edition [27]) and H_1/10 (in the fourth edition [28]) for irregular waves. Indeed, for some authors H_1/10 is extremely conservative, and for others, such as Jensen et al. [65], H_1/20 is reported to be a more suitable parameter. Despite H_S is, even nowadays, the most extended parameter for the characterization of irregular waves, it does not give enough information about the highest waves of the series which, in fact, are the most likely to be responsible for armor damage development. Other studies discussing the adequacy of alternative wave height descriptors are highlighted below:

(1)

Teisson [66] presented in 1990 a statistical approach for characterizing the duration of extreme storms and its consequences on breakwater damage. In the study it was stated that “to select H_S as design wave height in Hudson formulae assumes that the associated storm will last for only 10 minutes: this choice could lead to an under estimation of breakwater design.” Alternatively, the expression H_D = 1.18 H_S t^0.095 was proposed for the calculation of the design wave height. Teisson tried to relate regular and irregular wave effects on stability. Furthermore, by an integrated theoretical approach, he developed a step-by-step methodology for calculating cumulated damage during storms, assuming that a storm can be described by a sequence of significant wave heights steps with a certain duration each.

(2)

Vidal et al. [67] suggested in 1995 that the wave height parameter should contain information about the distribution of the highest waves, the length of the time series and the number of times it is recycled to achieve a given degree of damage. The longer the test is, the higher waves are likely to attack the structure, i.e., damage evolution after testing two time series with a certain H_S and a duration t will be different compared to testing a unique time series with the same H_S and duration of 2t. Based on numerical simulation, they reported variations in the damage parameter exceeding 50% when considering H_S after testing JONSWAP spectra with the same H_S and T_p but different random seeds. They also confirmed that the highest waves were related to wave groupness. An H_n parameter, directly related to test duration and suitable also for breaking conditions (i.e., where wave height distribution during storms can depart from Rayleigh’s due to non-linearity), was proposed for a better characterization of wave-damaging energy. This new parameter was expected to facilitate the comparison of stability results obtained in different investigations, including those carried out with regular waves. Initially, H_n=H₁₀₀ was proposed, but in Vidal et al. [58,59] this parameter was adjusted to H₅₀ after comparing the datasets from Thompson and Shuttler [41], Losada and Giménez-Curto [44] and VdM [50]. In addition, a new formula was developed based on VdM equations. The lack of consistency of the Rayleigh distribution for breaking conditions was also accounted for by other authors, such as Battjes and Groenendijk [68], which suggested a Weibull distribution for damage models in shallow foreshores after a spectral analysis, or Méndez and Castanedo [69], which provided a model for the depth-limited distribution of the highest waves in a sea state.

(3)

Jensen et al. [65] tested both regular and irregular waves for identifying a wave height parameter within the irregular waves corresponding to the wave height of a regular series with a similar damage level on the structure. In line with Vidal et al. [67], they ended up in 1996 with an H_n parameter, but in this case a n-value of approximately 250 was found to be more suitable.

(4)

Medina [70] claimed in 1996 for non-stationary stochastic models as more adequate for modeling real waves. He defined five conditions for any rational armor damage model to properly take into account the storm duration, such as damage must necessarily increase with the duration of the storm under random wave attack in deep water conditions. A wave-to-wave exponential model accomplishing these conditions was proposed. The model depends on the number of waves and introduces an asymptotic maximum damage to the armor layer under a constant regular wave attack. It also depends on a mean damage parameter consisting of the number of regular waves causing 63% of the maximum asymptotic damage, which is linked to the concept of mean lifetime of the structure. Medina compared the results with the models from Teisson and Vidal (based on different assumptions but accomplishing most of the five aforementioned conditions) and applied the new method to a real case: the partial failure of Zierbana breakwater (Port of Bilbao, Spain) under construction in February 1996. In Gómez-Martín and Medina [71], the wave-to-wave exponential model was slightly modified and the mean damage parameter was found to be dependent on Iribarren’s number. They also designed a neural network (NN) applicable to random waves in non-stationary conditions, finding that the estimation of accumulated armor damage using both wave-to-wave exponential method and NN model showed a good agreement to damage observations.

2.5. The 90s: Probabilistic Approaches and Damage Progression Models. Melby and Kobayashi’s Model

After the intensive research on stability of rubble mound breakwaters during the 70s and 80s, there was a new attempt in the 90s to compile and check the most relevant models on armor stability:

(1)

In 1991, Pfeiffer [72] compared Hudson [25], Hedar [47], Losada and Giménez-Curto [44], and VdM [50] formulae with prototype data from Burns Harbor breakwater (Indiana).

(2)

In 1992, the final report from PIANC [73] provided a review on random wave’s armor stability models. In the same year, Koev [74] studied statistically a homogeneous set of 21 armor layer stability formulae developed under regular waves, and proposed a regression model valid for 0.04 ≤ H/L ≤ 0.1 and 1.1 ≤ cotα ≤ 20.

(3)

In 1994, Vidal et al. [75] revised the available methodologies for the calculation of the armoring hydraulic stability, both for breakwater heads and trunks. Formulations and design recommendations for berm breakwaters, low-crested breakwaters and conventional breakwaters were also included.

At the same time, some clues about physical modeling and laboratory techniques in coastal engineering were established based on the experience gathered by the scientific community over 60 years:

(1)

Hughes [76] presented in 1993 a publication of reference in physical modeling. It discussed the principles of dimensional analysis, scale effects, and similitude criteria including specific similitude requirements for different coastal structures. Also, considerations about movable-bed models, generation of gravity waves in laboratory, and a discussion on laboratory measurements and data analysis were included. Scale effects were further addressed in 2004 by Tirindelli et al. [77], who focused on wave impacts, run-up, overtopping, structure deformation, porous flow, and flow forces on plants and organisms.

(2)

Davies et al. [78] summarized in 1994 the different methods for damage measurement, describing the different techniques available for this purpose. Using experimental data, a comparison between damage measured with a profiler and damage defined by stone counting was carried out, with good agreement for low levels of damage. In addition, the sliding failure of the armor layer was investigated.

(3)

Burcharth et al. [79] suggested in 1999 a methodology for scaling core material in small scale rubble mound breakwater models. As it was first identified by Broderick and Ahrens [9], this kind of experiments can be subjected to significant scale effects when the flow type through the model core is different than in prototype: in a Froude small-scale model the flow through the core is usually laminar whereas in a full-scale core the flow is mainly turbulent. Indeed, Hegde and Srinivas [80] demonstrated experimentally that as core porosity is increased the stability is also increased. De Jong [81] ratified the obtention of lower values of damage after scaling the core according to Burcharth’s methodology. Additional information on core permeability and damage, together with prototype data, can be found in Reedijk et al. [82].

Similarly to the effects after the effort made in the 70s to compile the available formulae and to provide some design recommendations, the new attempt to get together armor stability models and laboratory techniques coincides with a prolific decade, in which two main conceptual and methodological progresses were made: the consolidation of the probabilistic approach and the origin of the first damage progression models.

At first, the introduction of the probabilistic approach faced a controversial acceptance in a profession where experience and expert criteria were fundamental regarding breakwater design. Some of the reasons of this initial skepticism in favor of deterministic methodologies based on safety factors were (1) the negative connotation of the term "probability of failure" (which can be referred euphemistically by the opposite concept "reliability"), (2) the lack of confidence due to uncertainty in the calculation of the strict failure bound (see Figure 1), and (3) the limited availability of probabilistic information from the different factors affecting the structural response, as well as from the structural response itself.

In spite of the difficulties, the probabilistic approach, widely extended nowadays, was materialized in the design codes. As an example, reliability methodologies were included in the European PROVERBS [83], the Spanish ROM 0.0-01 [84], or the North-American Coastal Engineering Manual [1]. Further investigations on the probabilistic character of armor stability in the 90s are given below:

(1)

Wang and Peene [85] possibly headed in 1990 the first attempt on the development of a fully probabilistic model of rubble mound armor stability, based on the stochastic nature of both wave forces and resistance forces. After pointing out that “the random nature of the resistant force offered by the armor blocks has not been seriously addressed at all,” they examined the behavior of interlocking resistance and the random nature of breakwater stability through laboratory experiments. A theoretical probabilistic model containing six random variables was proposed, which behaviors were not all known at that time, not even hitherto. Instead of computing wave loading in a conventional way, they calculated the resistance of the armor layer by pulling out the units with a motorized lift line (static stability test), recording the force history with a load cell. They applied a modified version of the Kolmogorov-Smirnov D-test on five data sets (with about 120 samples each) with tetrapods and dolosse, after testing different bed slopes, pull-out directions and locations, placement methods and unit sized. They concluded that the resistance of tetrapods and dolosse could be treated as a random variable with a log-normal distribution. Furthermore, comparing these results with the ones carried out with stones, they cast a doubt on the stabilizing effects of interlocking properties of artificial units.

(2)

Carver and Wright [86] pointed out in 1991 the random variations in the stability response of stone-armored rubble mound breakwaters after carrying out stability tests with depth limited irregular waves. They concluded that “repeat testing is a must,” indeed repeating each spectrum six or seven times. Also, they registered how the lower stabilities occurred at the lowest values of h/L in shallow water, i.e., regarding the longer wave periods.

(3)

Medina [70], as detailed in Section 2.4, recommended in 1996 non-stationary stochastic models for being more adequate for modeling real waves.

(4)

Burcharth [87] published in 1997 the Chapter Reliability-Based Design of Coastal Structures as part of a book which summed up the advances in coastal and ocean engineering. Different sources of uncertainty were identified, together with a probabilistic methodology for single failure mode probability analysis, including formulations and examples with Level III methods, Level II methods and FORM (first-order reliability method). Additionally, further examples of probabilistic methodologies applied to breakwater design can be found in Castillo et al. [88,89], Mínguez et al. [90], Tørum et al. [91], or Gouldby et al. [92]. Burcharth [87] also set a discussion on the probability analysis of failure mode systems, typically faced using fault trees. Some practical examples on this topic can be found in van Gent and Pozueta [93], who studied the rear-side stability of rubble mound structures after being overtopped, or Campos et al. [94], who addressed the effects of a fuse parapet failure in other failure modes of a caisson breakwater using Monte Carlo simulations.

(5)

Similar to Wang and Peene [85], Hald [8] and Hald and Burcharth [95] also chose an alternative approach away from correlating wave parameters directly to armor layer stability. For this purpose, they studied the wave-induced loading by means of a force transducer connected to an 8mm steel rod attached to a single stone of average size made in coated plastic foam. The rest of the armor layer was made of conventional natural stones. Largest forces were found to take place in a normal direction and upslope, and a dimensionless force model of the normal and the peak force was proposed as a function of wave parameters. A log-normal distribution was found to be suitable for describing the limit pullout forces. From the force model, they derived a stability model in 1998 based on a lifting criterion, obtaining comparable scatter with respect to the equations from Hudson [25] and VdM [50].

The evolution toward probabilistic approaches in the 90s revealed deeper comprehension on the wave-structure interaction and, consequently, it happened together with an improvement in the conception of hydraulic stability models. Traditional hydraulic stability formulae were aimed to characterize the initiation of movement, which was useful for estimating the design weight of armor units but failed to give information about the evolutionary behavior of rubble mounds. Some of these models, such as the formulae from Iribarren [12] or Hudson [25], were complemented with tabulated or graphic information relating damage percentage with coefficients such as the aforementioned H/H_D=0. These approaches are helpful for design purposes, but they are not time dependent and they assume starting from a non-damaged structure. In a context where, not only the design was meant to be addressed, but also the evolution of the structure regarding maintenance strategies and useful-life total costs, hydraulic "static" stability concept moved on toward damage progression models. The latter are aimed to predict the evolution of rubble mound’s geometry by means of a quantitative damage parameter, usually the dimensionless erosion area (S).

Probably the first damage progression model was due to VdM [50] in 1988. As it was detailed in Section 2.2, these widely known formulae include the possibility of being applicable to a probabilistic design, although the way of adapting the equation for this purpose presents some shortcomings from a probabilistic point of view. For example, the scatter is all focused on one parameter, instead of designing a fully probabilistic formulation with an analytic CDF (cumulative distribution function). As it was stated by some authors, VdM’s formulae are not valid for the long-term, as they tend to overestimate damage for more than 7000 waves. Also, despite providing reasonable estimates of the stability number for two specified damage levels (damage initiation and exposure of the filter layer), they failed to yield predicted damage levels with accuracy.

Further damage progression models proposed in the 90s were already mentioned and are summarized in Table 2. The formulation of Teisson [66] was developed using a step-by-step methodology assuming that a storm can be described by a sequence of steps with a certain duration each. Kaku et al. [56] proposed an exponential model with similarities between the stability number and the Shields parameter used in sediment transport. Medina [70] and Gómez-Martin and Medina [71] also suggested exponential models based on a wave-to-wave approach. However, probably one of the main contributions to damage progression formulations was due to Melby and Kobayashi, referred as M&K in this paper. Their first model was published in 1998 [96,97] and re-formulated in 1999 [98] for allowing non-zero initial damage values:

$${\left[\overline{S}\left(t\right)\right]}^{1/b}={\left[\overline{S}\left(t_n\right)\right]}^{1/b}+{\left(a_s{N}_s^5\right)}^{1/b}\frac{t-t_n^{}}{T_m^{}}; t_n\le t\le t_{n+1}$$

(9)

where $\overline{S}\left(t_n\right)$ is the known damage at the time t_n, N_S is the stability number based on the highest one-third wave heights from a zero-upcrossing analysis, T_m is the mean period, b is an empirical coefficient introduced for long-duration tests, and a_s is related to breakwater slope, permeability and an empirical coefficient derived from the tests. Therefore, it is an iterative damage progression model that allows the calculation of the damage at the instant t_n+1 based on the damage level at the instant t_n and the incident wave conditions between t_n and t_n+1 represented by a constant value of H_S and T_m.

The experiments of M&K, with a limited validity range, were conducted at the WES. A complete description of them, together with an analysis of incipient stability of armor units, damage measurement and damage definition, is detailed in the technical report of Melby [99]. Mean and standard deviation of different damage descriptors were used to characterize the trends, variability and ranges of the damaged profile descriptors. The damage variability ${\sigma }_S$ tends to increase with the mean damage value $\overline{S}$. However, the relative variability defined as $\sigma /\overline{S}$ decreases when increasing $\overline{S}$, i.e., the higher the damage is, the lower is the relative error of the estimation.

Contrary to other authors who maintained that damage evolution reaches an equilibrium, M&K observed that damage keeps on increasing with a lower rate. Indeed, Equation (9) does not have an asymptotical trend and, thus, a large number of small waves between two storms would theoretically provoke a damage increase. For that reason, they recommend setting up a critical stability number so that the damage would not increase beyond this lower limit. This concept is similar to the equilibrium damage level introduced by Kaku et al. [56]

The damage progression model from M&K is especially helpful in a lifecycle analysis because it allows engineers to balance initial costs with expected maintenance costs. Thus, it is also useful for rehabilitation or maintenance plans in already-built breakwaters. However, it presents a number of shortcomings, mainly pointed out by the authors themselves [100]. The results were reported to be conservative for most applications because they were based on severely breaking waves, a relatively steep beach slope and relatively impermeable core. However, damage initiation is unpredicted in all series studied and fails to yield accurate predictions for low damage levels. Furthermore, the model is aimed just at the mean damage evolution and, consequently, it does not provide information about the probability density function (PDF) of each estimated value of damage, which is needed for a more precise probabilistic analysis of damage progression.

2.6. Researches on Breakwater’s Damage in the XXI Century. Castillo et al.’s Model

The XXI century started with the publication of some aforementioned design codes with a high impact all over the world: the European PROVERBS in 1999 [83], the Spanish ROM 0.0-01 in 2001 [84] and the North-American CEM in 2002 [1]. Initiatives such as the European project on berm breakwater structures (MAS2-CT94-0087 [101]) also contributed to enhance the knowledge on the breakwaters response. There were also compiling efforts of the state-of-the-art such as the one focused on berm breakwaters from the Working Group 40 of PIANC [102]. In addition, the International Standard Organization published in 2007 [103] a new standard, ISO 21650 "Actions from Waves and Currents on Coastal Structures," which became the first one in coastal engineering.

During this recent period, some of the existing models were reassessed. The already referred studies from Van Gent et al. [61], Gómez-Martín and Medina [71], Vidal et al. [58,59], Mertens [60] or Verhagen and Mertens [62] are some examples of this trend. Besides, some concepts were re-analyzed. For example, Wolters and Van Gent [104] and Van Gent [105] addressed the effect of oblique wave attack on cube and rock armored rubble mound breakwaters in 2010 and 2014 respectively. They obtained results in agreement with the ones from Losada and Giménez-Curto [45] for gravity armor units, i.e., the stability of steep slopes under oblique wave attack is not worse than for normal incidence.

Additionally, new approaches were proposed. In 2003, Tomasicchio et al. [106] developed a conceptual model for evaluating armor stone abrasion and compared their results with the ones from Latham [107]. Eslami and Van Gent [108] presented in 2010 one of the few studies relating wave overtopping and rubble mound stability with combined loading of waves and current. Esteban et al. [109] assessed in 2012 the stability of rubble mound breakwaters against tsunami attacks, with experimental and real data. Medina et al. [110] studied in 2010 the influence of armor units’ placement on armor porosity and hydraulic stability, analyzing different placement grids of cubes and Cubipod^® by means of pressure clamps. Similar experiments were discussed by Pardo et al. [111] in 2013 after accounting for the three armor randomness indexes (ARIs) introduced by Medina et al. [112]. These experiments were measured with a laser scanner in order to quantify the randomness in the placement of cube and Cubipod^®. Further information on the influence of initial placement can be found in Yagci and Kapdasli [113,114], Gürer et al. [115], Van Buchem [116], or Medina et al. [117]. Recently, Marzeddu et al. [118] analyzed different approaches for taking into account wave storm representation on damage measurements. Also, Clavero et al. [119] proposed a methodology to improve breakwater design and to assimilate the data from different wave flumes after analyzing the bulk wave dissipation in the armor layer of rock and cube armored small-scale models.

It is also remarkable the advances in numerical models reproducing the physics of the main processes involved (see Farhadzadeh et al. [120] or García et al. [121] for some examples of the latest contributions). Although addressing the structural stability at the scale of individual units is still computationally overdemanding, studies such us Matsumoto et al. [122] have already been able to model relatively complex scenarios by computational fluid dynamics. Soft computing techniques, such as artificial neural networks (ANN) [123,124,125,126], genetic programming [127], or extreme learning machine (ELM) models [128], were also applied for predicting damage levels. They also offer the possibility of optimizing temporal and economical costs of small-scale physical modeling [129]. Probably the most ambitious project of this kind is the one from the International Breakwater Directory, aimed to expand the CLASH (acronym for “Crest Level Assessment of coastal Structures by full scale monitoring, neural network prediction and Hazard analysis on permissible wave overtopping”) overtopping database for including information about armor damage in rubble mound breakwaters and seawalls [7].

However, the lack of a stochastic damage progression model fully designed from a probabilistic approach, together with some damage accumulation modeling difficulties reported by M&K when re-formulating their model [98], motivated Castillo et al. to investigate on the basis for building consistent stochastic models. Instead of developing an empirical or semi-empirical formulation based on a set of experimental data, Castillo et al. proposed in 2012 [130] a dimensionless stochastic damage progression model with general validity, avoiding the selection of easy-to-use mathematical functions, which were replaced by those resulting from a set of properties to be satisfied. The model was proven to have a normal distribution. It was built after applying dimensional analysis (using Π Buckingham’s theorem), compatibility conditions and the central limit theorem. It is expressed in terms of the following cumulative distribution function (CDF):

$$F_{D^{\ast }\left(t^{\ast }\right)}\left(x\right)=\mathsf{\Phi }\left(\frac{{\left(x-\gamma \right)}^{1/b}-{\mu }_0-k{t}^{\ast }}{\sqrt{{\sigma }_0^2+r{t}^{\ast }}}\right)$$

(10)

where γ and b are breakwater dependent, k and r include wave action and μ₀ and σ₀ depend on the initial conditions. Note that damage is not normal according to the model, but transformed damage, i.e., (D*-γ)1/b, is instead.

An initial calibration of the model was carried out in Campos et al. [131,132]. The experimental results, with limited validity, pointed out the lack of repeatability of damage progression, as it was stated by other authors such us Benedicto et al. [133] or Van Gent et al. [134]. This not only confirms the random nature of damage, but also points out the need of reproducibility of the relatively common damage accumulation experiments. Reproducibility is necessarily linked to standardization. In this sense, a lack of standards in damage initiation and progression experiments is perceived by the authors in line with the conclusions from other authors such us Losada and Giménez-Curto [44] or Ota et al. [135].

As a conclusion from this section, it is expectable that main challenges in the near future for rubble mound designers and managers will be related to develop or calibrate decision-making tools capable of accurately characterizing the structural response, together with its uncertainties, during the whole useful life. Not only the knowledge gathered over the years would play a role, but also the huge possibilities opened by the increase in computing capacity and the affordability and accuracy of monitoring tools such as global navigation satellite systems [136] or photogrammetric restitutions using drones. Risk-based design and management, including the possible effects of climate change and the upgrading of already built structures [137,138,139], would definitely require a deeper comprehension about the spatial and temporal evolution of the stochastic nature of damage.

3. Conclusions

In this paper, a historical overview of damage initiation and progression in rubble mound breakwaters is presented, highlighting the most relevant approaches and main hydraulic stability formula and damage progression models up to date. The amount of literature on this topic is quite extensive, which has to do with the importance of hydraulic stability in the design and maintenance of this kind of coastal structures as well as with the great challenge that implies the characterization of damage initiation and damage progression rates. This challenge is related to the following summing-up points:

(1)

Hydraulic instability of armor layers is a complex process because of the stochastic nature of both wave loading, initiation of movement, and damage progression. The highly non-linear flow over the slope, involving wave breaking, together with the variable shape of armor units and their random placement deals necessarily with a probabilistic concept of damage initiation and progression, as it was pointed out by many authors.

(2)

The concept of “damage” in a rubble mound breakwater, understood as the partial or total loss of its functionality, is subjected to different interpretations. Not every structure is designed under the same functional requirements and there is a wide variety of typologies. In addition, not every structure presents the same fragility. For instance, what is considered as damage in single-layered structures might not be equivalent in multi-layered structures. Also, each type of armor unit presents a singular behavior against wave action. These facts complicate the correlation between quantitative damage descriptors and qualitative damage levels such as the ones defined in Vidal et al. [140].

(3)

The parameterization and measurement of damage has been just slightly addressed in this paper due to the extensive information available. Indeed, Part II of this review [6] is aimed entirely on this topic. The dimensionless erosion area (S) from Broderick and Ahrens [9] is the most widely damage index nowadays. However, it seems that there is not a worldwide standard on how to measure it, which is crucial for reproducibility and for a consistent comparison between the results from different laboratories. Not only that, regarding the random nature of damage and the frequent accomplishment of damage initiation/progression tests in coastal laboratories, this kind of experiments needs to be reproducible. Thus, a concise methodology ideally agreed by the scientific community and shared worldwide may be helpful.

(4)

As exposed in Section 2.4, an adequate selection of wave action parameters is crucial for correlating damage just with the waves of the irregular train directly responsible of the hydraulic instability of the armor layer.

(5)

Damage has a spatial component that cannot be completely addressed with the solely characterization of the well-known dimensionless erosion area (S). Most studies aimed to develop empirical or semi-empirical damage formulations are conducted on wave flumes, using physical profilers or visual counting for characterizing damage. However, the recent advances and affordability of scanning and photogrammetric techniques [141] allow nowadays a more complete analysis of the geometrical evolution of the armoring both in wave tank’s models and in full-scale prototypes. Also, the increase in the resolution of global navigation satellite systems and other technological advances such as drones, can help to foster damage monitoring in real structures.

(6)

Indeed, taking into account that the structural response is a random variable, its characterization is necessarily linked, not only to repeating experiments in laboratory, but also to monitoring prototypes. Only in this way the models can be properly calibrated.

As it was mentioned before, the future of rubble mound design is linked to risk-based approaches and advanced management strategies, which would require a deep comprehension about the spatial and temporal evolution of damage during the whole useful life. Damage progression probabilistic models, such as the one from Castillo et al. [130], offer the framework to make possible this evolution and, at the same time, highlight the need of enhancing the knowledge of the singular behavior of each structure. To accomplish this, full-scale monitoring and standardization will presumably play a key role in the near future.

Table 1. Review of hydraulic stability formulae (extended from Hald [8]). See Appendix B for the list of symbols.

RESEARCHER	FORMULA	RESEARCHER	FORMULA
de Castro (1933) [4]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left\{\frac{{\rho }_w}{K}{\left(\cot \alpha +1\right)}^2\sqrt{\cot \alpha -\frac{2}{S_r}}\right\}}^{1/3}$	Rybtchevsky (1964) [23]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K\frac{H}{L}\right)}^{1/3}\cos \alpha \sin \alpha$
Iribarren (1938) [5]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=\left(\cos \alpha -\sin \alpha \right){\left(\frac{{\rho }_w}{K}\right)}^{1/3}$	Iribarren(1965) [11]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=\left(\tan \varphi \cos \alpha -\sin \alpha \right){\left(\frac{1}{K}\right)}^{1/3}$
Mathews(1948)	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left\{K\frac{H}{T}{\left(\cos \alpha -0.75\sin \alpha \right)}^2\right\}}^{1/3}$	Metelicyna (1967) [24]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\cos \left(23º+\alpha \right)$
Epstein, Tyrrel (1949) [16]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\left(\tan \varphi -\tan \alpha \right)$	SPM (1973) [27]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}={\left(K_D\cot \alpha \right)}^{1/3}$
Hickson, Rodolf (1951) [17]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K\frac{H}{T}\right)}^{1/3}\tan \left(45°-\frac{\alpha }{2}\right)$	Losada, Gim.-Curto(1979) [44]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\gamma }_w^{1/3}{\left\{K_1\left(\xi -{\xi }_0\right)\exp \left(K_2\left(\xi -{\xi }_0\right)\right)\right\}}^{-1/3}$
Larras (1952) [13]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\left(\cos \alpha -\sin \alpha \right)\frac{\sinh \left(4\pi h/L\right)}{\left(2\pi H/L\right)}$	SPM (1984) [28]	$\frac{H_{1/10}}{\mathsf{\Delta }{D}_{n50}}={\left(K_D\cot \alpha \right)}^{1/3}$
Hudson, Jackson (1953) [18]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\frac{\left(\tan \varphi \cos \alpha -\sin \alpha \right)}{\tan \varphi }$	Hedar (1986) [47]	$\frac{H_b}{\mathsf{\Delta }{D}_{n50}}={\left(\frac{6}{\pi }\right)}^{1/3}\frac{K_2{f}_1\left(\gamma \right)\cos \alpha }{K_1\left(\frac{h_b}{H_b}+0.7\right)\left(\tan \varphi +2\right)}$ $f_1\left(\gamma \right)=permeability function$
Beaudevin (1955) [14]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left\{K\left(\frac{\cot \alpha -0.8}{1.12-0.15\cot \alpha }\right)\right\}}^{1/3}$	Medina and McDougal (1990) [55]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=\frac{1.86}{1.27}\sqrt{\frac{2}{\ln N_w}}{\left(K_D\cot \alpha \right)}^{1/3}$
Hudson(1958) [22]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K_D\cot \alpha \right)}^{1/3}$	Wang and Peene (1990) [85]	$\frac{\tilde{w}}{W_D}=\frac{{\left(f_c\right)}_n^3{\left[C_S/C_V\right]}_n^3{H}_n^3}{R_n^3{D}_{n50}^3{\left[f\left(\theta \right)\right]}_n}$ $\begin{array}{l}see Wang and Peene(1990)\\ for the list of symbols\end{array}$
Goldschtein, Kononenko (1959) [21]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K_{}{\tan }^{1.83}\alpha \right)}^{1/3}$	Koev (1992) [74]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(\left(\frac{\cot {\left(\alpha \right)}^{K_2}}{K_1}\right){\left(\frac{H}{L}\right)}^{K_3}\right)}^{1/3}$
SN-92-60(1960) [22]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K\frac{H}{L}\sqrt{1+{\cot }^3\alpha }\right)}^{1/3}$	Hald and Burcharth (2000) [95]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}={\left(\frac{1}{K_1}\frac{1+{\sin }^2\alpha }{{\tan }^{K_2}\alpha }\right)}^{\frac{1}{K_3}}$
Svee(1962) [20]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\cos \alpha$	-

Table 1. Review of hydraulic stability formulae (extended from Hald [8]). See Appendix B for the list of symbols.

RESEARCHER	FORMULA	RESEARCHER	FORMULA
de Castro (1933) [4]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left\{\frac{{\rho }_w}{K}{\left(\cot \alpha +1\right)}^2\sqrt{\cot \alpha -\frac{2}{S_r}}\right\}}^{1/3}$	Rybtchevsky (1964) [23]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K\frac{H}{L}\right)}^{1/3}\cos \alpha \sin \alpha$
Iribarren (1938) [5]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=\left(\cos \alpha -\sin \alpha \right){\left(\frac{{\rho }_w}{K}\right)}^{1/3}$	Iribarren(1965) [11]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=\left(\tan \varphi \cos \alpha -\sin \alpha \right){\left(\frac{1}{K}\right)}^{1/3}$
Mathews(1948)	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left\{K\frac{H}{T}{\left(\cos \alpha -0.75\sin \alpha \right)}^2\right\}}^{1/3}$	Metelicyna (1967) [24]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\cos \left(23º+\alpha \right)$
Epstein, Tyrrel (1949) [16]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\left(\tan \varphi -\tan \alpha \right)$	SPM (1973) [27]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}={\left(K_D\cot \alpha \right)}^{1/3}$
Hickson, Rodolf (1951) [17]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K\frac{H}{T}\right)}^{1/3}\tan \left(45°-\frac{\alpha }{2}\right)$	Losada, Gim.-Curto(1979) [44]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\gamma }_w^{1/3}{\left\{K_1\left(\xi -{\xi }_0\right)\exp \left(K_2\left(\xi -{\xi }_0\right)\right)\right\}}^{-1/3}$
Larras (1952) [13]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\left(\cos \alpha -\sin \alpha \right)\frac{\sinh \left(4\pi h/L\right)}{\left(2\pi H/L\right)}$	SPM (1984) [28]	$\frac{H_{1/10}}{\mathsf{\Delta }{D}_{n50}}={\left(K_D\cot \alpha \right)}^{1/3}$
Hudson, Jackson (1953) [18]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\frac{\left(\tan \varphi \cos \alpha -\sin \alpha \right)}{\tan \varphi }$	Hedar (1986) [47]	$\frac{H_b}{\mathsf{\Delta }{D}_{n50}}={\left(\frac{6}{\pi }\right)}^{1/3}\frac{K_2{f}_1\left(\gamma \right)\cos \alpha }{K_1\left(\frac{h_b}{H_b}+0.7\right)\left(\tan \varphi +2\right)}$ $f_1\left(\gamma \right)=permeability function$
Beaudevin (1955) [14]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left\{K\left(\frac{\cot \alpha -0.8}{1.12-0.15\cot \alpha }\right)\right\}}^{1/3}$	Medina and McDougal (1990) [55]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=\frac{1.86}{1.27}\sqrt{\frac{2}{\ln N_w}}{\left(K_D\cot \alpha \right)}^{1/3}$
Hudson(1958) [22]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K_D\cot \alpha \right)}^{1/3}$	Wang and Peene (1990) [85]	$\frac{\tilde{w}}{W_D}=\frac{{\left(f_c\right)}_n^3{\left[C_S/C_V\right]}_n^3{H}_n^3}{R_n^3{D}_{n50}^3{\left[f\left(\theta \right)\right]}_n}$ $\begin{array}{l}see Wang and Peene(1990)\\ for the list of symbols\end{array}$
Goldschtein, Kononenko (1959) [21]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K_{}{\tan }^{1.83}\alpha \right)}^{1/3}$	Koev (1992) [74]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(\left(\frac{\cot {\left(\alpha \right)}^{K_2}}{K_1}\right){\left(\frac{H}{L}\right)}^{K_3}\right)}^{1/3}$
SN-92-60(1960) [22]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}={\left(K\frac{H}{L}\sqrt{1+{\cot }^3\alpha }\right)}^{1/3}$	Hald and Burcharth (2000) [95]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}={\left(\frac{1}{K_1}\frac{1+{\sin }^2\alpha }{{\tan }^{K_2}\alpha }\right)}^{\frac{1}{K_3}}$
Svee(1962) [20]	$\frac{H}{\mathsf{\Delta }{D}_{n50}}=K_{}^{1/3}\cos \alpha$	-

Table 2. Review of damage progression’s models. See Appendix B for the list of symbols.

RESEARCHER	FORMULA
Van der Meer (1988) [51]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=\{\begin{array}{l}6.2{\xi }_0^{-0.5}{P}^{0.18}{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}\hspace{1em}\to for plunging waves\\ {\xi }_0^P\sqrt{\cot \alpha } P^{-0.13}{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2} \to for surging waves\end{array}$
Teisson (1990) [66]	$N_d\left(t\right)=A{\left({\displaystyle \sum H_{S,i}^{B/C}\mathsf{\Delta }{t}_i}\right)}^C$
Kaku et al. (1991) [56]	$S=S_e\left(1-e^{-K\cdot N_w}\right)$
Medina (1996) [70]	$D\left(H,T,N_w\right)=D_0\left(H,T\right)\left[1-e^{\frac{-N_w}{n63\%}}\right]$
Melby and Kobayashi (1998) [96]	$\overline{S}\left(t\right)=\overline{S}\left(t_n\right)+a_s{N}_s^5\left(\frac{t^b-t_n^b}{T_m^b}\right); t_n\le t\le t_{n+1}$
Melby and Kobayashi (1999) [98]	${\left[\overline{S}\left(t\right)\right]}^{1/b}={\left[\overline{S}\left(t_n\right)\right]}^{1/b}+{\left(a_s{N}_s^5\right)}^{1/b}\frac{t-t_n^{}}{T_m^{}}; t_n\le t\le t_{n+1}$
Melby and Hughes (2003) [3]	$N_m=\{\begin{array}{l}0.5{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}{P}^{0.18}{\left(\cot \alpha \right)}^{0.5}\hspace{1em}\to for plunging waves\\ 0.5{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}{P}^{0.18}{\left(\cot \alpha \right)}^{0.5-P}{s}_m{}^{-P/3} \to for surging waves\end{array}$
Van Gent et al. (2003) [61]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=1.75{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}\sqrt{\cot \alpha }\left(1+\frac{D_{n50,core}}{D_{n50}}\right)$
Gómez-Martín and Medina (2004) [71]	$D\left(H,Ir,N_w\right)=D_0\left(H,Ir\right)\left[1-2^{\frac{-N_w}{n50\%}}\right]$
Vidal et al. (2006) [59]	$\begin{array}{l}\frac{H_{50}}{\mathsf{\Delta }{D}_{n50}}=\{\begin{array}{l}4.44P^{0.18}{S}^{0.2}{\xi }_0^{-0.5}\hspace{1em} \to for {\xi }_0<{\xi }_{mc} and \cot \alpha \le 4\\ 0.716\sqrt{\cot \alpha }{P}^{-0.13}{S}^{0.2}{\xi }_0^P \to for {\xi }_0\ge {\xi }_{mc} or \cot \alpha \ge 4\end{array}\\ {\xi }_{mc}={\left(6.2P^{0.31}\sqrt{\tan \alpha }\right)}^{\frac{1}{P+0.5}}\end{array}$
Castillo et al. (2012) [130]	$F_{D^{\ast }\left(t^{\ast }\right)}\left(D\right)=\mathsf{\Phi }\left(\frac{{\left(D-\gamma \right)}^{1/b}-{\mu }_0-k{t}^{\ast }}{\sqrt{{\sigma }_0^2+r{t}^{\ast }}}\right)$

Table 2. Review of damage progression’s models. See Appendix B for the list of symbols.

RESEARCHER	FORMULA
Van der Meer (1988) [51]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=\{\begin{array}{l}6.2{\xi }_0^{-0.5}{P}^{0.18}{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}\hspace{1em}\to for plunging waves\\ {\xi }_0^P\sqrt{\cot \alpha } P^{-0.13}{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2} \to for surging waves\end{array}$
Teisson (1990) [66]	$N_d\left(t\right)=A{\left({\displaystyle \sum H_{S,i}^{B/C}\mathsf{\Delta }{t}_i}\right)}^C$
Kaku et al. (1991) [56]	$S=S_e\left(1-e^{-K\cdot N_w}\right)$
Medina (1996) [70]	$D\left(H,T,N_w\right)=D_0\left(H,T\right)\left[1-e^{\frac{-N_w}{n63\%}}\right]$
Melby and Kobayashi (1998) [96]	$\overline{S}\left(t\right)=\overline{S}\left(t_n\right)+a_s{N}_s^5\left(\frac{t^b-t_n^b}{T_m^b}\right); t_n\le t\le t_{n+1}$
Melby and Kobayashi (1999) [98]	${\left[\overline{S}\left(t\right)\right]}^{1/b}={\left[\overline{S}\left(t_n\right)\right]}^{1/b}+{\left(a_s{N}_s^5\right)}^{1/b}\frac{t-t_n^{}}{T_m^{}}; t_n\le t\le t_{n+1}$
Melby and Hughes (2003) [3]	$N_m=\{\begin{array}{l}0.5{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}{P}^{0.18}{\left(\cot \alpha \right)}^{0.5}\hspace{1em}\to for plunging waves\\ 0.5{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}{P}^{0.18}{\left(\cot \alpha \right)}^{0.5-P}{s}_m{}^{-P/3} \to for surging waves\end{array}$
Van Gent et al. (2003) [61]	$\frac{H_S}{\mathsf{\Delta }{D}_{n50}}=1.75{\left(\frac{S}{\sqrt{N_w}}\right)}^{0.2}\sqrt{\cot \alpha }\left(1+\frac{D_{n50,core}}{D_{n50}}\right)$
Gómez-Martín and Medina (2004) [71]	$D\left(H,Ir,N_w\right)=D_0\left(H,Ir\right)\left[1-2^{\frac{-N_w}{n50\%}}\right]$
Vidal et al. (2006) [59]	$\begin{array}{l}\frac{H_{50}}{\mathsf{\Delta }{D}_{n50}}=\{\begin{array}{l}4.44P^{0.18}{S}^{0.2}{\xi }_0^{-0.5}\hspace{1em} \to for {\xi }_0<{\xi }_{mc} and \cot \alpha \le 4\\ 0.716\sqrt{\cot \alpha }{P}^{-0.13}{S}^{0.2}{\xi }_0^P \to for {\xi }_0\ge {\xi }_{mc} or \cot \alpha \ge 4\end{array}\\ {\xi }_{mc}={\left(6.2P^{0.31}\sqrt{\tan \alpha }\right)}^{\frac{1}{P+0.5}}\end{array}$
Castillo et al. (2012) [130]	$F_{D^{\ast }\left(t^{\ast }\right)}\left(D\right)=\mathsf{\Phi }\left(\frac{{\left(D-\gamma \right)}^{1/b}-{\mu }_0-k{t}^{\ast }}{\sqrt{{\sigma }_0^2+r{t}^{\ast }}}\right)$

Table 3. Timeline of some of the most relevant models on armor layer’s stability, together with the parameters/properties accounted for and the innovations introduced by each of them. See Appendix B for the list of symbols.

RESEARCHER	Wave Action					Breakwater				Others
	R	I	T	Ir	Nw	h	a	AT	P	D0	DPM	CDF
Iribarren (1938) [5]
Mathews (1948)
Larras (1952) [13]
Hudson (1958) [22]
SPM (1973) [27]		H1/3
Losada & G.C. (1979) [44]
SPM (1984) [28]		H1/10
Hedar (1986) [47]		Hb
VdM (1988) [51]		H1/3
Teisson (1990) [66]		H1/3
Kaku et al. (1991) [56]		H1/3
Medina (1996) [70]		Hi
M & K (1998, 1999) [96,98]		H1/3
G.M.& Medina (2004) [71]		Hi
Van Gent et al. (2003) [61]		H1/3
Vidal et al. (2006) [59]		H50
Castillo et al. (2012) [130]		H

Table 3. Timeline of some of the most relevant models on armor layer’s stability, together with the parameters/properties accounted for and the innovations introduced by each of them. See Appendix B for the list of symbols.

RESEARCHER	Wave Action					Breakwater				Others
	R	I	T	Ir	Nw	h	a	AT	P	D0	DPM	CDF
Iribarren (1938) [5]
Mathews (1948)
Larras (1952) [13]
Hudson (1958) [22]
SPM (1973) [27]		H1/3
Losada & G.C. (1979) [44]
SPM (1984) [28]		H1/10
Hedar (1986) [47]		Hb
VdM (1988) [51]		H1/3
Teisson (1990) [66]		H1/3
Kaku et al. (1991) [56]		H1/3
Medina (1996) [70]		Hi
M & K (1998, 1999) [96,98]		H1/3
G.M.& Medina (2004) [71]		Hi
Van Gent et al. (2003) [61]		H1/3
Vidal et al. (2006) [59]		H50
Castillo et al. (2012) [130]		H