Numerical Investigations on the Instability of Boulders Impacted by Experimental Coastal Flows

Abstract

Coastal boulders transported inland by marine hazards, such as tsunamis and storms, are commonly found worldwide. Studies on the transport process of coastal boulders contribute to the understanding of a wide range of phenomena such as high-energy flow events, fluid-structure interaction, and coastal sediments. Consequently, it is crucial to understand how boulders move, but even more important to determine the instability condition for boulder transport. The hydrodynamic formulas including drag and lift coefficients are widely used to predict the incipient motion of boulders while few studies are conducted to evaluate the capability of these formulas. Recently, a series of laboratory experiments carried out at the Hydraulic Engineering Laboratory (Italian acronym LIDR) of the University of Bologna, Italy, revealed that boulders can start moving when the flow height and flow velocity are lower than the theoretical threshold computed by hydraulic formulas. In this paper, we use a numerical shallow water model to reproduce these freely available laboratory data with the aim of testing the capability of the model in capturing the main evolution of the process, and of casting new light on the instability condition of coastal boulders.

1. Introduction

Coastal boulders are regarded as the consequence of extreme marine events, and many imply tsunami events [1,2]. Thus, the term “tsunami boulder” was used as a substitution of coastal boulder, though other high-energy events, i.e., storms [3] and typhoons [4], are reported as triggers of boulder mobilization. The destabilization and transport of boulders (in mountain and coastal areas) in interaction with various geophysical flows, including the aforementioned marine hazards, is essential for understanding geophysical flows dynamics [5], for assessing natural hazards [6,7], and also for interpreting sedimentary evidence [8,9].

The definition of boulders comes from the classification of sediments where “particles” of a length ranging from 0.25 to 4.1 m [10] or even larger, are defined as boulders. According to sedimentological criteria, boulders can be classified as exotic and imbricated boulders [11]. The imbrication of boulders is also interpreted as the consequence of mass-transport processes, such as debris flows [9]. Therefore, the existence of boulders distributed in different zones is the possible sign of geophysical-flow activity [12]. From the perspective of the triggering mechanism, the initialization of boulders transport by geophysical flows has been rarely investigated, while the impact forces of various mass flows have been studied numerically [13,14] and experimentally [15,16], concentrating on the risk assessment and hazard mitigation.

The generic term of boulder is used here to denote a solid block that interacts with the water flow under different conditions. The conventional hydraulic formulas can be applied to analyze the incipient motion of a boulder. By considering the limit equilibrium state of a block acted by forces (i.e., weight, friction, drag, and lift forces) and torques, one can obtain the minimum flow conditions needed for the block destabilization [17,18]. The use of hydrodynamic formulas is straightforward while strong uncertainties on the analytical models are widely acknowledged. Numerical models [19,20] and parallel laboratory tests were implemented to study boulder transport [21]. The essence of these numerical models was to describe the flow-boulder interaction by coupling a shallow water model for flow dynamics and an analytical model for the block.

Although gradually more sophisticated studies have been carried out to figure out the dynamic process of boulder transport, investigations are still needed to verify the accuracy and capability of the existing numerical models and to better understand the basic physical process. On the basis of the experiments on the incipient motion of boulders impacted by water flows conducted by Bressan et al. [22], the aim of this paper is to explore the performance of the shallow water equations in predicting the flow properties (depth and velocity) and to improve our understanding of boulders mobilization. The investigation focuses on the properties of the impacting flow with regard to the instability condition of the boulders, and the numerical results are compared with the laboratory experiments. The numerical model is first calibrated against the experimental data of pure flow and then applied to the data of a water flow impacting a boulder. We anticipated that we will find that the classical shallow water model is capable of capturing the main features of the experimental flows while more complicated models should be developed to improve the accuracy especially in the simulation of the flow velocity. Considerations on the capabilities of our numerical model to predict block mobilization are done in the last section of the paper.

2. Laboratory Experimental Setup

In this paper we refer to the experiments carried out at the Hydraulic Engineering Laboratory of the University of Bologna, Italy [22], whose data are freely available from the site (ftp://137.204.48.34). In those experiments, simple blocks of different geometries, densities, and shapes located on a homogeneous and gentle slope plane were impacted by various flows that were generated by quickly opening a water tank (2.0 m long and 0.5 m wide), i.e., by using the same flow-source technique as Nandasena and Tanaka [21]. The bore was running along the flat part of the flume and then climbing up a dry slope where it impacted the block. In this way, the experiments simulated a tsunami bore hitting a coastal boulder. The whole flume was 11 m long and 0.5 m wide. A sketch of the laboratory setup is displayed in Figure 1.

In contrast to Nandasena and Tanaka [21], who set the initial water heights in the tank at a high value to reduce the effects of turbulence in the flow along the slope, Bressan et al. [22] increased the initial water height gradually in the basin starting from a small value with the aim to better identify the instability condition of the block. Notice that the flume remained wet when the experiments were conducted. This last feature has been taken into account in the numerical simulations illustrated in this paper by considering two different water heights in the initial configuration, i.e., by including the initial water level in the tank h₀, which is the chief element to determine the energy of the flow, as well as the height of the thin residual water layer h₁ in the flat section of the flume. Along the slope, a dry bed condition was assumed.

In all the experiments, Bressan et al. [22] placed a set of instruments consisting of three water level sensors (WL), three Doppler velocimeters (DVn), and one 100 fps camera to record the information. DV1 was set up with the acoustic beam aligned with the flume basis and was suitable to record the along-flume component of the flow velocity. The velocity of the flow impacting the block was measured by combining DV2 and DV3 data, installed 2 cm ahead of the block. The time when the block started moving and the corresponding flow height were determined by analyzing the images recorded by the camera. All relevant details can be found in Bressan et al. [22].

3. Numerical Model

3.1. Governing Equations

By means of a depth-average technique, the equations of mass and momentum balance for the fluid are reduced to the shallow water equations (SWEs), which have been widely applied to problems of geophysics and environmental engineering, such as the dynamic evolution of flow-like avalanches [23], the dam break problem [24], and the propagation of tsunamis [25]. Moreover, a laboratory dam-break experiment with water flowing over a triangle hump was set up by the European project EU CADAM (European Union Concerted Action on Dam Break Modeling) to provide a benchmark that was later widely used to validate the classical SWEs with the additional hydraulic friction term [26]. The governing SWEs for the experimental flows in a 2D channel can be given the form:

$${\partial }_{t}h+{\partial }_x\left(hu\right)+{\partial }_y\left(hv\right)=0$$

(1)

$${\partial }_t\left(hu\right)+{\partial }_x\left(h{u}^2+\frac{gcos{\theta }_x{h}^2}{2}\right)+{\partial }_y\left(huv\right)=hgsin{\theta }_x-\gamma hu-hgcos{\theta }_x{\partial }_{x}b$$

(2)

$${\partial }_t\left(hv\right)+{\partial }_x\left(huv\right)+{\partial }_y\left(h{v}^2+\frac{gcos{\theta }_y{h}^2}{2}\right)=hgsin{\theta }_y-\gamma hv-hgcos{\theta }_y{\partial }_{y}b$$

(3)

where h is the depth-averaged flow depth, b(x,y) is the basal topography, u and v are the depth-averaged velocities along the x and y directions, θ_x and θ_y are the local slope angles of the channel along the x and y directions. The frictional terms in the equations are γhu and γhv, where $\gamma =\frac{g{n}^2\sqrt{{\left(hu\right)}^2+{\left(hv\right)}^2}}{h^{7/3}}$ and n is the Manning’s coefficient. Typically, the range of n is from 0.01 to 0.06 s m^−1/3, depending on the roughness of the bottom [25].

In a 1D channel the above equations simplify to:

$${\partial }_{t}h+{\partial }_x\left(hu\right)=0$$

(4)

$${\partial }_t\left(hu\right)+{\partial }_x\left(h{u}^2+\frac{gcos{\theta }_x{h}^2}{2}\right)=hgsin{\theta }_x-\gamma hu-hgcos{\theta }_x{\partial }_{x}b$$

(5)

3.2. Numerical Algorithm

The numerical scheme we adopt here to solve Equations (4) and (5) is the central Nessyahu–Tadmor (NT) scheme [27], combining a first-order Lax–Friedrichs scheme with a piecewise linear reconstruction. The staggered cell average values are computed at the interfacing break-points. The method has the advantage of the simplicity given by a Riemann-solver-free approach and has been proven to be adequate in solving shallow water models [28,29]. The adopted scheme solves the system of equations on a fixed Cartesian grid, identifying shocks as those regions where large space gradients of the variables are found. In the numerical simulations, the gentle-slope is treated as a slope with a negative angle –$\theta$ and the block is implemented as a basal topography feature. The implementation of the numerical algorithm presented in this paper has been verified against analytical solutions, experimental data, and a real landslide case [30].

4. Hydraulic Formulas for Instability Conditions

The theoretical instability conditions for a block can be analytically derived from the expressions of the forces and torques acting on the block. The block starts moving as soon as the balance condition either of forces or torques is broken. The laboratory experiment setup can be naturally simplified through the sketch shown in Figure 2.

Following the analysis made by Bressan et al. [22], the instability inequality for the torque can be expressed as:

$$D{b}_D+L{b}_L=0.5{\rho }_w{v}^2\left(b_D{C}_D{A}_{wfs}+b_L{C}_L{A}_{wbs}\right)>b_{W}W$$

(6)

where L and D denote the lift and drag forces respectively. Further, W is the weight expressed as $W=gV\left({\rho }_s-{\rho }_w{V}_w/V\right)$, ${\rho }_w$ is the water density, C_L is the lift coefficient, A_wbs is the area of the wet bottom surface, A_wfs is the area of the wet front surface, $v$ is the depth-averaged water velocity, C_D is the drag coefficient, V is the boulder volume, V_w is the wet volume, and θ is the slope angle. With b_L, b_D, and b_W representing the force arms of the lift and drag forces and of the effective weight respectively, the minimum flow velocity needed to break the torque equilibrium is given by:

$$v^2>\frac{2b_{W}W}{{\rho }_w\left(b_D{C}_D{A}_{wfs}+b_L{C}_L{A}_{wbs}\right)}=v_M^2$$

(7)

To deduce the instability inequality for the forces, let us introduce the following dynamic pressures and their corresponding velocities:

$$P_k=\frac{{\rho }_w}{2}{v}_k^2=\frac{W_p+\mu W_n}{C_D{A}_{wfs}+\mu C_L{A}_{wbs}}$$

(8)

$$P_D=\frac{{\rho }_w}{2}{v}_D^2=\frac{W_p}{C_D{A}_{wfs}}$$

(9)

$$P_L=\frac{{\rho }_w}{2}{v}_L^2=\frac{W_n}{C_L{A}_{wbs}}$$

(10)

$$P_k-P_D=\mu W\frac{C_D{A}_{wfs}cos\theta -C_L{A}_{wbs}sin\theta }{C_D{A}_{wfs}\left(C_D{A}_{wfs}+\mu C_L{A}_{wbs}\right)},P_L-P_k=W\frac{C_D{A}_{wfs}cos\theta -C_L{A}_{wbs}sin\theta }{C_D{A}_{wbs}\left(C_D{A}_{wfs}+\mu C_L{A}_{wbs}\right)}$$

(11)

Here $W_p=Wsin\theta$, $W_n=Wcos\theta$ and $\mu$ is the static friction coefficient. With the experimental configuration $tan\theta <\frac{C_D{A}_{wfs}}{C_L{A}_{wbs}}$, we can obtain the inequality chain of velocities v_D² < v_k² < v_L² [22]. The force equilibrium is broken if driving forces are higher than resistance forces. The friction force has the expression:

$$F_s=max\left(\mu \left(W_n-L\right),0\right)$$

(12)

When $F_s\ne 0$, the block moves if the drag force exceeds the sum of the along-slope component of the weight and of $F_s$, i.e.,

$$D>W_p+\mu \left(W_n-L\right)$$

(13)

from which one gets:

$$v^2>\frac{2\left(W_p+\mu W_n\right)}{{\rho }_w(\mu C_L{A}_{wbs}+C_D{A}_{wfs})}=v_k^2$$

(14)

Eventually, the instability condition for the flow velocity can be obtained by combining the expressions (7) and (14), i.e.,

$$v^2>min\left(v_M^2,v_k^2\right)$$

(15)

If boulders are rectangular-based prisms, one can easily obtain the areas of the front wet bottom surface $A_{wbs}$ and of the wet front surface $A_{wfs}$ in case of full and partial submersion. In this paper we will focus on a specific set of tests by Bressan et al. [22], namely those involving a block of 2.95 × 3.10 × 3.00 cm (length, width, and height). Considering the experiment setting displayed in Figure 1, the corresponding theoretical instability profile is given in Figure 3. The adopted densities of the water flow and of the block are 1000 kg/m³ and 1900 kg/m³ respectively. The lift coefficient is set to C_L = 0.178 as given by Nandasena and Tanaka [21]. As for the value of the drag coefficient, we note that the choice of C_D is empirical and recent studies use a value close to 2 [21,31,32]. However, we notice further that the inverse analysis of impact pressure of snow avalanches suggests that C_D can range between 2 and 6 [33]. Here, based on the Bressan et al. [22] study, where it was shown that adopting values of C_D between 0.8 and 2.0 underestimates the mobilization capability of the experimental flow, we present theoretical instability curves with two values of the drag coefficient, that are C_D = 2 and C_D = 6. The friction coefficient $\mu$ of the block is 0.5 measured in the laboratory. From Figure 3, one can observe that v_k is always smaller than v_M, if the curves are computed with the same C_D value. Therefore, all theoretical instability graphs in the following sections will only refer to v_k.

5. Numerical Investigations

The model formed by the Equations (4) and (5) is the classical 1D shallow water model with the additional Manning’s friction term, which is solved through the method given in Section 3.

In the laboratory tests, the bore flows along the flume and then climbs up the dry slope hitting the block located there. The flat sector and slope beds are made of different materials and thus, in principle, are characterized by different friction coefficients [26,34]. Based on a number of numerical back-analyses, we found that for the flat sector of the flume the value n = 0.0125 s m^−1/3 provides the best fit with the observed flow. As for the incline, we found it more convenient to use values of the Manning’s coefficient ranging from n = 0.0125 s m^−1/3 to n = 0.0225 s m^−1/3.

We compared the laboratory data with numerical results for two different series of experiments, i.e., in tests where no block is placed in the slope, and in tests where the measured flow impacts a block. Each test refers to initial water levels in the tank covering the range h₀ = 0.15–0.28 m. In our simulations, we assume that there is an initial residual water level h₁ also in the flat part of the flume, according to experimental records.

The simulations were conducted with a space grid size set to Δx =0.01 m to treat pure water flows with no block, and to Δx = 0.005 m to handle water flows impacting a block. In the latter case, the block was simplified as a fixed topography, i.e., a relief, and the fix boundary condition was enforced at the node representing the block when the height of the impacting flow was lower than the height of the block.

5.1. Flows without Block

To illustrate our results, it is convenient to compare data and simulations separately for the phase when the bore advances along the flat flume and the phase when it climbs up the slope. In the first part, we present the propagation times of the flow deduced from the records of the bore transit at the water level gauges that are installed at the beginning and in the middle of the flat section of the flume (see Figure 1). In the second part, we consider data taken along the flume and just in front of the block. More precisely, we make use of the flow velocities determined in positions where the Doppler velocimeters are installed: DV1 makes measurement in position P1, while DV2 and DV3 are used together to obtain experimental values in position P2. In addition, the camera images provide data on the flow height in position P2 (see Figure 1).

5.1.1. Bore Propagation Times Measured by the WL Records

The time arrivals of the bore at the WL sensors can be inferred from the recorded time series. Notice that measuring the time elapsed between the exact moment of the dam-break and the bore arrival at the WL gauges may be affected by a bias due to uncertainty on the exact time of the manual opening of the gate. It is therefore preferable to trace the transit times from one WL to the others. If we denote by t₁, t₂, and t₃ the instants when the bore passes by the gauges WL1, WL2, and WL3, the time intervals Δt₁ = t₂ − t₁ and Δt₂ = t₃ − t₁ represent the propagation times of the bore and can be used to compare experimental and numerical data, which is done in the plots of Figure 4.

From the graphs, one can observe that the experimental values for Δt₁ and Δt₂ referring to the same level of water in the tank show a significant spread (in the order of 0.1–0.2 s) from one test to another, which is likely the effect of random errors and also of different levels h₁ of the residual water left on the flume. The numerical simulations mainly capture both the spread and the decreasing trend of the propagation time with increasing h₀, which is indicative of an increasing bore velocity. It is worth noting that the discrepancy between the observed and computed propagation times is unsatisfactory in the upper graph of Figure 4. The bad performance of the numerical model in computing Δt₁, which is the propagation time taken very close to the gate, may derive from the “long-wave” approximation that constitutes the fundament of the SWEs. Looking at the numerical values for Δt₁, the computed bore moves slower than the experimental bore when h₀ is small and vice-versa for larger h₀. The data presented here refer to simulations with Manning’s friction coefficient equal to 0.0125 s m^−1/3 but using other values does not help improve the fit. On the other hand, it is very important to outline the very good agreement between data and calculations shown in the plot of Δt₂, which means that the numerical code provides reliable values of the flow front velocity at the WL3, and therefore that WL3 is sufficiently far from the source region to allow a safe application of the SW approximation.

5.1.2. Flow Height at Position P2

In this section, we present experimental and numerical time series of the flow height at the position P2, with time origin taken at the flow arrival time. More precisely, we provide graphs representing how the flow height changes in P2 during the whole process where the flow first runs along the slope up to a maximum distance and then comes back as a reverse flow. The simulations are carried out within the limits of friction coefficients n_min = 0.0125 s m^−1/3 and n_max =0.0225 s m^−1/3 for the slope and with various h₁ values: h₁ = 1, 2, 3, 4, 5, 6, 7, 8 mm. For a better representation of the results, instead of plotting the curves corresponding to all possible combinations of the parameters n and h₁, we consider the two areas swept by all curves sharing the minimum (n_min ) and the maximum (n_max) value of the Manning’s coefficient and denote them by HN_m and HN_M, respectively. In Figure 5, we plot distinctly the intersection of such areas HNI (HN_m ∩ HN_M) and the relative complements HN_mM (HN_m − HN_M) and HN_Mm (HN_M − HN_m).

The time interval of the process is captured quite well by the numerical model. Based on the width of the two curves, it can be seen that the heights of h₁ has some influence on the flow height when the initial basin water level h₀ is low, and that this influence gradually decreases. Comparing these curves with experimental data, one may observe that the model underestimates the flow height when the initial basin height is small and, vice-versa, it overestimates it when the basin level is high. As for the influence of n, considering that the overall trend of HN_mM agrees better with the data (dashed lines) in the upslope and downslope part of the curves, one may conclude that the choice of n = 0.0125 s m^−1/3 is better in capturing the height evolution profiles at P2.

5.1.3. Peak Velocity at P1 and P2

The flow formed by the gate opening moves with the maximum velocity at the bore front, which is a region affected by high turbulence (see Bressan et al. [22]). Such a region is hard to catch properly in the calculations since the velocity appears as a step function. The peak velocities taken at the DV1 (position P1) and along the slope (position P2) are compared with the observed data under various numerical configurations in Figure 6a,b, respectively. To reduce turbulence noise, the experimental data have been smoothed by a cubic spline approximation. Three different values of friction coefficients have been considered for the slope (n₁ = 0.0125 s m^−1/3, n₂ = 0.0175 s m^−1/3, and n₃ = 0.0225 s m^−1/3) and for each friction coefficient, different h₁ values have been implemented, which allows us to obtain peak velocity mean values with associated error bars, as shown in Figure 6.

By looking at the error bars, one can deduce that peak velocities are sensitive to the variation of h₁, no matter the change of the basin water level. For the peak velocity at P1, the numerical model generally overestimates the flow velocity. Instead, the model provides reasonable peak velocity at P2 when the initial basin water level is higher than 0.24 m.

5.1.4. Froude Number at Position P2

A flow and its impact power can be characterized by the Froude number as was underlined by the observations of snow avalanches by Sovilla et al. [33]. Here, the Froude number Fr, defined as the ratio v/(gh)^1/2 (v, h, and g being the flow velocity, height, and the gravity acceleration), is used to compare numerical and experimental data. The heights and velocities of the arrival flows are computed at position P2. Curves of Froude numbers vs. time are given in Figure 7. It may be seen that in all cases, Fr is a monotonic decreasing function of time, with the maximum value at the initial time, and that it passes from supercritical to subcritical values. This is due to the fact that when the bore arrives at the measuring point, its front has the minimum height and the maximum velocity, and therefore the corresponding Fr value is very high. In general, our numerical model overestimates the Froude number when h₀ is smaller than 0.2 m, while it locates it at a reasonable range for higher values of h₀. Notice that hereafter we will disregard the bore front region where the flow height is extremely low, i.e., smaller than 0.003 m since, as already noted in the previous section, experimental and numerical velocity data may be unreliable here.

5.2. Flows with Block

To model a condition where the water flow interacts with a block, we represent the block as a topographical relief and the impacting flow property is investigated by numerically monitoring the grid node placed in front of the block at a distance of 2 cm. The velocity is set to zero in all the nodes along the block surface. Considering the experiments of the quasi-cube block with dimensions 2.95 × 3.10 × 3.00 cm (length, width, height), the block topography is simplified as a cubic hump, 2.95 cm long and 3.00 cm high.

The time intervals between the instant of the flow arrival at the block and the instant of the instability onset are extracted from the experimental data and are shown in

Table 1. Time interval (in seconds) between the flow arrival time and the time of block instability.

h0 (cm)	15	16	17	18	19	20
Case 1	stable	0.51	stable	0.17	0.19	0.12
Case 2	stable	stable	0.34	0.20	0.15	0.16
Case 3	stable	1.02	0.21	0.18	0.13	0.09
Case 4		stable	0.30	0.22	0.11	0.08
Case 5		stable	0.37	0.18	0.10	0.12
Case 6		stable	0.22	0.16	0.17	0.11
Case 7		stable	0.22	0.25	0.13	0.14

. There are three experimental cases for the configuration of h₀ = 15 cm and seven cases for the other configurations. Considering the data, one can distinguish three different dynamic conditions, i.e., block stability (h₀ = 15 cm), block probable stability (h₀ = 16 and 17 cm), and block instability (h₀ > 17 cm).

Based on these data and on data plotted in the instability graph that will be shown later, one can conclude that the present numerical model combined with the hydraulic formula related to a higher drag coefficient (C_D = 6) might provide reasonable results to predict the minimum flow conditions for the block mobilization. A further relevant consideration is that the block destabilization, if occurs, takes place in the very first moments of the flow attack, with the delay interval getting shorter and shorter as the initial level of the water in the tank increases. Mobilization occurs much before the instant when the flow reaches the peak height (which occurs around 2 s after the flow arrival as shown in Figure 5). This suggests that to provide very accurate predictions on block stability or instability any numerical model needs to reproduce the time series of the bore transit especially in the very beginning, that is soon after the arrival of the bore front.

5.2.1. Flow Height Profiles

The experimental settings cover the range of the experimental flow from the situation where the flow is unable to mobilize the block to the situation where the flow breaks the block balance and moves it. In this section we compare numerical and experimental data under different configurations. The flow information is extracted from the node placed 2 cm ahead of the block in the numerical simulations.

From Figure 8, it can be seen that the introduction of the block on the flow path modifies the flow height in front of the block, especially when the flow is not strong (namely for smaller values of h₀). We note further that the numerical model overestimates the flow height in the initial (0.5 s long) phase mainly because of the existence of the fix 1D boundary that represents the block. Interestingly, the experimental results computed without (NWB) and with the block (NWTB) are shown between the curves. Also noticeable is that, when h₀ = 0.17 m and h₀ = 0.18 m, NWB is capable of capturing the evolution of the flow height quite well (see Figure 8c,d). This ability is proven also for higher values of h₀, whose graphs are however omitted here. These results show that the influence of the block on the flow height is negligible when the initial basin level is sufficiently large to mobilize the block soon.

5.2.2. Flow Conditions in Front of the Block

Experimental and numerical data of the flow measured in P2, that is the position 2 cm back from where the block will be placed, are displayed on the theoretical instability graph shown in Figure 9. For each single case, the diagram is built by plotting the values of the normalized flow velocity square vs. the normalized flow height only for the various instants during the bore advancement. The corresponding data related to the backflow are omitted. Each case is then represented by a set of points (numerical and experimental) corresponding to different times. In a sense, the time can be thought to evolve from left to right in the diagram while the flow height is increasing. The points in the leftmost region of the graph (smallest flow height) refer therefore to the initial stage of the flow, whereas the points with the largest flow height refer to the central stage of the flow (usually around 2 s, see Figure 5) where the flow attains the highest peak value. Considering the nearly-monotonic shape of all curves, that are almost flat at the beginning and then exhibit an increasing negative gradient, one deduces that in general (i) the impacting flow arrives with high velocity and low height, (ii) the flow velocity keeps nearly stable while the flow height gradually increases and then (iii) the flow decelerates rapidly with the end of the run-up process.

In the diagram we compare experimental and numerical results without a block and with a block. For the shown numerical runs the Manning’s coefficient is chosen as n₁ = 0.0125 s m^−1/3. Using data of Table 1, we are able to identify the flow properties exactly at the block destabilization instant. These are shown as INSTB dots (instability dots) in Figure 9.

Considering the numerical curves, we see that the cases with and without block differ only slightly for all values of h₀, the discrepancy being higher for lower values. Even the experimental data show a difference between the block-yes and block-no cases that tend to decrease with the increase of h₀, but the discrepancy is much larger especially for the smallest values of h₀ (see e.g., graphs (a) and (b) of Figure 9). By comparing simulations and laboratory data, one notes that numerical curves are always external and embrace the observations, which means that they overestimate the flow velocity in the leftmost section of the graph and overestimate the flow height when velocities are small.

Looking at the instability graph, one reaches the conclusion that our results predict that the numerical flow will make the block unstable when h₀ is higher than 0.18 m (C_D = 2), or higher than 0.16 m (C_D = 6). On the other hand, experimental data show that the onset of the block instability (marked with the solid dots INSTB in Figure 9) occurs in several cases within the stability region of the diagram. This is especially true for h₀ = 0.16 m and h₀ = 0.17 m, while for higher values of h₀ the destabilization takes place at the boundary of the instability domain.

It is further worth noting that the block is mobilized when the flow height approximates the threshold of 4–5/10 of the boulder height (except for the data in h₀ = 0.16 m), observation that has been also reported by Bressan et al. [22]. This may indicate that the flow height can significantly contribute to some hydraulic forces when it reaches such a limit.

6. Discussion and Conclusions

The processes of boulders transport are important to study sedimentology, hazard assessment, and geophysical flow evolution. On the basis of freely available data by Bressan et al. [22] of the incipient motion of a block acted by experimentally generated flows, a 1D shallow water numerical model was developed and used to explore the model capability and to improve our understanding of the processes. We analyzed experimental data that refer to two conditions, namely of bores propagating on a slope and bores impacting a block placed on that slope. The experimental flows were generated by opening a gate of a water tank with initial water level in the tank sweeping a range sufficiently large to cover a reasonable interval of flow velocities and heights. The numerical model was found to reproduce quite well the propagation time of the bore in the flat portion of the flume, with discrepancies only found in the flume sector closer to the gate, that is in the region of bore formation, which however do not affect the analysis. In the slope, we found that the numerical model reproduces well the overall time series of the bore height from the first arrival to the final backflow phase. However, the numerical results show a slight overestimation of the bore height, larger for more energetic flows, that is for the ones generated by higher initial water levels. This overestimation persists even by changing the bottom friction coefficients along the slope in the model within reasonable limits. The reproduction of the height of flows impacting the block allows interesting considerations. To notice first is that both experimental and numerical data show that the influence of the block on the flow height evolution in the position just in front of the block is relevant for less energetic flows and tends to decrease when flows are generated with higher initial water levels. The explanation is simple, since more energetic flows are able to destabilize the block soon after hitting it, and after block mobilization they behave like in cases without block. Interestingly, when the flow is really perturbed by the presence of the block, i.e., in the first stage of the block attack for lower values of h₀, empirical heights lie between heights simulated with block and without block. In a sense, it shows that the perturbation induced by the numerical block is overestimated, which is very likely because of the use of a 1D model, where the block is represented by a fix topography. Indeed, this implies that the flow is bound either to stop against the block or to flow above, but it is not allowed to go around the obstacle like it instead happens in the laboratory experiments.

As for the ability of the model to reproduce the experimental flow velocity on the slope, we found that results are less satisfactory than for the flow height. This is due to (1) the large spread of experimental data affected by high-frequency turbulence that imposes smoothing to get manageable velocity curves of the bore front, and also to (2) the fact that, as the simulated results show, the velocity profile along the bore is the highest in the front region and decreases gradually when moving away from the front. The consequence is that, like in a shock wave, the grid nodes in the neighborhood of the front pass from being dry (with zero velocity) to being wet (with high velocity) in a very short time. Though this situation is managed by the code and is no cause of numerical instabilities, the accuracy of the results may be less than desired.

With the above considerations in mind, let us focus on what we learned from the application of our numerical model to the laboratory data. The instability diagrams of the normalized flow square speed vs. the normalized flow height are useful since they allow one to compare simulated and experimental data locating them in the instability/stability domains identified by theoretical hydraulic formulas. When one puts all data related to a single case flow, that is to a single bore passage, on that graph, one sees that these data are arranged along a curve, which has a typical shape with a monotonic trend. It can be seen as a trajectory for a representative point of the flow in the phase space ($v_n^2,h_n$). If we restrict to the phase of bore advancement, such point moves rightward starting from the left end of the trajectory, because when the bore front attacks the block its velocity is high whereas its height is very small. On the bi-log graphs these trajectories start almost horizontal from the vertical axis (which corresponds to a stage of almost constant bore speed and increasing bore height), then exhibit an increasing inclination until the horizontal axis is almost normal (which corresponds to the stage of maximum inundation height and vanishing velocity). The two end points, that is the intersections of the trajectory with the diagram axes, can be used to identify the trajectory. These kind of trajectories are associated with the way the bore is formed, that is through the dam-braking mechanism realized by manually releasing a gate. Increasing the level of water in the tank has the effect of producing more energetic bores that on the graph follow trajectories having end points more distant from the graph origin, so that one can establish a one-to-one relation between the level of water in the tank h₀ and the pair of end points. All trajectories start on the stable part of the domain. If h₀ is large enough and correspondingly the trajectory left end point is high enough, the trajectory may intersect the theoretical stability curve, which implies that the flow is able to mobilize the block. On the other hand, if h₀ is too small the trajectory does not cross the stability curve at all, and the flow remains always in the stability region, being theoretically unable to displace the block. With this understanding, it is clear that the generation mechanism can produce flow trajectories that sweep a relevant area of the stability graphs but cannot enable one to fully explore the stability region. In particular, the stability domain characterized by low flow velocity and large flow height could be unreachable by such flow trajectories. To cover that region, the only strategy will be to increase the level h₀, but this way the trajectory will start from an end point on the velocity axis high enough to enter the instability region, and the block will be mobilized. Therefore, to produce slow and high flow, one should consider other generation mechanisms, such as fluid circulation produced by hydraulics pumps where velocity and flow height are not correlated; one possibility could be to remove the water from the flume at the upper end of the slope and to re-input it on the other side of the flume at the beginning of the flat bottom sector, creating therefore a water circulation with controlled fluxes and volumes.

As for the block stability, we can highlight two relevant observations. First, we note that numerical trajectories computed with and without the block do not differ too much in the phase space of the system but are systematically above the set of points that describe the experimental trajectories. Second, looking at the experimental points in the phase space marking the onset of block mobilization, we may point out that all of them are well below the theoretical instability curves. This latter feature was already noted by Bressan et al. [22]. Their conclusions were that probably the theoretical curves are based on an over-simplified hydraulic model (not considering flow height that do not completely submerge the block), but more importantly that their experimental data were affected by high-frequency fluctuations due to turbulence, so that instability could have been caused by temporary velocity peaks, representing large deviations from the mean flow velocity. If one were able to plot these deviations in the stability graphs, the experimental points will appear shifted upward, that is closer to or even within the theoretical block instability region. Our numerical results confirm the validity of those conclusions. Indeed, the shallow water model cannot account for the flow turbulence, and in case of turbulent flow is only able to compute mean flow properties, but with a caveat. Turbulence is known to subtract energy to the mean flow since the total available energy should feed both features. On the other hand, the only way to subtract energy to the numerical mean flow is by means of bottom stresses. Our results show that changing the bottom friction in the flume does not modify significantly the computed flow, nor, as a consequence, the flow trajectory in the phase space. Therefore, bottom friction cannot be used to reduce the mean flow energy in a shallow water model as a substitute of a turbulence dissipation mechanism. This explains why, considering cases with the same initial water potential energy (i.e., cases with the same initial water level in the tank), numerical trajectories lie always above the corresponding experimental points. Further, it is worth to highlight, that, though with the limitations mentioned above, if we group results according to the initial water level, and we assess the instability from the intersection of the numerical with the theoretical curves, the predicted block instability respects almost fully the experimental data if we use the curve obtained with C_D = 6, that is block stability for h₀ = 15 cm, near instability for h₀ = 16 cm and block mobilization for h₀ larger than 16 cm.

As concluding remarks, we can state that combining laboratory experiments and the numerical analyses provides a very powerful strategy to interpret physical phenomena, and in our case to improve our understanding of the mobilization of boulders by hydraulic flows. In our study we found that simulations and laboratory data confirm that the destabilization of boulders, occurring when the mean flow dynamics in principle is unable to cause it, can be explained by the relevant role of temporary fluctuations of the flow velocity, that are able to move the representative point of the flow in the phase space of the system upward, from the block stability into the block instability region. It is relevant to add that our conclusions are based on the analysis of data of only one of the types of blocks taken into account by Bressan et al. [22], and further, that, due to the limitation of the numerical model, the instability of blocks is interpreted only in terms of the 1D interaction between the flow and the block. The complex interaction mechanisms taking place in a 3D system will be studied with more advanced numerical techniques in our further research. However, in spite of the formal limitations, we think that our present findings might have a general validity and can also be applied to many other block shapes.