Eight separate experimental models were investigated in this study representing three different geometrical setups: placed riprap layer on a ramp (RP), half a dam with placed (HP) or dumped (HD) riprap on a rockfill shoulder, and full dam rockfill dam models protected with placed riprap (FP). One of the half dam models (HPT) had a support at the downstream toe of the riprap layer. These models have been well introduced in Dezert et al. (2022a). All these physical model tests were carried out in the hydraulic laboratory at NTNU. The models were built in a flume that was 25 m in length, 1 m in width, and 2 m in height. Within the flume, an elevated platform of 0.35 m was created using aluminum boxes so as to minimize the effects of backwater. A geotextile was layered on this platform so as to replicate friction at the toe. Similarly, geotextile was layered on the ramp for the model RP. Froude similarity was used to design this 1:10 scale conceptual model with a dam height of 1 m. In accordance with the existing protective layers on Norwegian dams, the riprap stones were placed at an angle of β ≈ 60 ° with respect to the platform on the slopes and at an angle of β ≈ 90 ° at the crest section of the dam. Both the shell (supporting rockfill) material and the material used within the protective layer (filter and riprap) were dimensioned in consonance with Norwegian guidelines for the construction of embankment dams by the NVE (NVE, 2012). The gradation curves for the material used in the model structures were downscaled by a ratio of 1:10. Nonetheless, the material was aligned with the coarser interval of the gradation curves as the finest material (<0.5 mm) had to be excluded due to technical limitations within the laboratory setup. In order to have a sectional view of the dam model, four glass panes (1 m width and 2 m height) were merged into the flume structure. This transparent section of the flume was used to record the progression of the dam breach. The overall sectional and planar view of the model setup has been illustrated in Figure 2.

The placed riprap layer on a ramp model (RP) consists of having a simple filter layer and riprap layer placed on a ramp representing the downstream half of the dam. A geotextile was placed on top of the ramp and a 0.1 m thick layer of filter was placed upon it. The filter layer consisted of stones with a density of ρ_s,Filter = 3,050 kgm^-3. The thickness of the filter layer was in accordance with the NVE guidelines (NVE, 2012). The riprap consisted of rhyolites, an igneous rock with a density of ρ_s,Riprap = 2,710 kgm^-3, and an average mass of 0.24 kg (Hiller et al., 2018). The median riprap stone diameter was d₅₀ = 0.057 m. The value was calculated as d₅₀ = (abc)^1/3, where a,b, and c represent the longest, intermediate, and shortest axis of the stone respectively (Dezert et al., 2022a; Ravindra et al., 2021).

In order to gain a comprehensive assessment of rockfill dam stability during throughflow and overtopping events the half dam models were developed adding a rockfill shoulder and placing the filter and riprap layers on the downstream rockfill slope (Ravindra and Sigtryggsdóttir, 2021). This would provide a better insight into the interactions between different components of the dam and pinpoint critical elements and locations where dam failure could originate. Initially, these half dams were built without any riprap protection but instead incorporated varying toe configurations which have been investigated by Kiplesund et al. (2021). Further studies amalgamating filter and riprap layers are presented in the NVE report (Ravindra and Sigtryggsdóttir, 2021). These studies included both unsupported placed and dumped riprap which has been further analyzed in this study. To simplify the model design the central core and filter layer were represented by an impervious aluminum element. The core element has a height of 0.8 m, thus having 0.2 m of shell material above it (dam height 1 m). This was in agreement with the requirements of Norwegian regulations for dams with central moraine cores. The shell material consists of graded rockfill material with a density of ρ_s,Shell = 2,720 kgm⁻³ while the filter layer consists of stones with a density of ρ_s,Filter = 2,900 kgm⁻³. The same riprap stones as the previous study were used for this setup.

The shoulder of the dam was constructed in layers of 0.1 m thickness. After the material had been laid, a tamper was used to uniformly compact the layer. The tamper weighing 4.54 kg was dropped ten times in an overlapping pattern from a height of approximately 0.1 m onto the layer surface (Ravindra and Sigtryggsdóttir, 2021). The half dam with toe support (HPT) model setup is similar to the setup of the half dams without a toe. The only difference is the inclusion of a fixed-toe support. Load cells were also incorporated with the fixed-toe support in order to measure the loads at the toe. The same shell material was used while the density of the riprap stones used was ρ_s,Riprap = 2,600 kgm⁻³, and the filter material density was ρ_s,Filter = 3,050 kgm⁻³ (Dezert et al., 2022a; Dezert and Sigtryggsdóttir, 2024).

As a continuation of the previous tests, full dam profiles were investigated next. For these models, the impermeable element was represented by a 1 mm thick rubber membrane. The black styrene butadiene rubber sheet was taped and sealed on the bottom and the wall side using Wurth duct tape while transparent Stokvis joining tape was used on the glass pane side. Furthermore, the upstream shoulder was also constructed. The same shell material was used with a density of ρ_s,Shell = 2,720 kgm⁻³. The shell material alone had a mass of 4,500 kg (Kiplesund et al., 2023).

For each setup, the inflow during the breach was not set. The models were subjected to sequential overtopping with incremental increases in discharge which was held constant for a fixed time interval until complete failure was reached for a critical inflow discharge. This information has been summarized in Table 1. The models were placed sufficiently upstream from the pump outlet to ensure calm flow conditions in the upstream section of the dam models.

Table 1 provides a summary of the test models investigated in this paper. Two full dam models (FP) had a pilot channel with bottom and top width W_c,bb and W_c,tt respectively and height H_c (Dezert et al., 2024). A pilot channel was constructed on the crest of the dam on the glass pane side to induce the breach on that side where it can be both observed and documented. However, the pilot channel did not influence the breach initiation in the protected dams.

A visual representation of the models summarized in Table 1 has been presented in Figure 3. The rubber sheet provided sufficient sealing and met predefined criteria such as acceptable impact of the watertight membrane, realistic phreatic surface, ease of construction, uniformity across tests, etc. (Kiplesund et al., 2023). Furthermore, the shell material which has been used for these models has been tested with varying toe configurations to study the throughflow in rockfill dams, there have been both physical and numerical studies (Kiplesund et al., 2021; Smith et al., 2021). These studies provided satisfactory results with regard to pore pressure distribution and realistic phreatic line. Even though the gradation was conducted with the minimum permissible particle size set at 0.5 mm, the gradation curve was determined based on data derived from real rockfill dams and checked against criteria in the NVE guidelines (NVE, 2012; Kiplesund et al., 2021). Moreover, Darcian or linear flow is rarely encountered in real-world applications with flow through rockfill (Leps, 1973). Linear Darcy flow theory is only applicable to flow through small grains around 0.5 mm in size (Wilkins, 1955). Due to characteristics such as porosity, particle shape, particle size, roughness, and tortuosity of rockfill material, flow within rockfill structures usually deviates from the linear flow regime, resulting in high velocities (Siddiqua et al., 2011).

Figure 4 provides a summary of the PIV workflow, detailing various steps involved in the process, from data acquisition to post-processing and analysis. The boxes in the blue dotted region in Figure 4 are steps that were carried out within PIVlab (Thielicke and Stamhuis, 2014). The workflow has been assigned numerical identifiers to each step. These step numbers will be referenced in the subsequent methodology section, where each step will be described in detail.

The cameras have been set up so that they record the entire breaching process (Step 1). Typically, this entails documenting the process from the initial filling of the flume reservoir upstream of the model dam all the way through to the complete breach of the model embankment dam. As this entire process takes quite an extended amount of time the files are considerable in size (>20 GB). Hence, for this study, the files were trimmed (Step 2) and only the time intervals of interest were selected (the main breach). Within this interval, all frames available were taken into account and analyzed in PIVlab. PIVlab is a MATLAB toolbox presented in a GUI form, designed for calculating velocities within image data.

After the time interval of interest was selected, a mask was applied in PIVlab to exclude some of the areas (Step 3). These excluded areas mostly include the elevated platform, the columns, the aluminum core when present, and most of the backdrop area behind the dam model. As these objects did not move from one frame to another, it was not necessary to include them in the PIV analysis. Figure 5 presents the mask being applied to model FP_01.

The images were then pre-processed using the Contrast Limited Adaptive Histogram Equalization (CLAHE) technique (Step 4). This technique improves the contrast of the images while overcoming the limitations of Adaptive Histogram Equalization (AHE) where it adds noise as a result of over-amplification at analogous regions of the image. The CLAHE technique does not take the histogram of the entire image and instead takes the histograms of smaller tiles which is crucial as there is no assurance that the complete image will be uniformly exposed (Westerweel, 1993). Using this technique substantially increases the likelihood of identifying valid vectors between the images (Shavit et al., 2007).

Next, a PIV algorithm was chosen under the PIV settings. The most important part of the PIV analysis is the cross-correlation. The cross-correlation is used to generate the most probable spatial displacement of the particles within the interrogation areas of successive singly exposed images (Adrian and Westerweel, 2010; Thielicke, 2014). The discrete cross-correlation function used for this purpose is Equation 1 (Huang et al., 1997): C m , n = ∑ i ∑ j P I P 1 i , j P I P 2 i − m , j − n (1)

Where PIP ₁ and PIP ₂ (Particle Image Pattern) correspond to the image masks in the first and second images of the single exposed image pair. This correlation matrix shows where the particles are most likely to have been displaced from PIP ₁ to PIP ₂ based on the location of the intensity peak in matrix C (m,n) (Huang et al., 1997).

This equation can be solved using two general methods, either in the frequency domain or in the spatial domain (Gonzalez and Wintz, 1987; Willert and Gharib, 1991). Evaluation in the frequency domain is done by a Fast Fourier Transform (FFT) algorithm while the approach on the spatial domain is done using the Direct Cross Correlation (DCC) technique (Adrian, 1991). PIVlab uses Discrete Fourier Transform (DFT) which is calculated using the FFT algorithm (Step 5).

DCC supports the use of interrogation areas of different sizes (Stamhuis, 2006). According to (Huang et al., 1997) DCC is able to generate more precise calculations as compared to the standard DFT method where both the random and systematic errors are reduced significantly. However, the main drawback with regard to the DCC technique is that it has a very high computational cost. This aspect is most significant when the interrogation area considered is large (Soria, 1996). This disadvantage can be solved using the DFT technique which uses interrogation areas of equal size (Raffel et al., 2007). This technique can however result in the increase of background noise which is a result of the loss of information because of identical interrogation areas. This inevitably decreases the accuracy of the results as it becomes difficult to identify the intensity peaks. To overcome this drawback and also maintain a minimum computational cost, multiple passes can be run on the same data set (Westerweel et al., 1997). To further improve this, an algorithm was proposed (Scarano and Riethmuller, 1999) that enhances the dynamic range and spatial resolution of the vector map. This is done by implementing successive interrogation passes where the interrogation grid is refined with every step (Step 6). The initial step uses a larger interrogation area and is thus able to read higher displacements, in the subsequent steps this area is reduced (in this setup the area was reduced by 50%, and three passes were implemented). The interrogation area started from 64px, the second pass was set at 32px, and the final pass at 16px. A smaller stop window would have given a better resolution however, caution was exercised not to use too small of an interrogation area as this would notably increase noise, hence generating inaccurate correlations. Furthermore, a smaller stop window would significantly increase the computational time. Subsequently, an individual analysis was performed for each frame (Step 7).

It is best to follow the above step with calibration before any data validation is done (Step 8). Up until this step, the units that have been used in PIVlab are pixels per frame. When calibrating it is possible to select an image separately or simply use a frame from the video recorded. Calibration can be applied using the known distance between two points and the time step between each frame. The coordinate system is established by defining the positive and negative directions for both x and y (Step 9). In this configuration, the flow direction is from left to right, hence, the positive x direction was chosen as right (streamwise direction), and the positive y direction was downwards (gravity direction).

Next, the data needs to be validated in order to filter any possible erroneous vectors. This was initially done through velocity-based validation (Step 10). Here, either the velocity limits could be set or all vectors of all the frames could be displayed in a scatter plot and discard the anomalies. In this study, the interest is with regard to the vectors representing the riprap. PIVlab performs a correlation-based analysis that identifies the average displacement of patterns of pixels within an interrogation window between successive frames. If one needs the precise center of each stone tracked individually, one may need to use an object-tracking algorithm. As the water flows from left to right in the flume, the riprap can move in the same direction and vertically up and down. However, it is very unlikely that the riprap moves in the opposite direction to the flow. However, there were vectors pointing in the opposite direction as there was turbulent flow around the riprap stones. Furthermore, as the water splashes on some of the stones it bounces instantaneously in the opposite direction and the software sometimes detects this and generates vectors. Using the scatter plot, all the vectors in the opposite direction to the flow were filtered out. This was followed by filtering based on image-based validation (Step 11). During this process there are two options, firstly filter low contrast and secondly filter bright objects. Then adjusting these threshold values one can suppress areas that have low contrast and vectors in bright areas. This is of significance because certain regions exhibited shadows, mainly due to the presence of columns supporting the glass panes. In contrast, other areas were exceptionally bright, owing to the reflection of floodlights off the glass panes. Additionally, there were areas characterized by the presence of white water, a consequence of turbulence created by water splashing onto the riprap stones. It is also possible to interpolate any of the missing vectors within each frame. However, it was decided not to use the interpolated values and rely solely on the calculated values as the reliability of interpolated values is somewhat degraded when compared to calculated values.

Subsequently, specific regions were delineated from which AVP values were obtained (Step 12). Four different areas were drawn: 1) downstream riprap only, 2) downstream riprap and filter, 3) crest riprap only, and 4) crest riprap and filter. Taking into consideration all the velocity vectors, PIVlab calculated the AVP for the selected area for each frame. These AVP values of each frame were saved in a text file (Step 13) which was later exported to MS Excel and MATLAB for analysis (Step 14). The velocimetry data pattern calculated using the PIV software is an average for the selected area. This area is static and does not follow the protective layer (i.e., the riprap stones) dynamically and only gives the reading of the movement that occurs within this selected area.

In order to establish where the breach initiates, a section of the riprap layer was selected from the top and bottom of the downstream slope section. The chosen segment measured approximately 20–32 cm in length. It was not possible to consistently select an area of the same length for every model. This was because in some cases the toe sections of the models were beyond the field of view of the cameras used (Table 2) or in other cases the columns holding the glass panes obstructed the visualization of the area of interest. Figure 5 presents a visual representation of the selected sections for model FP_01.

As illustrated in Figure 6, there were some limitations of this set-up as the camera failed to capture the most downstream toe section. These limitations have been summarized in Table 2.

The Table 3 provides a summary of the observations made on the breach initiation analysis and the initial purpose behind the model tests.

When looking into the results presented in Table 3 it can be seen that for all unsupported models with placed riprap, the breach initiation took place simultaneously in the top and bottom sections (Figure 7) whereas, for the unsupported models with dumped riprap, the initiation was observed in both the top and bottom sections, but not simultaneously (Figure 7). The breach initiation at the downstream toe indicates the point where the hydrodynamic forces surpass the static frictional force at the dam’s toe, between the riprap stones and the geotextile. In the case of unsupported models with placed riprap, the riprap layer experienced a sliding failure as noted in previous studies (Ravindra et al., 2020; Dornack, 2001; Siebel, 2007). It is important to recognize that the origin point on the horizontal axis does not correspond to the true zero point, which signifies the actual start of data collection (Figures 7, 8, 10, 11). This is the reason why it's referred to as “duration” instead of “time”.

Ongoing movements of the riprap and filter material are noted in the case of RP_01 (on the ramp) in comparison to the full dam FP_03 which includes the shell material. In the absence of this frictional resistance against the shell, the riprap and filter material on the ramp undergoes significant rearrangement. This was further established by looking at the AVP values for the protective layer before the main breach occurred Figure 8.

The data has been truncated to enhance graphical clarity and facilitate comparison. For instance, in Figure 8, the data has been adjusted to facilitate comparison, while Figure 7 provides a close-up view of the initiation points across all studied models. These initiation points denote the juncture where the graph’s overall trend begins to ascend, occurring approximately around the 6-s mark in Figure 8.

When looking into the time period preceding the main breach the time-averaged AVP value for FP_03 is 0.008 m/s whereas the time-averaged AVP value for the RP_01 model on a ramp without any shell material is 0.020 m/s. Figure 8 illustrates the AVP before the main breach where one can clearly see the higher values for the model on the ramp as compared to the full dams with the shell material. This may be due to the fact that there is less frictional resistance for the RP_1 model where the shell material is absent. Another reason could be the contrast in throughflow within the dam structure. Previous studies have shown that the dynamic throughflow within the dam structure can have a significant effect on the structural stability of a dam (Dezert et al., 2022a; Morán and Toledo, 2011). The flow within the shell material will increase the internal pore pressure. This in conjunction with the hydraulic drag and lift forces exacerbates the stability of the protective layer.

When examining Figure 9, which illustrates the data in Table 4, it becomes evident that there exists a linear relationship between the normalized peak AVP and the normalized inflow during the breach. Higher peak values are associated with higher inflows during the breach. However, one notable exception is readily illustrated in Figure 9. Model HPT_01 has the lowest peak values even though it has the highest breach inflow. The main reason behind this is that this peak AVP is an average for the selected section as highlighted in Section 2.2 (in this case the entire protective layer) and the material in the bottom section adjacent to the toe support remained stationary. V refers to the peak AVP value, while Vmax refers to the maximum peak AVP value for all the models presented. Correspondingly, Q refers to the inflow during the breach, while Qmax refers to the maximum inflow during the breach for all models presented.

Another aspect that was investigated was the impact of the toe support. This retaining element clearly had a significant influence on the movement in the most downstream sections of the riprap layer. Figure 10 offers a comparative visualization of the most downstream section of the riprap layer for two models with placed riprap, one with the toe support (HPT_01) and the other without (HP_01).

The physical modeling revealed that the placed riprap models exhibit significantly greater resistance as compared to dumped riprap models. Hence, able to withstand higher overtopping discharges. This is further illustrated in Table 4 and Figure 9. Moreover, the placed riprap models experienced a sliding failure while the dumped riprap models underwent progressive erosion. Additionally, for the slope and friction angle of the material examined in this study, no major movement was observed solely due to throughflow. This variation in construction philosophies, namely placed and dumped riprap also leads to a distinct variation in the behavior of the protective layer when it fails under overtopping conditions. Figure 11 visually presents the results obtained from the analysis of the entire protective layer.

As mentioned in the set-up and methodology section it should be noted that AVP values were obtained through a two-point calibration using the measured distance between two columns on the side of the flume. Although this study focuses on the AVP rather than the absolute velocity values, previous studies have shown that such calibration generates satisfactory results where the deviation between peak velocities from PIV analysis and smart stones ranged between 0.01 and 0.06 m/s (Ravindra et al., 2020).

The interlocking stones in placed riprap create a bearing structure and this structural layer seems to slide as a single entity in the absence of a rigid toe support. However, the unsupported models with dumped riprap lack such interlocking characteristics. Hence, these structures underwent progressive erosion where the riprap layer was removed in stages (Dezert et al., 2022a) where the precise location of breach initiation was inconclusive, as occurrences have been observed both at the top and bottom for the slope and friction angle examined in this study. Conversely, the results of the present study indicate that the supported model with placed riprap exhibited a different point of breach initiation i.e. at the top of the riprap layer. As the toe support provides stabilizing resistance the riprap layer is able to withstand higher hydrodynamic forces. While experiencing these forces, they are subjected to compaction, creating a gap at the transition between the crest and the downstream section which has been observed by Hiller et al. (2018). The results of this study exemplify this matter. Furthermore, past studies (Ravindra et al., 2021; Dezert and Sigtryggsdóttir, 2024; Dezert et al., 2022b) revealed that the deformation models with placed riprap and toe support underwent closely mimic the Euler buckling behavior of a long column which is pinned at one end and free at the other. The results of this study support this statement as the results reveal that the breach movement initiates at the top section (free end) of the protective layer while the movement of the bottom section is confined by the toe support. Figure 7 illustrates this graphically where there are insignificant velocity readings at the bottom section of the riprap layer as compared to the top which has significant velocity readings.

Continuous displacements of the riprap and filter material are revealed when it comes to RP_01 as compared to FP_03. As mentioned in [(Ravindra et al., 2020)] the sustained effect of the hydraulic drag and uplift forces causes the stones to rearrange continuously and create a more compact structure. Although this effect can be identified in models such as FP_01 and FP_03 as well, it is not as significant as RP_01. Just as there are frictional forces between the filter layer and riprap interface (Ravindra et al., 2020), there are also similar frictional elements between the filter and the shell material. This friction between the filter and the shell material seems to minimize the movement of the filter as well as the adjacent riprap layer. As RP_01 is built on a ramp (with geotextile), there is no shell material to provide the required and more realistic frictional resistance. When looking at Figure 9, it can be seen that RP_01 follows a similar evolutionary pattern to that of the other models. The reasoning behind this could be because the geotextile could distribute the loads more evenly thus mitigating its influence on the movement relative to the friction angle alone. The interaction between the filter and geotextile very much mimics the resistance that would arise naturally between the filter and shell interface. Nonetheless, this resistance is marginally inferior as it is not possible to perfectly replicate the micromechanical interlocking of the rock particles as illustrated in Figure 8. This brings forth the importance of including the shell for more realistic modeling.

The primary stabilizing agent or resisting force on the filter and riprap material is the frictional force (Wörman, 1993). This further expounds on the fact that the greater displacement over a unit time for RP_01 is due to the lower frictional resistance in the absence of shell material. Furthermore, the destabilizing drag force is significantly larger than the seepage force when the dam has been overtopped and the riprap (surface stones) experience only 50% of the seepage forces in comparison to the internal stones (Martin and Aral, 1971).

In the case of model HPT_01 the toe support restricts the movement of the protective layer significantly serving its purpose. Hence, even though the top section of the protective layer has been considerably displaced, the AVP for the entire section is quite low (Table 4) as the bottom section has undergone minor displacement.

This analysis also reveals the effectiveness of the supporting toe structure. When looking into the AVP values of the riprap layer at the top and bottom of the protective layer it can be seen that there is significant movement, particularly at the bottom section for structures with placed riprap without the supporting toe. In the case of placed riprap with a supporting toe, this movement at the bottom of the riprap layer is minor in comparison as can be observed in Figure 10.

The results revealed that models with placed riprap have a single velocity peak whereas models with dumped ripraps have double peaks in the velocity pattern as seen in Figure 11. These double peaks were a result of progressive surface erosion where the dumped riprap rolled down near the toe of the dam and piled up creating a restricting structure (like a small weir) backing up the water. When the flow is large enough, it breaks through this structure, creating a secondary peak as can be seen in Figure 12. In structures featuring placed riprap, the protective layer remains intact until the breach discharge reaches a significant level. This breach discharge value is typically higher than that for structures with dumped riprap (Dezert et al., 2022a). At that point, the entire structure, along with its protective layers, is removed in a single peak. This peak is typically higher than either of the double peaks seen in structures with dumped ripraps.